Dynamical excision boundaries in spectral evolutions of binary black hole spacetimes

We represent a general time-dependent map parameter λ(t) as a polynomial in time with a piecewise constant Nth derivative:

$$\begin{equation} \lambda (t) = \sum _{n=0}^{N} \frac{1}{n!}(t-t_i)^n\lambda ^n_i, \qquad {\rm for}\quad t_i \le t < t_{i+1}, \end{equation} \tag{ 11 }$$

where for each time interval t_i ⩽ t < t_{i + 1} the quantities $\lambda ^n_i$ are constants.

At the beginning of each new time interval t_i, we set the constants $\lambda ^n_i$ in equation (11) as follows. First for n = N we set $\lambda ^N_i=U(t_i)$, where U(t) is the control signal defined in detail below. For n < N we set $\lambda ^n_i= {\rm d}^n\lambda /{\rm d}t^n|_{t=t_i}$, where the derivative is evaluated at the end of the previous time interval. In this way all the derivatives of λ(t) except the Nth derivative are continuous across intervals. The goal will be to compute the control signal U(t) so as to drive the map parameter λ(t) to some desired behavior. Before we describe how to compute the control signal U(t), we first discuss the control error Q, which will be used in the computation of U(t).

To appropriately define the control error Q, one must answer the question of how a small change in the map parameter corresponds to a change in the observed variables. If the control error is defined to be too large, then the controller will consistently overshoot its target, potentially leading to unstable behavior; conversely, if the control error is defined to be too small, then the controller may never be able to reach its target value.

If there exists a target value of λ(t), call it λ_target, that does not depend on the map but may depend on other observable quantities in the system (call them A, B,...), then we define

$$\begin{equation} Q = \lambda _{\rm target}(A,B,\ldots ) - \lambda . \end{equation} \tag{ 12 }$$

The goal is to drive Q to zero and thereby drive λ to λ_target.

If instead, as often happens for nonlinear systems, the target value of λ depends on λ itself, even indirectly, then we define Q differently using a generalization of equation (12): we require that λ attains its desired value when Q → 0, and we require that

$$\begin{equation} \frac{\partial Q}{\partial \lambda } = -1 + \mathcal{O}(Q). \end{equation} \tag{ 13 }$$

The primary motivation for this condition is its anticipated use in relating time derivatives of Q to those of λ (see equation (18), below). Note that in either equations (12) or (13), Q could in principle be multiplied by an arbitrary factor; however, if this were done, that factor would need to be taken into account in the computation of the control signal U(t) below. Without loss of generality we assume no additional scaling.

In the case of several map parameters λ_a(t) with corresponding Q_a, where a is an index that labels the map parameters, the desired value of some λ_a may depend on a different map parameter λ_b. In this case, we generalize equation (13) and require that each Q_a satisfy

$$\begin{equation} \frac{\partial {Q_a}}{\partial {\lambda _b}} = - \delta _{ab} + \mathcal{O}(Q), \end{equation} \tag{ 14 }$$

where δ_ab is a Kronecker delta. This criterion ensures that we can treat each λ_a independently when all control errors are small. A way of understanding equation (14) is to consider a set of $Q^{\prime }_a$ that are obtained via equation (13) without regard to coupling between different λ_a. Then a set of Q_a satisfying equation (14) can be obtained by diagonalizing the matrix $\partial {Q^{\prime }_a}/\partial {\lambda _b}$. In the remainder of this section we assume that if there are multiple map parameters λ_a, the corresponding control errors Q_a satisfy equation (14). We therefore drop the a indices and write equations for U(t) in terms of a single Q(t) satisfying equation (13) that represents the control error for a single λ(t).

The control error Q(t) is a function of several observables. In the case of the mapping functions that are designed to move and distort the excision boundaries to follow the motion and shapes of the apparent horizons, Q(t) is some function of the current position or shape of one or more apparent horizons. The precise definition of Q(t) is different for each map parameter, and depends on the details of how each map parameter couples to the observables. We will discuss the control error Q for each of the map parameters in section 4. But it is not necessary to know the exact form of the control error in order to compute the control signal U(t); it suffices to know only that U(t_i) is equal to $\lambda ^N_i$ in equation (11), and that the control error Q obeys equation (13).

We now turn to the computation of the control signal U(t). There is some flexibility in the control law determining U(t), so long as key feedback mechanisms are in place (as shown in section 2). In SpEC, we use either a standard PID controller,

$$\begin{equation} U(t) = a_0\int Q(t)\,{\rm d}t + a_1\,Q(t) + a_2\frac{{\rm d}{Q}}{{\rm d}t}, \end{equation} \tag{ 15 }$$

or a special PD (proportional-derivative) controller,

$$\begin{equation} U(t) = \sum _{k=0}^K a_k\frac{{\rm d}^{k}{Q}}{{\rm d}t^{k}}, \end{equation} \tag{ 16 }$$

where typically K = 2.

We set the constants a_k so that the system damps Q to zero on some timescale τ_d that we choose. We assume that τ_d is longer than the interval t_{i + 1} − t_i defined in equation (11), so that we can approximate Q(t) and λ(t) as smooth functions rather than as functions with piecewise constant Nth derivatives. Under this assumption, we write

$$\begin{equation} U(t) = {\rm d}^N\lambda /{\rm d} t^N. \end{equation} \tag{ 17 }$$

We also write

$$\begin{eqnarray} {\rm d}Q/{\rm d}t &= (\partial Q/\partial \lambda ) ({\rm d}\lambda /{\rm d}t), \nonumber \\ &= -{\rm d}\lambda /{\rm d}t. \end{eqnarray} \tag{ 18 }$$

In the first line we have neglected the time dependence of other parameters besides λ that enter into Q under the assumption that the control system timescale is shorter than the timescales of the quantities that we want to control. In the second line we have used equation (13) and we have assumed that Q is small. Similarly we write

$$\begin{eqnarray} {\rm d}^2Q/{\rm d}t^2 &= \frac{\partial Q}{\partial \lambda } \frac{{\rm d}^2\lambda }{{\rm d}t^2} +\frac{\partial ^2 Q}{\partial \lambda ^2} \bigg(\frac{{\rm d}\lambda }{{\rm d}t}\bigg)^2, \nonumber \\ &= \frac{\partial Q}{\partial \lambda } \frac{{\rm d}^2\lambda }{{\rm d}t^2}, \nonumber \\ &= -\frac{{\rm d}^2\lambda }{{\rm d}t^2}, \end{eqnarray} \tag{ 19 }$$

where in the second line we have retained only terms of linear order in dQ/dt (and therefore in dλ/dt). For higher derivatives we continue to retain only terms linear in Q and its derivatives, so from equation (17) we obtain

$$\begin{eqnarray} U(t) &= {\rm d}^N\lambda /{\rm d} t^N,\nonumber \\ &= -{\rm d}^NQ/{\rm d} t^N. \end{eqnarray} \tag{ 20 }$$

For the PID controller and N = 2, combining equations (15) and (20) yields

$$\begin{equation} -\frac{{\rm d}^2 Q}{{\rm d} t^2} = a_0\int Q(t)\,{\rm d}t + a_1\,Q(t) + a_2\frac{{\rm d}{Q}}{{\rm d}t}. \end{equation} \tag{ 21 }$$

If we choose $a_0=1/\tau _d^3$, $a_1=3/\tau _d^2$, and a₂ = 3/τ_d, then the solution to equation (21) will be exponentially damped on the timescale τ_d,

$$\begin{equation} Q\propto {\rm e}^{-t/\tau _d}. \end{equation} \tag{ 22 }$$

The same exponential damping holds for the PD controller, equation (16), for appropriate choices of the parameters a_k. For K = 2 and N = 3, the parameters a_k are identical to those in the PID case above.

The PID controller, equation (15), is computed by measuring the control error Q, its time integral, and its time derivative. The PD controller, equation (16), may require multiple derivatives of Q(t) depending on the order K. Generally only Q, and not its derivatives or integrals, is available in the code at any given time step. The simplest way to compute the derivatives of Q is by finite differencing in time, and the simplest way to compute the integral is by a numerical quadrature.

We measure the control error Q at time intervals of length τ_m, where we choose τ_m < t_{i + 1} − t_i. The measuring time interval τ_m is then used as the time step for finite difference stencils and quadratures.

The measured Q is typically a function of apparent horizon locations or shapes, and this measured Q includes noise caused by the finite resolution of the evolution, and the finite residual and finite number of iterations of the apparent horizon finder. Taking a numerical derivative of Q amplifies the noise, and then this noise is transferred to the control signal via equations (15) or (16), and then to the map. If the noise amplitude is too large, the control system will become unstable. The PID controller generally handles noise better than the PD controller for two reasons. First, each successive numerical derivative amplifies noise even further, so using only one derivative instead of two (or more) results in a more accurate control signal. Second, the inclusion of an integral term acts to further smooth the control signal.

In some cases, however, the noise in Q can be problematic even for the PID controller. In these cases we implement direct averaging of the control error in one of two ways: (1) we perform a polynomial fit of order N to the previous M measurements of the control error, where M > N, or (2) we perform an exponentially-weighted average, with timescale τ_avg, of all previous control error measurements and their derivatives and integrals. The latter is our preferred method, which we describe in detail in appendix A.

In the previous section we have identified four timescales relevant for each control system. The first is the damping timescale τ_d; this describes how quickly the control error Q(t) falls to zero, and therefore how quickly the map parameter λ(t) approaches its desired value. The second timescale is the control interval Δt_i := t_{i + 1} − t_i, which represents how often the Nth derivative of the function λ(t) in equation (11) is updated. The third is the measurement timescale τ_m, which indicates how often the control error Q is measured. The fourth is the averaging timescale τ_avg, which is used to smooth the control error Q (and its derivatives and integrals) for use in computing the control signal U(t). These timescales are not all independent; for example we have assumed Δt_i < τ_d in deriving equation (21), and we have assumed τ_m < Δt_i so that we can obtain smooth measurements of the derivatives of Q.

Because binary black hole evolutions are nonlinear dynamic systems, we adjust the damping timescale, τ_d, throughout the simulation. We then set the three other timescales, Δt_i, τ_m, and τ_avg, based on the current value of the damping timescale τ_d, as we now describe.

The timescale τ_d should be shorter than the timescale on which the physical system changes; otherwise, the control system cannot adjust the map parameters quickly enough to respond to changes in the system. But if the timescale τ_d is too small, then the measurement timescale τ_m on which we compute the apparent horizon must also be small, meaning that frequent apparent horizon computations are needed; this is undesirable because computing apparent horizons is computationally expensive. We would like to adjust τ_d in an automatic way so that it is relatively large during the binary black hole inspiral, decreases during the plunge and merger, and increases again as the remnant black hole rings down. In the canonical language of control theory, we would like to implement 'gain scheduling' [19], tuning the behavior of the control system for different operating regimes.

We do this as follows: for all map parameters (except the ℓ = 0 component of the horizon shape map, which is treated differently; see section 5.3), the damping timescale is a generic function of Q and $\dot{Q}$, i.e., the error in its associated map parameter and its derivative. Whenever we adjust the control signal U(t) at interval t_i, we also tune the timescale in the following way:

$$\begin{equation} \tau _d^{i+1} = \beta \tau _d^i, \end{equation} \tag{ 23 }$$

where typically

$$\begin{equation} \beta = \left\lbrace \begin{array}{@{}ll@{}}0.99, & {\rm if }\ \dot{Q}/Q > -1/2\tau _d \ {\rm and }\ |Q|\ {\rm or }\ |\dot{Q}\tau _d| > Q_t^{\rm Max} \\ 1.01, & {\rm if }\ |Q| < Q_t^{\rm Min}\ {\rm and }\ |\dot{Q}\tau _d| < Q_t^{\rm Min} \\ 1, & {\rm otherwise.} \end{array}\right. \end{equation} \tag{ 24 }$$

Here $Q_t^{\rm Min}$ and $Q_t^{\rm Max}$ are thresholds for the control error Q. The idea is to keep $|Q|<Q_t^{\rm Max}$ so that the map parameters are close to their desired values, but to also keep $|Q|>Q_t^{\rm Min}$ because an unnecessarily small |Q| means unnecessarily small timescales and therefore a large computational expense (because the apparent horizon must be found frequently). For binary black hole simulations where the holes have masses M_A and M_B, we find that the following choices work well:

$$\begin{equation} Q_t^{\rm Max} = \frac{2\times 10^{-3}}{M_A/M_B + M_B/M_A} \end{equation} \tag{ 25 }$$

$$\begin{equation} Q_t^{\rm Min} = \frac{1}{4}Q_t^{\rm Max}. \end{equation} \tag{ 26 }$$

Once we have adjusted the timescale τ_d for every control system, we then use these timescales to choose the times t_{i + 1} in equation (11) at which we update the polynomial coefficients $\lambda _i^n$ in the map parameter expression λ(t),

$$\begin{equation} \Delta t_i := t_{i+1} - t_i = \alpha _d\min (\tau _d), \end{equation} \tag{ 27 }$$

where typically α_d = 0.3, and the minimum is taken over all map parameters (except for the ℓ = 0 component of the horizon shape map). This ensures that the coefficients $\lambda _i^n$, are updated faster than the physical system is changing, and faster than the control system is damping. For α_d too large, we find that the control system becomes unstable.

We also use the timescale τ_d to choose the interval τ_m at which we measure the control error. For many map parameters, the associated control error is a function of apparent horizon quantities, which is why we desire τ_m to be as large as possible. But a certain number of measurements are needed for each Δt_i so that the control signals, defined in equations (15) and (16), are sufficiently accurate and the control system is stable. We choose

$$\begin{equation} \tau _m = \alpha _m \Delta t_i, \end{equation} \tag{ 28 }$$

where α_m is typically between 0.25 and 0.3. In other words, we measure the control error three or four times before we update the control signal. This also ensures that the averaging timescale is greater than the measurement timescale, as we typically choose τ_avg ∼ 0.25τ_d.

The scaling map, $\mathcal {M}_{\rm Scaling}$, causes the grid to shrink or expand (and the excision boundaries to move respectively closer together or farther apart) in the inertial frame, thereby allowing the grid to follow the two black holes as their separation changes.

This map transforms the radial coordinate with respect to the origin (i.e., the center of the outermost sphere in figure 4), such that the region near the black holes is scaled uniformly by a factor a and the outer boundary is scaled by a factor b,

$$\begin{equation} R \mapsto aR + (b-a) R^3/R^2_{\rm OB}. \end{equation} \tag{ 41 }$$

Here R is the radial coordinate, R_OB is the radius at the outer boundary, a is a parameter that will be determined by a control system, and b is another parameter that will be determined empirically.

In [1], the scaling map was simply xⁱ↦a(t)xⁱ, which is recovered by equation (41) if b = a. In that case, as the black holes inspiral together and a decreases, the outer boundary of the grid decreases as well. For long evolutions, the outer boundary decreases so much that we can no longer extract gravitational radiation far from the hole. The addition of b to the map alleviates this difficulty, allowing the motion of the outer boundary to be decoupled from the motion of the holes. We choose b by an explicit functional form

$$\begin{equation} b(t) = 1-10^{-6} t^3/(2500+t^2), \end{equation} \tag{ 42 }$$

which is designed to keep the outer boundary from shrinking rapidly, but to allow the boundary to move inward at a small speed, so that zero-speed modes are advected off the grid (and thus need no boundary condition imposed on them). We choose a cubic function of time in equation (42) because we have found that at least the first two time derivatives of b(t) must vanish initially, or else a significant ingoing pulse of constraint violations is produced at the outer boundary.

We control the scale factor a of this map using the error function

$$\begin{equation} Q_a = a({\rm d}x^0 - 1), \end{equation} \tag{ 43 }$$

where

$$\begin{equation} {\rm d}x^i = \frac{\hat{\xi }^i_A-\hat{\xi }^i_B}{C^0_A - C^0_B}. \end{equation} \tag{ 44 }$$

We assume that the separation vector $C^i_A-C^i_B$ in the grid frame is parallel to the x-axis. The idea is that the distorted-frame separation of the horizon centers along the x-axis is driven to be the same as the separation of the excision boundary centers.

To show that the control system for a obeys equation (13), we consider the change in Q_a under variations of a with all other maps held fixed, and with the inertial-coordinate centers $\bar{\xi }^i_H$ of the horizons held fixed⁴. This means that the distorted-frame centers of the horizons $\hat{\xi }^i_H$ appearing in equation (44) change with variations of a. The second term in equation (41) is small when evaluated near the horizon where (R/R_OB)² ≪ 1, so we can write the action of the scaling map as

$$\begin{equation} \bar{\xi }^i_H = a\hat{\xi }^i_H, \end{equation} \tag{ 45 }$$

and therefore under variations δa,

$$\begin{equation} 0=\delta \bar{\xi }^i_H = \hat{\xi }^i_H \,\delta a + a\,\delta \hat{\xi }^i_H. \end{equation} \tag{ 46 }$$

Taking variations of equation (43), using equation (46) to substitute for $\delta \hat{\xi }^i_H$, and noting that the excision-boundary centers $C^i_H$ are constants and do not vary with a, we obtain

$$\begin{equation} \delta Q_a = -\delta a, \end{equation} \tag{ 47 }$$

so that equation (13) is satisfied.

To verify that the more general decoupling equation, equation (14), is satisfied for the scaling map, we must show, in addition to equation (47), that the quantity dx⁰ defined in equation (44) is invariant under the action of the other maps, at least in the limit that all the control errors Q are small. We consider each map in turn. The translation map moves both apparent horizons together, so it leaves dx⁰ invariant. Changes in the rotation map parameters will change dx⁰ only by an amount proportional to the control errors of the rotation map (see equations (53) and (54) below). The skew map (below) leaves the centers of the excision boundaries invariant. Because the intent of the rotation, translation, and scaling maps is to drive the centers of the apparent horizons toward the centers of the excision boundaries, this means that the skew map changes dx⁰ by an amount proportional to the control errors of the rotation, translation, and scaling maps. Finally, the shape and CutX maps connect the distorted and grid frames, so they cannot affect dx⁰.

The rotation map, $\mathcal {M}_{\rm Rotation}$, is a rigid 3D rotation about the origin that tracks the orbital phase and precession of the system,

$$\begin{equation} x^i \mapsto \mathcal {R}^i_j x^j, \;\;\; \mbox{where} \;\; \mathcal {R} = \left( \begin{array}{@{}ccc@{}}\cos \vartheta \cos \varphi & -\sin \vartheta & \cos \vartheta \sin \varphi \\ \sin \vartheta \cos \varphi & \cos \vartheta & \sin \vartheta \sin \varphi \\ -\sin \varphi & 0 & \cos \varphi \end{array}\right). \end{equation} \tag{ 48 }$$

The pitch and yaw map parameters (φ, ϑ) are controlled so as to align the line segment connecting the apparent horizon centers with the distorted-frame x-axis. Note that the map parameters (φ, ϑ) are functions of time, and are not to be confused with the polar coordinates (ϕ_H, θ_H) centered about each excision boundary. Then the control error is given by

$$\begin{equation} Q_{\varphi } = -\frac{\hat{\xi }^2_A-\hat{\xi }^2_B}{\hat{\xi }^0_A-\hat{\xi }^0_B} \end{equation} \tag{ 49 }$$

$$\begin{equation} Q_{\vartheta } = \frac{\hat{\xi }^1_A-\hat{\xi }^1_B}{\big(\hat{\xi }^0_A-\hat{\xi }^0_B\big) \cos \varphi }. \end{equation} \tag{ 50 }$$

In the case of extreme precession where φ → π/2, these equations are insufficient because Q_ϑ diverges. Our solution is to use quaternions, which avoid this singularity (for a complete discussion, see [21]).

The control errors Q_ϑ and Q_φ obey equation (14). To show this, consider variations of ϑ and φ with other maps held fixed, and with the inertial-coordinate centers $\bar{\xi }^i_H$ of the horizons held fixed. The rotation map, equation (48), implies that under these variations,

$$\begin{equation} \delta \bar{\xi }^i_H = 0 = \mathcal{R}^i_j \delta \hat{\xi }^j_H + (\delta \mathcal{R})^i_j \hat{\xi }^j_H. \end{equation} \tag{ 51 }$$

Multiplying this equation by $\mathcal{R}^{-1}$ we obtain

$$\begin{equation} \delta \hat{\xi }^i_H = -(\mathcal{R}^{-1} \delta \mathcal{R})^i_j \hat{\xi }^j_H, \end{equation} \tag{ 52 }$$

which yields

$$\begin{eqnarray} \fl \frac{\partial \hat{\xi }^i_H}{\partial \varphi } &= \mathcal {A}^i_j\hat{\xi }^j_H, \;\;\; \mbox{where} \;\; \mathcal {A} := -\mathcal {R}^{-1}\frac{\partial \mathcal {R}}{\partial \varphi } = \left(\begin{array}{@{}rrr@{}}0&0&-1\\ 0&0&0\\ 1& 0&0 \end{array}\right) \end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray} \fl \frac{\partial \hat{\xi }^i_H}{\partial \vartheta } &= \mathcal {B}^i_j\hat{\xi }^j_H, \;\;\; \mbox{where} \;\; \mathcal {B} := -\mathcal {R}^{-1}\frac{\partial \mathcal {R}}{\partial \vartheta }= \left(\begin{array}{@{}ccc@{}}0&\cos \varphi &0\\ -\cos \varphi &0&-\sin \varphi \\ 0&\sin \varphi &0 \end{array}\right). \end{eqnarray} \tag{ 54 }$$

We can now verify equation (14) for the special case where the indices a and b in equation (14) are either ϑ and φ. This result is obtained in a straightforward way by differentiating equations (49) or (50) with respect to ϑ or φ, and substituting equations (53) or (54).

In addition, the control errors Q_φ and Q_ϑ are independent (to leading order in the control errors) of changes in the parameters of the other maps that make up $\mathcal M$. Variations of equations (49) and (50) with respect to the scaling map parameter a are zero, because both the numerator and denominator of Q_φ and Q_ϑ scale in the same way with a. Similarly, Q_φ and Q_ϑ are independent of the translation map, since both apparent horizons are translated by the same amount. The skew map can change Q_φ and Q_ϑ, but only by an amount proportional to control errors, because the skew map leaves $C^i_H$ invariant and other maps ensure that $\hat{\xi }^i_H$ are close (within a control error) to $C^i_H$. The shape and CutX maps cannot affect Q_φ and Q_ϑ because those maps connect the grid and distorted frames (and therefore they cannot change the distorted-frame horizon centers $\hat{\xi }^i_H$).

The translation map, $\mathcal {M}_{\rm Translation}$, transforms the Cartesian coordinates, xⁱ, according to

$$\begin{equation} x^i \mapsto x^i + f(R) T^i, \end{equation} \tag{ 55 }$$

where f(R) is a Gaussian centered on the origin with a width set such that f(R) falls off to machine precision at the outer boundary radius, and Tⁱ(t) are translation parameters that are adjusted by a control system.

The translation map moves the grid to account for any drift of the 'center of mass' of the system (as computed assuming point masses at the apparent horizon centers) in the inertial frame. This drift can be caused by momentum exchange between the black holes and the surrounding gravitational field [22, 23], by anisotropic radiation of linear momentum to infinity, or by linear momentum in the initial data. This is the third map (the other two are rotation and scaling) that drives the centers of the apparent horizons toward the centers of the excision boundaries. The apparent horizon centers $\hat{\xi }^i_A$ and $\hat{\xi }^i_B$ represent six degrees of freedom: one is fixed by the scaling map, two by rotation, and three by translation.

The control errors will be more complicated than for the other control systems because translation and rotation do not commute. We define the control errors $Q_T^i$ for each of the translation directions i = 0, 1, 2 as

$$\begin{equation} Q_T^i = a \mathcal{R}^i_j \big[\hat{\xi }^j_B + \mathcal{P}^j_k \big(\hat{\xi }^k_A-\hat{\xi }^k_B\big)\big], \end{equation} \tag{ 56 }$$

where $\mathcal{R}$ is the rotation matrix in equation (48), and $\mathcal{P}$ is some matrix yet to be determined, which may depend on the rotation parameters (φ, ϑ) and on the constants $C^i_H$, but which may not depend on the translation parameters Tⁱ. This control error must have the property that $Q_T^i=0$ when $\hat{\xi }^i_H=C^i_H$, so our first restriction on $\mathcal{P}$ is that it satisfies

$$\begin{equation} 0 = C^i_B + \mathcal{P}^i_j\big(C^j_A-C^j_B\big). \end{equation} \tag{ 57 }$$

To check equation (14), we note that near the black holes we can neglect the last term in equation (41) and we can use f(R) ∼ 1 in equation (55), so that the apparent horizon centers in the inertial and distorted frames are related by

$$\begin{equation} a\mathcal {R}^i_j\hat{\xi }^j_H = \bar{\xi }^i_H - T^i. \end{equation} \tag{ 58 }$$

Inserting equation (58) into equation (56), we obtain

$$\begin{equation} Q_T^i = \bar{\xi }^i_B - T^i - (\mathcal{R}\mathcal{P}\mathcal{R}^{-1})^i_j \big(\bar{\xi }^j_A-\bar{\xi }^j_B\big). \end{equation} \tag{ 59 }$$

The second term in equation (59) is the only term that depends on the translation parameter Tⁱ, so we have

$$\begin{equation} \frac{\partial Q_T^i}{\partial T^j} = -\delta _{ij} \end{equation} \tag{ 60 }$$

When varying map parameters, the inertial-frame horizon centers remain fixed, so the only other term in equation (59) that depends on map parameters is the last term, which depends on (ϑ, φ) because of the rotation matrices and because of the (ϑ, φ) dependence in $\mathcal{P}$. We can therefore write the variation of $Q_T^i$ with respect to (ϑ, φ) as

$$\begin{equation} \frac{\partial Q_T^i}{\partial W} = -\frac{\partial }{\partial W}(\mathcal{R}\mathcal{P}\mathcal{R}^{-1})^i_j \big(\bar{\xi }^j_A-\bar{\xi }^j_B\big), \end{equation} \tag{ 61 }$$

$$\begin{equation} = -\frac{\partial }{\partial W}(\mathcal{R}\mathcal{P}\mathcal{R}^{-1})^i_j a \mathcal{R}^j_k \big(\hat{\xi }^k_A-\hat{\xi }^k_B\big), \end{equation} \tag{ 62 }$$

where W stands for either ϑ or φ, and where in the last line we have used equation (58). To obey equation (14), $\partial Q_T^i/\partial W$ must either be zero, or on the order of a control error Q. For i = 1, 2, the quantity $(\hat{\xi }^i_A-\hat{\xi }^i_B)$ is proportional to the control error of the rotation map, equations (49) and (50), so these terms can be neglected in equation (62). However, for i = 0, $(\hat{\xi }^i_A-\hat{\xi }^i_B)$ is not proportional to a control error; instead $(\hat{\xi }^0_A-\hat{\xi }^0_B)$ is driven to a constant finite value of $C^0_A-C^0_B$ by the scaling control system. Therefore, in order to satisfy equation (14), we require

$$\begin{equation} -\frac{\partial }{\partial W}(\mathcal{R}\mathcal{P}\mathcal{R}^{-1})\mathcal{R} \left( \begin{array}{@{}c@{}}1\\ 0\\ 0 \end{array} \right) = 0. \end{equation} \tag{ 63 }$$

We find that we can satisfy both equations (57) and (63) by choosing

$$\begin{equation} \mathcal{P} = \frac{1}{C^0_B-C^0_A} \left( \begin{array}{@{}ccc@{}}C^0_B & -C^1_B & -C^2_B\\ C^1_B & C^0_B + C^2_B \tan \varphi & 0 \\ C^2_B & -C^1_B \tan \varphi & C^0_B \end{array} \right). \end{equation} \tag{ 64 }$$

Here we have again assumed that the separation between the centers of the excision boundaries is parallel to the x-axis, i.e., $C^1_A=C^1_B$ and $C^2_A=C^2_B$.

In figure 4, a prominent feature of the domain decomposition is a plane (a vertical line in the two-dimensional figure) that is perpendicular to the grid-frame x-axis and lies between the two excision boundaries A and B. We call this plane the 'cutting plane'.

The skew map, $\mathcal {M}_{\rm Skew}$, acts on the distorted-frame coordinates, in which the cutting plane is still perpendicular to the x-axis. The skew map leaves the coordinates y, z unchanged, but changes the x-coordinate in order to give a skewed shape to the cutting plane, as shown in figure 5. Let $x^i_C$ be the intersection point of the line segment connecting the excision boundary centers and the cutting plane. The action of the skew map is defined as

$$\begin{equation} x^0 \mapsto x^0 + W \sum _{j=1,2} \tan [ \Theta ^j(t) ] \cdot \big(x^j-x^j_C\big), \end{equation} \tag{ 65 }$$

$$\begin{equation} x^j \mapsto x^j, \;\; \mbox{for} \;\; j=1,2 \end{equation} \tag{ 66 }$$

where W is a radial Gaussian function centered around $x^i_C$, and the angles Θ^j(t) are the time-dependent map parameters. For j = 1, 2 the parameter Θ^j(t) is the angle between the mapped and unmapped x^j-axis when projected into the $x^{3-j}=x^{3-j}_C$ plane. Note that the line intersecting $x^i_C$ and parallel to the x-axis is left invariant by the skew map⁵. The width of the Gaussian W is set such that W is below machine precision at the innermost wave extraction sphere. This implies that the spherical-shell subdomains used to evolve the metric in the wave zone will not be affected by the skew map.

Zoom In Zoom Out Reset image size

The purpose of the skew map is to (as much as possible) align the cutting plane with the surfaces of the apparent horizons in the region where the surfaces are closest to the cutting plane. We derive the control system in charge of setting the parameters Θ^j by the following condition: we demand that the angle between the mapped cutting plane and the x-axis at $x^i_C$ be driven to the (weighted) average of the angles at which the mapped apparent horizons intersect the same x-axis.

The input coordinates to the skew map are the coordinates in the distorted frame. For each horizon, we therefore measure an angle in the distorted frame as follows: let $x^{\rm Int}_H$ be the distorted-frame x-coordinate at which apparent horizon H intersects the line segment connecting the centers of the excision boundaries. For j = 1, 2, we calculate the normal to the surface at this intersection point. This normal is projected into the $x^{3-j}=x^{3-j}_C$ plane. We define $\Theta _H^j$ as the angle between the projected normal and the x-axis. Thus, projecting the normal into the y = const. plane gives $\Theta _H^z$ and vice versa (see figure 6). We then compute a weighted average of the $\Theta ^j_H$ such that the horizon closer to the cutting plane has a larger effect on the skew angles,

$$\begin{equation} \Theta ^j_{\rm Avg} = \frac{ w_A \Theta _A^j + w_B \Theta _B^j }{w_A W(C_A) + w_B W(C_B)}, \end{equation} \tag{ 67 }$$

where w_A, w_B are averaging weights,

$$\begin{equation} w_H = \exp \left[-\frac{x_C^0 - x^{\rm Int}_H}{x_C^0 - C^0_H}\right]. \end{equation} \tag{ 68 }$$

Here $C^i_H$ is the center of the excision boundary H.

Zoom In Zoom Out Reset image size

We measure $\Theta ^j_{\rm Avg}$ in the distorted frame, and in this frame the cutting plane is always normal to the x-axis. Therefore, thinking about the desired result in the distorted frame, we see that the control system for the skew map should drive $\Theta ^j_{\rm Avg}$ to zero. This leads us to consider the following control error for the skew angles: $Q^j_{\Theta }=\Theta ^j_{\rm Avg}$. Assuming that the function W is unity near the apparent horizons, we find that

$$\begin{equation} \frac{\partial \Theta ^j_{\rm Avg}}{\partial \Theta ^j} = -1+\mathcal{O}(Q), \end{equation} \tag{ 69 }$$

in agreement with equation (13). In deriving equation (69), it is helpful to observe that the partial derivative in equation (69) is taken with the inertial-frame apparent horizon held fixed. For a fixed inertial-frame horizon, the only map that can change the shape (as opposed to merely the center) of the distorted-frame horizon (and thus $\Theta ^j_{\rm Avg}$) is the skew map.

We do not use $Q^j_{\Theta } =\Theta ^j_{\rm Avg}$ for the entire evolution, however, because at early times, when the coordinate distance between the apparent horizons is larger than their combined radii, the skew map is not needed. Furthermore, the skew map can cause difficulties early during the run, especially during the 'junk radiation' phase when the horizons are oscillating in shape. For this reason, we gradually turn on the skew map as the black holes approach each other. This is done by defining a roll-on function g that is zero when the horizons are far apart, and one when they are close together. This roll-on function is defined as

$$\begin{equation} g = \frac{1}{2}\left[ 1-\tanh \left( 10\frac{x^{\rm Int}_A - x^{\rm Int}_B}{C^0_A-C^0_B}-5 \right) \right]. \end{equation} \tag{ 70 }$$

For values g < 10⁻³ the skew map is turned off completely; this is not strictly necessary, but it saves some computation. Given the function g, we define the control error as

$$\begin{equation} Q^j_{\Theta }=g\Theta ^j_{\rm Avg}-(1-g)\Theta ^j. \end{equation} \tag{ 71 }$$

This control error drives the skew angles to zero when the black holes are far apart and drives $\Theta ^j_{\rm Avg}$ to zero as they approach each other.

We define the shape map $\mathcal {M}_{\rm Shape}$ as:

$$\begin{equation} x^i \mapsto x^i\bigg(1-\sum _H\frac{f_H(r_H,\theta _H,\phi _H)}{r_H} \sum _{\ell m} Y_{\ell m}(\theta _H,\phi _H)\lambda ^H_{\ell m}(t)\bigg). \end{equation} \tag{ 72 }$$

The index H goes over each of the two excised regions A and B, and the map is applied to the grid-frame coordinates. The polar coordinates (r_H, θ_H, ϕ_H) centered about excised region H are defined by equations (32)–(34), the quantities Y_ℓm(θ_H, ϕ_H) are spherical harmonics, and $\lambda ^H_{\ell m}(t)$ are expansion coefficients that parameterize the map near excision region H. The function f_H(r_H, θ_H, ϕ_H) is chosen to be unity near excision region H and zero near the other excision region, so that the distortion maps for the two black holes are decoupled. Specifically, f_H(r_H, θ_H, ϕ_H) is determined as illustrated in figure 4. For excision region A in the figure, f_A(r_A, θ_A, ϕ_A) = 1 between the excision boundary and the magenta surface, it falls linearly to zero between the magenta and red surfaces, and it is zero everywhere outside the red surface. This means that the gradient of f_A(r_A, θ_A, ϕ_A) is discontinuous on the red surface, the magenta surface, and the blue surfaces in the figure. Because we ensure that these discontinuities occur on subdomain boundaries, they cause no difficulty with using spectral methods. Around excision region B, f_B(r_B, θ_B, ϕ_B) is chosen similarly.

In previous implementations of the shape map [3], the functions f_H(r_H, θ_H, ϕ_H) were chosen to be smooth Gaussians centered around each excision boundary rather than to be piecewise linear functions. We find smooth Gaussians to be inferior for two reasons. The first is that piecewise linear functions are easier and faster to invert (the inverse map is required for interpolation to trial solutions during apparent horizon finding). The second is that for smooth Gaussians, it is necessary to choose the widths of the Gaussians sufficiently narrow so that the Gaussian for excision region A does not overlap the Gaussian for excision region B and vice versa, so that the maps and control systems for A and B remain decoupled. However, decreasing the width of the Gaussians increases the Jacobians of the map, producing coordinates that are stretched and squeezed nonuniformly. We found that this form of 'grid-stretching' significantly increased the computational resources required to resolve the solution to a given level of accuracy. In other words, with smooth Gaussians we were forced to add computational resources just to resolve the large Jacobians.

A map very similar to equation (72) is also described in [4]. The difference compared to this work is the choice of f_H(r_H, θ_H, ϕ_H), which corresponds to a different choice of domain decomposition. The control system used to choose the map parameters $\lambda ^H_{\ell m}(t)$ in [4] is also different than what is described here.

We control the expansion coefficients $\lambda ^H_{\ell m}(t)$ of the shape map, equation (72), so that each excision boundary, as measured in the intermediate frame $(t,\hat{x}^i)$, has the same shape as the corresponding apparent horizon. In other words, we desire

$$\begin{equation} \frac{\hat{r}_{\rm AH}}{\langle \hat{r}_{\rm AH}\rangle } = \frac{\hat{r}_{\rm EB}}{\langle \hat{r}_{\rm EB}\rangle }, \end{equation} \tag{ 73 }$$

where for brevity we have dropped the index H that labels the excision boundary. Here the angle brackets mean averaging over angles, $\hat{r}_{\rm AH}(\hat{\theta },\hat{\phi })$ is the radial coordinate of the apparent horizon defined in equation (36), and $\hat{r}_{\rm EB}(\theta ,\phi )$ is the radial coordinate of the excision boundary, which from equation (72) can be written

$$\begin{equation} \hat{r}_{\rm EB} = r_{\rm EB} - \sum _{\ell m} Y_{\ell m}(\hat{\theta },\hat{\phi })\lambda _{\ell m}(t). \end{equation} \tag{ 74 }$$

Here r_EB is the radius of the (spherical) excision boundary in the grid frame. In deriving equation (74) we have used the relations $\mathcal{M}_{\rm CutX}=\boldsymbol {1}$, f(r, θ, ϕ) = 1, $\hat{\theta }=\theta$, and $\hat{\phi }=\phi$, which hold on the excision boundary.

Combining equations (36), (73), and (74) yields

$$\begin{equation} \frac{r_{\rm EB} - \sum _{\ell m} Y_{\ell m}(\hat{\theta },\hat{\phi }) \lambda _{\ell m}(t)}{r_{\rm EB} - Y_{00} \lambda _{00}(t)} = \frac{\sum _{\ell m} \hat{S}_{\ell m} Y_{\ell m}(\hat{\theta }, \hat{\phi })}{\hat{S}_{00} Y_{00}}, \end{equation} \tag{ 75 }$$

which we can satisfy by demanding that

$$\begin{equation} \lambda _{\ell m}(t)+ \hat{S}_{\ell m} \frac{r_{\rm EB} - Y_{00} \lambda _{00}(t)}{\hat{S}_{00} Y_{00}}=0, \quad \ell >0. \end{equation} \tag{ 76 }$$

Therefore, given an apparent horizon and given a value of λ₀₀, a control system can be set up for each (ℓ, m) pair with ℓ > 0, and the corresponding control errors are

$$\begin{equation} Q_{\ell m} = -\lambda _{\ell m}(t)- \hat{S}_{\ell m} \frac{r_{\rm EB} - Y_{00} \lambda _{00}(t)}{\hat{S}_{00} Y_{00}}, \quad \ell >0. \end{equation} \tag{ 77 }$$

Driving these Q_ℓm to zero produces an excision boundary that matches the shape of the corresponding apparent horizon.

Equation (77) determines λ_ℓm(t) only for ℓ > 0, and equations (73)–(77) can be satisfied for arbitrary values of the remaining undetermined map coefficient λ₀₀(t). Determination of λ₀₀(t) is complicated enough that it is described in its own section, section 5.

Note that equation (77) does not satisfy equation (14) because $\partial Q_{\ell m}/\partial \lambda _{00}=\hat{S}_{\ell m}/\hat{S}_{00}$, which does not vanish even if all the control errors are zero. For small distortions, this coupling between λ₀₀ and Q_ℓm seems to cause little difficulty. However, at times when the shapes and sizes of the horizons change rapidly (e.g., during the initial 'junk radiation' phase, after the transition to a single excised region when a common apparent horizon forms, and after rapid gauge changes), simulations using equation (77) exhibit relatively large and high-frequency oscillations in Q_ℓm that usually damp away but occasionally destabilize the evolution. Construction of a control system in which all Q_ℓm are fully decoupled will be addressed in a future work.

The map $\mathcal {M}_{\rm CutX}$ applies a translation along the grid-frame x-axis in the vicinity of the black holes, but without moving the excision boundaries (or the surrounding spherical shells) themselves. The action of the map is shown in figure 7.

Zoom In Zoom Out Reset image size

The goal of this map is to allow for a slight motion of the cutting plane toward the smaller excision boundary. This is important for binary black hole systems with mass ratio q ≳ 8. As such a binary gets closer to merger, the inertial-frame coordinate distance between each excision boundary and the cutting plane decreases. Eventually, this distance falls to zero for the larger excision boundary, producing a coordinate singularity in which the Jacobian of one of the other maps (often the shape map) becomes infinite. By pushing the cutting plane toward the smaller excision boundary, the map $\mathcal {M}_{\rm CutX}$ avoids this singularity. Even for evolutions in which the inertial-coordinate distance between the larger excision boundary and the cutting plane remains finite but becomes small, the map $\mathcal {M}_{\rm CutX}$ prevents large Jacobians from developing and thus increases numerical accuracy (because it is no longer necessary to add computational resources to resolve the large Jacobians).

Figure 8 shows the domain decomposition in the inertial frame for a binary black hole simulation with q = 8. The top panel shows the case without $\mathcal {M}_{\rm CutX}$, and the compressed grid near the larger excision boundary is evident. The lower panel shows the case with $\mathcal {M}_{\rm CutX}$, which removes the extreme grid compression. Looking at the Jacobian of the mapping from the inertial frame to the grid frame shortly before merger, we find that the infinity norm of the determinant of the Jacobian is twice as large in the case without $\mathcal {M}_{\rm CutX}$.

Zoom In Zoom Out Reset image size

The map $\mathcal {M}_{\rm CutX}$ is written as

$$\begin{equation} x^0 \mapsto x^0 + \rho (x^i) F(t), \end{equation} \tag{ 78 }$$

$$\begin{equation} x^j \mapsto x^j, \qquad j=1,2, \end{equation} \tag{ 79 }$$

where F(t) is adjusted by a control system.

We choose ρ(xⁱ) to be zero within either of the spherical shell regions, i.e., inside the magenta spheres around A and B in figure 7, and it is also zero outside the outer magenta sphere in the same figure. We set ρ(xⁱ) = 1 on the solid red boundaries in figure 7, and everywhere between the two solid red boundaries that enclose excision boundary B. Every other region on the grid is bounded by a smooth red boundary on one side and a smooth magenta boundary on the other; in these regions ρ varies linearly between zero and one. The solid blue boundaries are locations (in addition to the solid red and solid magenta boundaries) in which the gradient of ρ is discontinuous. Full details for the calculation of ρ can be found in appendix B.

As with the skew map, the map $\mathcal {M}_{\rm CutX}$ is inactive for most of the inspiral. Let $x^{\rm Exc}_H$ be the distorted-frame x-coordinate of the intersection of the excision boundary H with the line segment connecting the centers of the excision boundaries. This is similar to $x^{\rm Int}_H$ defined earlier; the difference is that $x^{\rm Int}_H$ refers to a point on the apparent horizon and $x^{\rm Exc}_H$ refers to a point on the excision boundary. The quantity $x^{\rm Exc}_H$ is time-dependent because it depends on the shape map.

The map $\mathcal {M}_{\rm CutX}$ is turned off completely as long as

$$\begin{equation} \frac{x^0_C-x^{\rm Exc}_H}{x^{\rm Exc}_H-C^0_H}\ge \frac{1}{2}, \end{equation} \tag{ 80 }$$

where $C^i_H$ are the excision boundary centers, and $x^i_C$ are the coordinates of the intersection of the cutting plane and the line segment connecting the centers of the excision boundaries, as introduced in section 4.4.

Let t₀ be the coordinate time at which the map $\mathcal {M}_{\rm CutX}$ is activated⁶. We designate the target position x = x_T of the cutting plane as

$$\begin{equation} x_{T} = \left. \frac{\Delta x_A \cdot x^{\rm Exc}_A + \Delta x_B \cdot x^{\rm Exc}_B}{\Delta x_A+ \Delta x_B } \right|_{t=t_0}, \end{equation} \tag{ 81 }$$

where

$$\begin{equation} \Delta x_H = \big| x^{\rm Exc}_H-C^0_H \big|, \quad H=A,B. \end{equation} \tag{ 82 }$$

Recall that $x^{\rm Exc}_H$ is time-dependent (as it is measured in the distorted frame), so we save the value of x_T as calculated at the time when the map $\mathcal {M}_{\rm CutX}$ is activated, rather than recalculating x_T at every measurement time.

Now that we have a target x_T, we could designate $x_{T}-x^0_C$ as the target value to the function F(t) by setting $Q_F=-F+x_{T}-x^0_C$ and have the control system drive Q_F to zero, thus driving the x-coordinate of the cutting plane to x_T. However, because we turn on the CutX control system suddenly at time t₀, we must be more careful. Turning on any control system suddenly will produce some transient oscillations, unless the control error and its relevant derivatives and integrals are all initially zero. In the case of the CutX map, which is turned on during a very dynamic part of the simulation when excision boundaries need to be controlled very tightly, these oscillations can prematurely terminate the run by, for example, pushing an excision boundary outside its accompanying apparent horizon.

To turn on $\mathcal {M}_{\rm CutX}$ gradually, we replace x_T by a new time-dependent target function T(t) that gradually approaches x_T at late times but produces a control error Q_F with Q_F = ∂Q_F/∂t = 0 at the activation time t = t₀. We start by estimating the time $t^{\rm XC}_H$ at which $x^{\rm Exc}_H$ will reach the cutting plane, where H is either A or B. This is done by the method described in appendix C. We then let t^XC be the minimum of $t^{\rm XC}_A$ and $t^{\rm XC}_B$. The smooth target function T(t) is then defined by

$$\begin{eqnarray} &&\fl T(t;t_0,x_T,t^{\rm XC}) = x_T + \exp \bigg[-\bigg(\frac{t-t_0}{0.3 \cdot t^{\rm XC}}\bigg)^2\bigg] \nonumber \\ &&\times \bigg [ x^0_C -x_{T} +F(t_0) + (t-t_0) \cdot \frac{\partial F(t) }{ \partial t}\bigg|_{t=t_0} \bigg ] . \end{eqnarray} \tag{ 83 }$$

Designating $T-x^0_C$ as the target for F(t)

$$\begin{equation} Q_F = T-F-x^0_C \end{equation} \tag{ 84 }$$

leads to

$$\begin{equation} Q_F\bigg|_{t=t_0} = 0, \quad \frac{\partial Q_F}{\partial t}\bigg|_{t=t_0} = 0, \end{equation} \tag{ 85 }$$

so at the activation time the control system does not produce transients. Furthermore, in the limit of small Q_F,

$$\begin{equation} x_C^0 + F(t)|_{t=t^{\rm XC}} \approx x_T, \end{equation} \tag{ 86 }$$

i.e., by the time the excision boundary would have touched the cutting plane (and formed a grid singularity), the smooth target function T(t) has approached x_T, and the cutting plane will have reached its designated target location, x_T. As the run proceeds, the behavior of the cutting plane is determined by the map $\mathcal {M}_{\rm CutX}$ while the motion of the excision boundaries is determined by gauge dynamics and the behavior of the other control systems. We continue to monitor the distance between the cutting plane and the excision boundaries, and if it is predicted to touch within a time less than τ_d/0.15, the $\mathcal {M}_{\rm CutX}$ control system is reset, constructing a new target function T(t; t₀, x_T, t^XC), where t₀ is the time of the reset, and x_T, t^XC are also recalculated based on the state of the grid at this reset time.

Each time a new target function T is constructed, the damping time τ_d of the $\mathcal {M}_{\rm CutX}$ control system is set to be t^XC/2.

The control systems responsible for $\mathcal {M}_{\rm CutX}$ and $\mathcal {M}_{\rm Shape}$ are decoupled, as $\mathcal {M}_{\rm CutX}$ controls the location of the cutting plane, leaving the excision boundary unchanged, while $\mathcal {M}_{\rm Shape}$ controls the shape of the excision boundary, leaving the location of the cutting plane unchanged. Recall that $\mathcal {M}_{\rm CutX}$ and $\mathcal {M}_{\rm Shape}$ define the mapping from the grid frame to the distorted frame. As the apparent horizons are found in the distorted frame, and the other maps only depend upon measurements of the horizons, $\mathcal {M}_{\rm CutX}$ and $\mathcal {M}_{\rm Shape}$ are decoupled from the other maps.

Controlling the size of the excision boundary is more complicated than simply keeping the excision boundary inside the apparent horizon. This is because black hole excision requires conditions on the characteristic speeds of the system, and if these conditions are not enforced they are likely to be violated.

The minimum characteristic speed at each excision boundary is given by

$$\begin{equation} v = -\alpha - \bar{\beta }^i \bar{n}_i - \bar{n}_i \frac{\partial \bar{x}^i}{\partial t}, \end{equation} \tag{ 87 }$$

where α is the lapse, $\bar{\beta }^i$ is the shift, and $\bar{n}_i$ is the normal to the excision boundary pointing out of the computational domain, i.e., toward the center of the hole. Here $(\bar{t},\bar{x}^i)$ are the inertial frame coordinates. The first two terms in equation (87) describe the coordinate speed of the ingoing (i.e., directed opposite to $\bar{n}_i$) light cone in the inertial frame, and the last term accounts for the motion of the excision boundary (which is fixed in the grid frame) with respect to the inertial frame.

In our simulations, we impose no boundary condition whatsoever at each excision boundary. Therefore, well-posedness requires that all of the characteristic speeds, and in particular the minimum speed v, must be non-negative; in other words, characteristics must flow into the hole. In practice, if v becomes negative, the simulation is terminated, because a boundary condition is needed, but we do not have one to impose. This can occur even when the excision boundary is inside the horizon. In this case, one might argue that if the simulation is able to continue without crashing (e.g. by becoming unstable inside the horizon), that the solution outside the horizon would not be contaminated. We have not explored this possibility. Instead, we choose to avoid this situation by terminating the code if negative speeds are detected.

Therefore, we would like to control λ₀₀ in such a way that v remains positive. We start by writing v in a way that separates terms that explicitly depend on $\dot{\lambda }_{00}$ from terms that do not. To do this we expand the derivative in the last term of equation (87) as

$$\begin{eqnarray} \frac{\partial \bar{x}^i}{\partial t} &= \frac{\partial \bar{x}^i}{\partial \hat{x}^a} \frac{\partial \hat{x}^a}{\partial t} \nonumber \\ &=\frac{\partial \bar{x}^i}{\partial \hat{t}} \frac{\partial \hat{t}}{\partial t}+ \frac{\partial \bar{x}^i}{\partial \hat{x}^j} \frac{\partial \hat{x}^j}{\partial t} \nonumber \\ &=\frac{\partial \bar{x}^i}{\partial \hat{t}}+ \frac{\partial \bar{x}^i}{\partial \hat{x}^j} \frac{\partial \hat{x}^j}{\partial t}, \end{eqnarray} \tag{ 88 }$$

where a in the first line of equation (88) is a four-dimensional spacetime index, and the last line of equation (88) follows from $\partial \hat{t}/\partial t = 1$. Inserting this into equation (87) yields

$$\begin{equation} v = -\alpha - \bar{\beta }^i \bar{n}_i - \bar{n}_i \frac{\partial \bar{x}^i}{\partial \hat{t}} - \hat{n}_i \frac{\partial \hat{x}^i}{\partial t}, \end{equation} \tag{ 89 }$$

where $\hat{x}^i$ is the frame that is obtained from the grid frame by applying the distortion map; see equation (30).

We then use equation (72) to rewrite this as

$$\begin{equation} v = -\alpha - \bar{\beta }^i \bar{n}_i - \bar{n}_i \frac{\partial \bar{x}^i}{\partial \hat{t}} + \hat{n}_i \frac{x^i}{r} \sum _{\ell m} Y_{\ell m}(\theta ,\phi )\dot{\lambda }_{\ell m}(t), \end{equation} \tag{ 90 }$$

where we have used the relations f(r, θ, ϕ) = 1 and $\mathcal{M}_{\rm CutX}=\boldsymbol {1}$ when evaluating the distortion map $\mathcal {M}_{\rm Distortion}$ on the excision boundary.

By combining all the terms that do not explicitly depend on $\dot{\lambda }_{00}$ into a quantity v₀, we obtain

$$\begin{equation} v = v_0 + \hat{n}_i \frac{x^i}{r} Y_{00} \dot{\lambda }_{00}. \end{equation} \tag{ 91 }$$

Thus, the characteristic speed v can be thought of as consisting of two parts: one part, v₀, that depends on the position and shape of the excision boundary and the values of the metric quantities there, and another part that depends on the average speed of the excision boundary in the direction of the boundary normal.

We now construct a control system that drives the characteristic speed v to some target speed v_T. This is a control system that controls the derivative quantity $\dot{\lambda }_{00}$, as opposed to directly controlling the map quantity λ₀₀. We choose

$$\begin{equation} Q = ({\rm min}(v)-v_T)/\langle -\Xi \rangle , \end{equation} \tag{ 92 }$$

where

$$\begin{equation} \Xi = \hat{n}_i \frac{x^i}{r} Y_{00}, \end{equation} \tag{ 93 }$$

the minimum is over the excision boundary, and the angle brackets in equation (92) refer to an average over the excision boundary. Note that Ξ < 0 because $\hat{n}_i$ points radially inward and xⁱ/r points radially outward; this means that $\dot{Q} = -\ddot{\lambda }_{00}$, in accordance with our normalization choice for a control system on $\dot{\lambda }_{00}$.

As in our other controlled map parameters, we demand that $\dot{\lambda }_{00}$ is a function with a piecewise-constant second derivative. It is then easy to construct λ₀₀ as a function with a piecewise-constant third derivative.

The control system given by equation (92) with a hand-chosen value of v_T has been used successfully [24, 25] in simulations of high-spin binaries. Figure 9 illustrates why characteristic speed control is crucial for the success of these simulations.

Zoom In Zoom Out Reset image size

The characteristic speed control described above has the disadvantage that it requires a user-specified target value v_T. If v_T is chosen to be too small, then small fluctuations (due to shape control, the horizon finder, or simply numerical truncation error) can cause the characteristic speed to become negative and spoil the simulation.

If v_T is chosen to be too large, the simulation can also fail. To understand why, recall that characteristic speed control achieves the target characteristic speed by moving the excision boundary and thus changing the velocity term in equation (87). So if v > v_T the control system moves the excision boundary radially inward, and if v < v_T the control system moves the excision boundary radially outward. If v_T is too large, the control system can push the excision boundary outward until it crosses the apparent horizon. This halts the evolution because the apparent horizon can no longer be found.

One way to prevent the excision boundary from crossing the horizon is to drive the excision boundary to some constant fraction of the horizon radius, or in other words, drive the quantity d/dt(Δr) to zero, where

$$\begin{equation} \Delta r = 1 - \frac{\langle \hat{r}_{\rm EB}\rangle }{\langle \hat{r}_{\rm AH}\rangle } \end{equation} \tag{ 94 }$$

is the relative difference between the average radius of the apparent horizon (in the intermediate frame) and the average radius of the excision boundary. Using equations (74) and (36), we can write

$$\begin{equation} \frac{{\rm d}}{{\rm d}t}\Delta r = \frac{\dot{\lambda }_{00}}{\hat{S}_{00}} + \frac{\dot{\hat{S}}_{00}}{\hat{S}_{00}} (1-\Delta r), \end{equation} \tag{ 95 }$$

and therefore a control system that adjusts $\dot{\lambda }_{00}$ to achieve d/dt(Δr) = 0 can be obtained by defining

$$\begin{equation} Q = \dot{\hat{S}}_{00}(\Delta r-1)-\dot{\lambda }_{00}. \end{equation} \tag{ 96 }$$

A slight generalization of this control system can be obtained by demanding that ${\rm d}/{\rm d}t(\Delta r) = \dot{r}_{\rm drift}$, where $\dot{r}_{\rm drift}$ is some chosen constant,

$$\begin{equation} Q = \dot{\hat{S}}_{00}(\Delta r-1)-\dot{\lambda }_{00}+\hat{S}_{00} \dot{r}_{\rm drift}. \end{equation} \tag{ 97 }$$

We have not found the control systems defined by equations (96) and (97) to be especially useful on their own. One drawback of these systems is that they do not prevent the minimum characteristic speed v from becoming negative. Instead, we use the horizon tracking control systems described here as part of a more sophisticated control system discussed in the next section.

Here we introduce a means of controlling λ₀₀ that combines the best features of characteristic speed control and horizon tracking. The idea is to continuously monitor the state of the system and switch between different control systems as the evolution proceeds.

At fixed intervals during the simulation that we call 'measurement times', we monitor the minimum characteristic speed v and the relative distance between the horizon and the excision boundary Δr. The goal of the control system is to ensure that both of these values remain positive.

Consider first the characteristic speed v. Using the method described in appendix C, we determine whether v is in danger of becoming negative in the near future, and if so, we estimate the timescale τ_v on which this will occur. Similarly, we also determine whether Δr will soon become negative, and if so we estimate a corresponding timescale τ_Δr. Having estimated both τ_v and τ_Δr, we use these quantities to determine how to control $\dot{\lambda }_{00}$. In particular, we switch between horizon tracking, equation (96), and characteristic speed control, equation (92), based on τ_v and τ_Δr.

We do this by the following algorithm, which favors horizon tracking over characteristic speed control unless the latter is essential. Assume that the control system for $\dot{\lambda }_{00}$ is currently tracking the horizon, equation (96), and that the current damping timescale is some value τ_d. If v is in danger of crossing zero according to the above estimate, and if τ_v < τ_d and τ_v < τ_Δr, we then switch to characteristic speed control, equation (92), we set v_T = 1.01v where v is the current value of the characteristic speed (the factor 1.01 prevents the algorithm from switching back from characteristic speed control to horizon tracking on the very next time step), and we reset the damping time τ_d equal to τ_v. Otherwise we continue to use equation (96), resetting τ_d = τ_Δr if τ_Δr < τ_d.

Now assume that the control system for $\dot{\lambda }_{00}$ is controlling the characteristic speed, equation (92). If v is not in danger of crossing zero, so that we no longer need active control of the characteristic speed, then we switch to horizon tracking, equation (96), without changing the damping time τ_d. If v is in danger of crossing zero, and if Δr is either in no danger of crossing zero or if it will cross zero sufficiently far in the future such that τ_Δr ⩾ σ₁τ_d, then we continue to use equation (92) with τ_d reset to min(τ_d, τ_v) so as to maintain control of the characteristic speed. Here σ₁ is a constant that we typically choose to be about 20.

The more complicated case occurs when both v and Δr are in danger of crossing zero and τ_Δr < σ₁τ_d. In this case, we have two possible methods by which we attempt to prevent Δr from becoming negative. We choose between these two methods based on the values of τ_Δr and τ_d, and also by the value of a quantity that we call the comoving characteristic speed:

$$\begin{equation} \fl v_c = -\alpha - \bar{\beta }^i \bar{n}_i - \bar{n}_i \frac{\partial \bar{x}^i}{\partial \hat{t}} + \hat{n}_i \frac{x^i}{r} \bigg[Y_{00}\dot{\hat{S}}_{00}(\Delta r-1)+ \sum _{\ell >0} Y_{\ell m}(\theta ,\phi )\dot{\lambda }_{\ell m}(t)\bigg]. \end{equation} \tag{ 98 }$$

This quantity is the value that the characteristic speed v would have if horizon tracking (equation (96)) were in effect and working perfectly. Equation (98) is derived by assuming Q = 0 in equation (96), solving for $\dot{\lambda }_{00}$, and substituting this value into equation (90). The reason to consider v_c is that for min(v_c) < 0 (which can happen temporarily during a simulation), horizon tracking is to be avoided, because horizon tracking will drive min(v) toward a negative value, namely min(v_c).

The first method of preventing Δr from becoming negative is to switch to horizon tracking, equation (96); we do this if min(v_c) > 0 and if τ_Δr < σ₂τ_d, where σ₂ is a constant that we typically choose to be 5. The second method, which we employ if σ₂τ_d < τ_Δr < σ₁τ_d or if min(v_c) < 0, is to continue to use characteristic speed control, equation (92), but to reduce the target characteristic speed v_T by multiplying its current value by some fraction η (typically 0.25). By reducing the target speed v_T, we reduce the outward speed of the excision boundary, and this increases τ_Δr. The reason we use the latter method for min(v_c) < 0 is that the former method, horizon tracking, will drive v toward v_c and we wish to keep v positive. The reason we use the latter method for σ₂τ_d < τ_Δr < σ₁τ_d is that in this range of timescales the control system is already working hard to keep v positive, therefore a switch to horizon tracking is usually followed by a switch back to characteristic speed control in only a few time steps, and rapid switches between different types of control make it more difficult for the control system to remain locked. The constants η, σ₁, and σ₂ are adjustable, and although the details of the algorithm (i.e., which particular quantity is controlled at which particular time) are sensitive to the choices of these constants, the overall behavior of the algorithm is not, provided 1 < σ₂ < σ₁.

The overall behavior of this algorithm is to cause $\dot{\lambda }_{00}$ to settle to a state in which both v > 0 and Δr > 0 for a long stretch of the simulation. Occasionally, when either v or Δr threatens to cross zero as a result of evolution of the metric or gauge, the algorithm switches between characteristic speed control and horizon tracking several times and then settles into a new near-equilibrium state in which both v > 0 and Δr > 0.

The actual value of Δr and of v in the near-equilibrium state is unknown a priori, and in particular this algorithm can in principle settle to values of Δr that are either small enough that control system timescales must be very short in order to prevent Δr from becoming negative, or large enough that excessive computational resources are needed to resolve the metric quantities deep inside the horizon. To prevent either of these cases from occurring, we introduce a constant correction velocity $\dot{r}_{\rm corr}$, a nominal target value Δr_target, and a nominal target characteristic speed v_target. We currently choose these to be $\dot{r}_{\rm corr}=0.005$, Δr_target = 0.08, and v_target = 0.08. These quantities are used only to decide whether to replace equation (96) by equation (97) when the algorithm chooses horizon tracking, as described below; these quantities are not used when the algorithm chooses characteristic speed control.

First we consider the case where Δr is too large during horizon tracking. If Δr > 1.1Δr_target, then we replace equation (96) with equation (97), and we use $\dot{r}_{\rm drift}=-\dot{r}_{\rm corr}$. We continue to use $\dot{r}_{\rm drift}=-\dot{r}_{\rm corr}$ until Δr < Δr_target, at which point we set $\dot{r}_{\rm drift}=0$.

The case where Δr is too small during horizon tracking is more complicated because we want to be careful so as to not make the characteristic speeds decrease, since it is more difficult to recover from decreasing characteristic speeds when Δr is small. In this case, if Δr < 0.9Δr_target, v < 0.9v_target, and 〈n^k∂_kv_c〉 > 0, then we replace equation (96) with equation (97), and we use $\dot{r}_{\rm drift}=v_{\rm corr}$. Here v_c is the comoving characteristic speed defined in equation (98) and the angle brackets refer to an average over the excision boundary. The condition on the gradient of v_c ensures that v_c will increase if Δr is increased; otherwise a positive $\dot{r}_{\rm drift}$ will threaten the positivity of the characteristic speeds. We continue to use $\dot{r}_{\rm drift}=v_{\rm corr}$ until either v > 0.9v_target, Δr > 0.9Δr_target, 〈n^k∂_kv_c〉 < 0, or if v will cross v_target or Δr will cross Δr_target within a time less than the current damping time τ_d (where crossing times are determined by the procedure of appendix C). When any of these conditions occur, we switch to $\dot{r}_{\rm drift}=0$, and we prohibit switching back to a positive $\dot{r}_{\rm drift}$ until either v > 0.99v_target or Δr > 0.99Δr_target, or until both v and Δr stop increasing in time. These last conditions prevent the algorithm from oscillating rapidly between positive and zero values of $\dot{r}_{\rm drift}$.

All that remains for a full specification of $\mathcal {M}_{\rm Ringdown}$ is to describe how parameters of the two discontinuous maps $\mathcal {M}_{\rm ShapeRD}$ and $\mathcal {M}_{\rm TranslationRD}$ are initialized at t = t_m. To do this, we first construct a temporary distorted frame with coordinates $\grave{x}^i$, defined by

$$\begin{equation} \bar{x}^i = \mathcal {M}_{\rm Translation} \circ \mathcal {M}_{\rm Rotation} \circ \mathcal {M}_{\rm ScalingRD} \grave{x}^i. \end{equation} \tag{ 103 }$$

The coordinates $\grave{x}^i$ are the same as the post-merger distorted coordinates except that the pre-merger translation map $\mathcal {M}_{\rm Translation}$ is used in equation (103). We then represent the common apparent horizon (which is already known in $\bar{x}^i$ coordinates) in $\grave{x}^i$ coordinates:

$$\begin{equation} \grave{x}^i_{\rm AH}(\grave{\theta }_C,\grave{\phi }_C,t) = \grave{x}^i_{\rm AH0}(t) + \grave{n}^i \sum _{\ell m} \grave{S}^C_{\ell m}(t) Y_{\ell m}(\grave{\theta }_C, \grave{\phi }_C). \end{equation} \tag{ 104 }$$

The expansion coefficients $\grave{S}^C_{\ell m}(t)$ and the horizon expansion center $\grave{x}^i_{\rm AH0}(t)$ are computed using the (inertial frame) common horizon surface plus the (time-dependent) map defined in equation (103). Here $(\grave{\theta }_C,\grave{\phi }_C)$ are polar coordinates centered about $\grave{x}^i_{\rm AH0}$, and $\grave{n}^i$ is a unit vector corresponding to the direction $(\grave{\theta }_C,\grave{\phi }_C)$. Equation (104) is the same expansion as equation (36), except here we write the expansion to explicitly include the center $\grave{x}^i_{\rm AH0}(t)$. In equation (104) we write $\grave{S}^C_{\ell m}(t)$ and $\grave{x}^i_{\rm AH0}(t)$ as functions of time because we compute them at several discrete times surrounding the matching time t_m. By finite-differencing in time, we then compute first and second time derivatives of $\grave{S}^C_{\ell m}$ and $\grave{x}^i_{\rm AH0}$ at t = t_m.

Once we have $\grave{S}^C_{\ell m}$ and its time derivatives, we initialize $\lambda ^C_{\ell m}(t_m) = -\grave{S}^C_{\ell m}(t_m)$, where $\lambda ^C_{\ell m}(t)$ are the map parameters appearing in equation (72), and the minus sign accounts for the sign difference between the definitions of equations (72) and (104). We do this for all (ℓ, m) except for ℓ = m = 0: we set $\lambda ^C_{00}(t_m)=0$, thus defining the radius of the common apparent horizon in the post-merger grid frame. For all (ℓ, m) including ℓ = m = 0 we set ${\rm d}\lambda ^C_{\ell m}/{\rm d}t\left.\right|{}_{t=t_m} = -{\rm d}\grave{S}^C_{\ell m}/{\rm d}t\left.\right|{}_{t=t_m}$, and similarly for the second derivatives.

Note that the temporary distorted coordinates $\grave{x}^i$ are incompatible with the assumptions of our control system, because the center of the excision boundary $\grave{x}^i_{\rm AH0}(t)$ in these coordinates is time-dependent. We therefore define new post-merger distorted coordinates $\hat{x}^i$ by

$$\begin{equation} \bar{x}^i = \mathcal {M}_{\rm TranslationRD} \circ \mathcal {M}_{\rm Rotation} \circ \mathcal {M}_{\rm ScalingRD} \hat{x}^i, \end{equation} \tag{ 105 }$$

where

$$\begin{equation} \hat{x}^i = \grave{x}^i-\grave{x}^i_{\rm AH0}(t). \end{equation} \tag{ 106 }$$

If we denote the map parameters in $\mathcal {M}_{\rm Translation}$ by $T_0^i$, the map parameters in $\mathcal {M}_{\rm TranslationRD}$ by Tⁱ, and if we denote $\mathcal {M}_{\rm Rotation}\circ \mathcal {M}_{\rm ScalingRD}$ by the matrix $M^i_j$, then equations (105) and (106) require

$$\begin{equation} T^i(t) = T_0^i(t) + M^i_j(t)\grave{x}^j_{\rm AH0}(t). \end{equation} \tag{ 107 }$$

Here we have assumed that f(R) appearing in equation (55) is unity in the vicinity of the horizon, so we can treat the translation map as a rigid translation near the horizon. We initialize Tⁱ and its first two time derivatives at t = t_m according to equation (107). Note that the distorted-frame horizon expansion coefficients $\grave{S}^C_{\ell m}$ and hence the initialization of $\lambda ^C_{\ell m}(t)$ are unchanged by the change in translation map because $\hat{x}^i$ and $\grave{x}^i$ differ only by an overall translation, equation (106), which leaves angles invariant.