Properties of Quasi-local mass in binary black hole mergers

Identifying a general quasi-local notion of energy-momentum and angular momentum would be an important advance in general relativity with potentially important consequences for mathematical and astrophysical studies in general relativity. In this paper we study a promising approach to this problem first proposed by Wang and Yau in 2009 based on isometric embeddings of closed surfaces in Minkowski space. We study the properties of the Wang-Yau quasi-local mass in high accuracy numerical simulations of the head-on collisions of two non-spinning black holes within full general relativity. We discuss the behavior of the Wang-Yau quasi-local mass on constant expansion surfaces and we compare its behavior with the irreducible mass. We investigate the time evolution of the Wang-Yau Quasi-local mass in numerical examples. In addition we discuss mathematical subtleties in defining the Wang-Yau mass for marginally trapped surfaces.


I. INTRODUCTION
The quasi-local definition of energy-momentum remains one of the major problems in classical general relativity [1,2]. The goal is to find appropriate notions of energy-momentum and angular momentum for finite, extended regions of spacetime. At spatial infinity and at null infinity, there are well-established concepts of energy and angular momentum. The energy-momentum defined by Arnowitt-Deser-Misner (ADM) [3] at spatial infinity measures the total energy in a spacetime, and it is conserved and shown to be positive [4,5]. The Bondi energy is measured at null-infinity and satisfies appropriate balance laws as gravitational radiation carries away energy and angular momentum [6]. Similarly, there are notions of quasi-local energy and angular momentum and balance laws applicable for black hole horizons [7][8][9]. In contrast, finding suitable analogous definitions for a finitely extended body or for an arbitrary region in spacetime is still under active research.
Finding appropriate quasi-local notions of energy momentum and angular momentum would be desirable for various reasons. For example, one might expect that gravitational waves emitted from a given region in spacetime would carry away energy thus leading to a corresponding decrease in the quasi-local mass. Such a link has been shown for the Bondi mass and for black hole horizons, but is still not available for general spacetime regions. Once fully understood, it could potentially allow us to infer detailed properties of dynamical spacetimes in * daniel.pook.kolb@aei.mpg.de † bowenzhao@bimsa.cn ‡ lars.andersson@bimsa.cn § badri.krishnan@ru.nl ¶ styau@tsinghua.edu.cn the strong field region from gravitational wave observations, such as from the merger of compact objects. On the mathematical side, it is likely that appropriate quasilocal notions of energy and angular momentum would play an important role in providing a full proof of the Penrose inequality. Similarly, as we shall discuss in this paper, quasi-local mass also plays an important role in the process of gravitational collapse and black hole formation via the hoop conjecture. We expect quasi-local mass to be a flux type integral on a closed space-like 2-surface Σ which bounds a space-like hypersurface Ω. Since Ω is not unique in the sense of being bounded by a given Σ, one would expect that a proper notion of quasi-local mass should not depend on which specific Ω is chosen. Restricting ourselves to vacuum spacetimes, we can enumerate some minimal requirements that any viable notion of quasi-local mass M (Σ) should satisfy [10]: • In flat Minkowski spacetime, M (Σ) should vanish.
• In a curved spacetime, the quasi-local mass should be non-negative.
• In the limit when Σ approaches a sphere at spacelike infinity on an asymptotically flat slice, or a cross-section of null infinity, the quasi-local mass must approach the ADM mass or the Bondi mass, respectively.
• When Σ is an apparent horizon, the quasi-local mass must be bounded from below by the irreducible mass of Σ, i.e. A Σ /16π, where A Σ is the area of Σ.
In this work we shall investigate properties of the quasilocal mass originally proposed by Wang & Yau [11]; see also [12]. There are several other proposed definitions of quasi-local mass, energy-momentum and angular momentum in the literature. Some notable ones are due to Bartnik [13], Hawking [14], and Penrose [15]; see [1,2] for a review. Based on a variational analysis of the action of General Relativity, Brown & York proposed a quasi-local energy arising as a boundary term in the Hamiltonian [16,17]; see also [10,12,18]. However, the Brown-York definition depends explicitly on a choice of spacelike hypersurface Ω that is bounded by the two-surface Σ under consideration. Specifically, the mean curvature of Σ as embedded in Ω appears in the Brown-York definition. Moreover, a specific choice of unit lapse and zero shift is needed in relating the Hamiltonian to the Brown-York mass. This rather arbitrary gauge-fixing is undesirable in general relativity studies. Furthermore, the Brown-York quasilocal mass can fail to be positive in general except for the time symmetric case [19]. On the other hand, there exist surfaces in Minkowski spacetime with strictly positive Brown-York mass. These undesirable features are resolved in Wang & Yau [20] by further including momentum information (second fundamental form in the time direction) in their definition. Indeed, Euclidean space can be regarded as the totally geodesic space-like hypersurface of zero momentum in Minkowski spacetime. While Brown & York defined their reference surface by an isometric embedding of Σ into 3-dimensional Euclidean space R 3 , Wang & Yau defined their reference surface by an isometric embedding into Minkowski space R 3,1 directly. The positivity proof of Wang-Yau quasilocal mass is given along with the definition [20]. The new definition is proven to recover the ADM mass at spatial infinity [21] and the Bondi mass at null infinity [22]. Further, the small sphere limit is proven to recover the stress-energy tensor at the limiting point for a spacetime with matter fields and is related to the Bel-Robinson tensor at higher orders for vacuum spacetime. Along the same line, they also give a quasi-local definition for angular momentum and center of mass [23], which are proven to be supertranslation invariant [24][25][26]. We will review the definition of Wang & Yau quasi-local mass below and compare with Brown & York when it is helpful.
Besides the requirements enumerated above, additional properties would be desirable when considering the dynamical aspects of general relativity. As mentioned earlier, for the Bondi mass at null infinity, the Bondi mass loss formula shows that gravitational waves carry away energy, leading to a decrease of the Bondi mass [6]. The flux of gravitational radiation is written as a surface integral over cross-sections of null infinity, and is manifestly positive. Similarly, restricting ourselves to black hole horizons and marginally trapped surfaces, similar balance laws with positive fluxes can be shown, leading to a physical process version of the area increase law [7][8][9]. Extending these considerations to a more general quasi-local setting would lead one to conjecture that the emission (or absorption) of gravitational radiation from a domain Ω could be written as a surface flux integral over Σ, directly related to the decrease (or increase) of the quasi-local mass. At present we do not have a well defined notion of such fluxes. As a first step in this direction, in this work we shall study the time evolution of Wang-Yau quasi-local mass, henceforth denoted as QLM, in the context of a binary black hole merger. This question is hard to answer analytically, and we resort instead to high precision numerical simulations of the full Einstein equations.
The plan for this paper is the following. The basics of the Wang-Yau QLM and its properties are introduced in Sec. II. We shall consider the head-on collision of two non-spinning black holes starting with time-symmetric initial data. The initial data and our numerical evolution scheme is described in III. Our numerical implementation for calculating the Wang-Yau QLM, and numerical convergence, are described in Sec. IV. The numerical results are presented in Sec. V in three steps. First, Sec. V A shows the results in the initial data, i.e. with time symmetry. As the evolution proceeds, the later time slices are no longer time-symmetric. Sec. V B shows results for non-time-symmetric slices and finally Sec. V C presents the time evolution of the QLM and also an exploration of the hoop conjecture in the context of the formation of the common horizon in a black hole merger. In the course of presenting the numerical results, it will be clear that there are mathematical subtleties in defining the QLM for a marginally trapped surface. This will be clarified mathematically in Sec. VI and will justify the various choices made in the numerical work. Finally, Sec. VII will discuss some implications of our results and suggestions for future work.

II. BASIC NOTIONS
The Wang-Yau quasi-local energy (QLE) associated with a suitable surface is defined through anchoring the surface intrinsic geometry while comparing the extrinsic geometry as embedded in the original spacetime N 4 versus that embedded in the flat Minkowski space R 3,1 . Given a spacelike two-surface Σ ⊂ N 4 with induced metric σ ab , let i 0 : Σ → R 3,1 be an isometric embedding into the Minkowski spacetime. Fixing a unit, future-pointing, timelike vector T 0 in R 3,1 , a one-to-one correspondence between vector fields in N 4 and those in R 3,1 is built through the 'canonical gauge' condition, where H and H 0 are mean curvature vectors of Σ ⊂ N 4 and i 0 (Σ) ⊂ R 3,1 , respectively, and ⟨ · , · ⟩ denotes the corresponding scalar product in N 4 and R 3,1 . The basis vectorsē α of N 4 , α ∈ {0, 1, 2, 3}, andě α of R 3,1 are chosen as follows. Letě 3 be the spacelike unit normal of i 0 (Σ) which is also perpendicular to T 0 . Letě 4 be the future-pointing, timelike unit normal that is perpendicular toě 3 . Then {ě 3 ,ě 4 } forms an orthonormal basis for the normal bundle of i 0 (Σ) ⊂ R 3,1 . The canonical gauge condition (1) picks uniquely a future-pointing, timelike unit normal of Σ,ē 4 . Thenē 3 is the spacelike normal of Σ that combined withē 4 gives an orthonormal basis for the normal bundle of Σ ⊂ N 4 . Given τ ∈ C ∞ (Σ), a generalized mean curvature for Σ is defined as where ∇ denotes the covariant derivative on Σ associated with σ ab , |∇τ | 2 = σ ab ∇ a τ ∇ b τ and we write αě 3 (∇τ ) = (αě 3 ) a ∇ a τ . The connection one-form αē 3 associated with the basis {ē 3 ,ē 4 } is defined as where Y ∈ T Σ and (4) ∇ denotes the covariant derivative in N 4 . Similarly, one can define αě 3 for the connection one-form associated with {ě 3 ,ě 4 } in R 3,1 as where (3,1) ∇ Y denotes the covariant derivative in R 3,1 . A generalized mean curvature for i 0 (Σ) is defined as The Wang-Yau quasi-local energy associated with τ is then defined as When the mean curvature vector H is spacelike, one can use H = −⟨H,ē 4 ⟩ē 4 + ⟨H,ē 3 ⟩ē 3 = pē 4 − kē 3 and its conjugate vector J = kē 4 − pē 3 to form an orthonormal basis for the normal bundle N Σ, {e H = − H |H| , e J = J |H| }. In terms of this mean curvature vector basis, and similarly for θ 0 in R 3,1 .
Solving the variational problem of minimizing the QLE with respect to the time function τ , one gets the Euler-Lagrange equation, called the optimal embedding equation (OEE), The minimum value of QLE is defined to be the Wang-Yau quasi-local mass where (see (4.4)-(4.5) in [27] for details) and The Wang-Yau QLM is defined for any closed spacelike surface Σ whose mean curvature vector is spacelike and where an admissible solution to the OEE (6) exists (see Definition 5.1 in [20] for admissible τ ). Note that if τ = const is admissible and solves the optimal embedding equation, it must be the global minimum of Wang-Yau quasi-local energy [28]. Substituting τ = const to (7), one sees that the Wang-Yau quasi-local mass reduces to the Liu-Yau mass in this case [10,12] If further Σ lies in a totally geodesic slice, the Wang-Yau quasi-local mass reduces to the Brown-York mass where k is the only nonzero component of H lying in the totally geodesic slice while k 0 is the mean curvature vector of i 0 (Σ) embedded in R 3 . Note that τ only appears through derivatives and we hence use τ = 0 and τ = const interchangeably. In black hole spacetimes, there is a particular set of surfaces of interest-the marginally outer trapped surfaces (MOTSs). These are used to study various aspects of black holes quasi-locally via the framework of isolated and dynamical horizons, respectively describing black holes in equilibrium and in dynamical situations (see e.g. [29][30][31][32][33]). In any given spacelike slice S ⊂ N 4 , the outermost MOTS is called the apparent horizon. Given a MOTS on a spacelike slice, the results of Andersson et al. [34][35][36] show the conditions under which it evolves smoothly. It is shown that when a MOTS is stable under outward deformations, then it will evolve smoothly. Recent numerical work has applied and further explored the stability of MOTSs and its implications for the time evolution [37][38][39][40][41][42].
A MOTS is defined as follows. Let Σ ⊂ N 4 be a closed spacelike surface and let ℓ + , ℓ − be two future directed null normal fields on Σ taken to point outward and inward, respectively. We fix the cross normalization by ⟨ℓ + , ℓ − ⟩ = −2, which still leaves a remaining freedom to scale for any positive function f . The outward and inward null expansions, denoted Θ + and Θ − respectively, are defined as where σ αβ = σ ab π α a π β b with π α a the projection onto the tangent bundle of Σ. Then, Σ is called a marginally outer trapped surface (MOTS) if Θ + = 0, an outer trapped surface if Θ + < 0 and an outer untrapped surface if Θ + > 0. A marginally trapped surface (MTS) is a MOTS with Θ − < 0. Note that although Θ + → f Θ + under (12), the signs of Θ ± are invariant and so these definitions are not affected.
We will often consider families of surfaces Σ s with constant Θ + = s ∈ R, called constant expansion surfaces (CESs). These do depend on the choice of ℓ ± , which we fix uniquely using the spacelike slice S ⊂ N 4 within which the family Σ s is constructed. Concretely, let Σ s ⊂ S, v the spacelike outward unit normal of Σ s in S and u the future timelike unit normal on S in N 4 . Then, we choose In terms of the null expansions Θ ± , the mean curvature vector H and its conjugate vector J can be expressed as Since ⟨H, H⟩ = −Θ + Θ − , the mean curvature vector becomes a null vector on a MOTS. If in addition the slice S is time symmetric, i.e. its second fundamental form vanishes, then σ αβ (4) ∇ α u β = 0 and thus Θ + = −Θ − implying that H and J both vanish on a MOTS. To characterize the quasi-local mass of a black hole region, we take the QLM of apparent horizons. However, the current definition of the Wang-Yau quasi-local mass assumes the surface mean curvature vector to be spacelike and hence does not apply to MOTSs. Therefore, one of our goals is to extend the definition of the Wang-Yau QLM (7) to a MOTS in time symmetry and to an MTS without time symmetry. Limiting ourselves to the case of axisymmetry and no angular momentum, we will argue that a suitable extension is With this extension, we can then investigate the time evolution of QLM during black hole collisions. As noted above if Σ lies in a totally geodesic slice, e.g. in the moment of time-symmetry, Wang-Yau QLM reduces to Brown-York mass. In this case, the above extension simply reduces to QLM = 1 8π k 0 , which is what one would expect for the Brown-York limit at minimal surfaces. Further extension of QLM to surfaces of timelike mean curvature vector, e.g. trapped surfaces inside event horizons, is certainly of great interest and will be studied elsewhere.

III. INITIAL DATA AND NUMERICAL EVOLUTION
We use Brill-Lindquist initial data [43], which solves the constraint equations of General Relativity with vanishing extrinsic curvature and vanishing scalar curvature, i.e. a time-symmetric slice in vacuum spacetime. The Riemannian three-metric is defined on R 3 \ {x 1 , . . . , x n } with n + 1 asymptotically flat ends, one at ∥x∥ → ∞ and n at the punctures x i . We restrict ourselves to the case n = 2, which describes a two-black-hole configuration. The three-metric can then be written as where δ ij is the flat metric and the conformal factor is We take the two punctures to be located on the z-axis at coordinates x A,B = (0, 0, ±d/2), respectively. The three ends at ∥x∥ → ∞, x A and x B , respectively, have ADMmasses For sufficiently large d, the slice S contains two separate black holes, each surrounded by a stable MOTS that contains either x A or x B . We shall call these the individual MOTSs Σ A and Σ B , respectively. If d becomes small enough, there exists a stable common MOTS Σ outer surrounding Σ A,B . In fact, as d passes through the value at which Σ outer appears, it is found that an unstable MOTS Σ inner forms together with Σ outer and "moves" inward as d is decreased. This is discussed in more detail elsewhere [44]. 1 The two common and two individual MOTSs for an equal mass configuration are shown in Fig. 1.
The numerical data for this initial slice are generated by the TwoPunctures [45] thorn of the Einstein Toolkit [46,47]. These data are evolved in time using an axisymmetric version of McLachlan [48], which in turn uses Kranc [49,50] for generating C++ code. This uses the BSSN formulation of the Einstein equations with gauge conditions chosen as the so-called 1 + log slicing and a Γdriver shift condition [51,52]. More details about our numerical simulation setup, including a convergence analysis, are described in [38].
Our analysis is based on two simulations, both starting from BL data. The first, referred to as Sim1, uses initial data with m B /m A = 2, d = 0.9 and the second, simulation Sim2, uses m B /m A = 1.6, d = 1. Both simulations were performed with different spatial grid resolutions to check the accuracy of our calculations. Results shown for Sim1 use a resolution 1/∆x = 720, which was evolved until simulation time t f = 6. For Sim2, we used a lower resolution of 1/∆x = 312 to extend the evolution up to time t f = 38.
The MOTSs and CESs are numerically found with high accuracy both in the analytical initial data as well as in slices produced by the Einstein Toolkit using the method in [44,53].
In general, the problem of locating a surface Σ s with expansion Θ + = s may have many solutions within a given slice S. For s = 0, this corresponds to the different MOTSs in S. By choosing suitable initial guesses for the numerical search, we can easily select which particular MOTS to find. As mentioned above, we focus here on the three stable MOTSs Σ outer , Σ A and Σ B , interpreted as the horizon of the merger remnant and the smaller and larger (in case of unequal masses) individual black holes. Choosing one of these MOTSs as initial guess, we construct CESs for s close to zero. Families of Σ s are then built by taking small steps in s, each time using the previous CES as initial guess for the next. CESs far from Σ outer in the nearly flat region of S are close to being spherical in our coordinates and so we can use coordinate spheres as initial guesses in this case.

IV. NUMERICAL METHOD FOR EVALUATING THE QUASI-LOCAL MASS
The strategy to solve the optimal embedding equation ∇ a j a = 0 is to consider ∇ a j a as a nonlinear operator L acting on τ , linearize that operator L and solve the linear problem multiple times, each time taking a small step towards a solution of the full nonlinear problem. This is also called the Newton-Kantorovich method [54,Appendix C].
Analytically linearizing L requires determining the explicit dependency of ∇ a j a on τ . We make use of axisymmetry, i.e. assuming the surface Σ, the embedding i 0 and τ are all axisymmetric, to simplify calculations. The time function τ then depends only on one parameter, say θ, increasing from one pole of Σ to the other, and the embedding i 0 can be expressed in terms of the intrinsic metric σ ab and τ ′ , where τ ′ = dτ dθ . Explicitly, for coordinates {y 1 = θ, y 2 = ϕ} on Σ, 0 < θ < π and 0 < ϕ < 2π, an axisymmetric ansatz for the embedding is Writing the two-metric as dσ 2 = P 2 dθ 2 + Q 2 sin 2 θ dϕ 2 , we then have where V = Q sin θ. Using this, we calculate ⟨H 0 , Furthermore, The respective terms in curved space are calculated differently using the null expansions, which we have in highly accurate form from the MOTS and CES finding process. In addition to (15), we use where Θ ± ,a = ∂Θ± ∂y a . We remark that Θ ± contains first derivatives of the 3-metric, which means that ∇ a j a contains third derivatives. In order to get numerically accurate results, we expand Θ ± into a set of basis functions and differentiate these directly. To linearize the operator Lτ = ∇ a j a , the above expressions are first inserted in turn into (8), (9) and ∇ a j a . Afterwards, we use SymPy [55] to symbolically differentiate L with respect to τ ′ , τ ′′ , τ (3) and τ (4) . In the end, the linearized operator we implement into our numerical code is of the form where ∆ is a scalar function on Σ. Starting with an initial guess τ 0 , usually τ 0 = 0, we perform steps τ i+1 = τ i + ∆ i , where ∆ i solves the linear equation which we solve using a pseudospectral method. In most cases, it took between 5 and 15 steps to converge up to numerical roundoff. Fig. 2 shows that the residual of the OEE (6) decreases exponentially with the resolution of τ , where the resolution is the number of basis functions used for the finite representation of τ .

A. The QLM in time-symmetric initial data
For a time-symmetric slice S, the mean curvature vector H of Σ lies in S. A MOTS therefore coincides with a minimal surface (k = −⟨H, v⟩ = 0, where, as before, v is the outward unit normal of Σ in S). Moreover, in time symmetry we have α H = −⟨ (4) ∇e J , e H ⟩ ≡ 0 since e H lies in S while e J is the normal to the totally geodesic slice S. And for τ = const, i 0 (Σ) ⊂ R 3 , by a similar argument, α H0 ≡ 0. Thus τ = const trivially solves the optimal embedding equation (6). Furthermore, τ = const is known to be the global minimum of the QLE provided that it solves the optimal embedding equation [28]. The Wang-Yau quasi-local mass reduces to the Brown-York mass m BY for any surface Σ in a moment of time symmetry.

On the monotonicity along geometric flows
It is well known that for some cases, the Brown-York mass exhibits a monotonically decreasing behavior [56]. As an example, consider a Schwarzschild black hole. In a time-symmetric slice, with the metric in isotropic coordinates the horizon lies at r = m/2. One can show that This is interpreted as negative gravitational field energy bringing down m BY as the surface approaches infinity [2]. In fact, this monotonicity property could be shown more precisely. Consider two mean convex surfaces in S, Σ i = ∂Ω i , i ∈ {1, 2}, with Ω 1 ⊂ Ω 2 and suppose there exists a geometric flow where B and B 0 are the second fundamental forms of Σ t ⊂ S and of i 0 (Σ t ) ⊂ R 3 , respectively. For a timesymmetric slice S, the scalar curvature R = 2T 00 and hence the R term can be interpreted as matter contribution, which in our case vanishes. The remaining |B 0 − B| 2 − (k 0 − k) 2 term is then supposed to characterize the pure gravitational field energy. Note that although the integrand |B 0 − B| 2 − (k 0 − k) 2 clearly depends on the foliation Σ t , the total integral does not. The work of Huisken & Yau [58] and later improvements [59][60][61][62][63] show that ends of an asymptotically flat Riemannian 3-manifold with positive ADM mass and nonnegative scalar curvature admit a unique canonical foliation through stable constant mean curvature (CMC) surfaces. The geometric flow is assured in this case with Σ t being CMC surfaces. Nonetheless, assume (31) holds true (in our case true numerically, see below), one can discuss the sign of the gravitation field energy term. For simplicity, take an orthonormal basis for Σ t , {e 1 , e 2 }, such that σ ab = δ ab . This basis is also an orthonormal basis for the isometric embedding i 0 (Σ t ). Then in this basis, where as a linear map in this ON basis, with D R 3 and D denoting covariant differentiation in R 3 and in Riemannian S, respectively. If B 0 − B can be chosen to be orientation preserving throughout a foliation Σ t , then it indicates negative gravitational field energy. Then m BY would be monotonically decreasing as the surfaces approach infinity. We suspect this is generally true for at least mean convex Σ but a proof is missing at the moment.
We examine the above equality (31) numerically. Fig. 3 shows the QLM for CESs with Θ + ≥ 0. Outside the apparent horizon, the QLM decreases monotonically with increasing distance to the horizon, whereas the expansion Θ + = k increases from 0 at first (right panel, upper part) and then drops back to 0 at infinity (right panel, lower part). The balance (31) is also verified (see Fig. 4) although we cannot prove for now that there exists a geometric flow among constant expansion (mean curvature) surfaces in our case.

Outer trapped surfaces in time symmetry
In time symmetry, we have Θ + = −Θ − and, since k < 0 everywhere for outer trapped surfaces (Θ + < 0). However, both the Wang-Yau QLE and the Brown-York mass implicitly make the assumption that k > 0 and hence do not apply to such surfaces. We argue that for k = −|H| < 0 in the time-symmetric case, the same formula works with k replaced by its absolute value. We first exemplify this with a time-symmetric slice of the Schwarzschild metric (30). Recall that there is an isometric inversion r → (m/2) 2 r that sends surfaces inside the horizon to the outside, reversing the sign of H while   Fig. 3, but continuing the families to the inside of the apparent horizon in both the BL data (blue) and in a Schwarzschild slice (orange). We show a comparison of the Brown-York mass (11) calculated either via 1 8π Σ k0 − |k| (solid lines) or via 1 8π Σ k0 − k (dotted lines). For Θ+ ≥ 0, the two definitions agree. In the Brill-Lindquist case, the CES family interpolates between the inner common MOTS Σinner and the apparent horizon Σouter via surfaces with Θ+ < 0. However, these latter surfaces intersect each other (Fig. 7). Monotonicity in this regime can therefore not be expected. See text for discussion. ted), and compared with the analogous Brill-Lindquist case (blue). Note that r denotes the isotropic radial coordinate. The QLM attains its maximum value 2m at the horizon while both the r → ∞ and r → 0 limit yield QLM = m, consistent with the interpretation of a time-symmetric slice in Schwarzschild as a wormhole connecting two identical, asymptotically flat regions. More specifically, for Θ + < 0, the surface lies in the other asymptotic region, where its outward normal is −v. If we use k < 0 as-is, then the BY-mass remains monotonic and diverges as ∥x − x i ∥ → 0 (i.e. approaching the other end x i ). In summary, a naive extension of Brown-York mass into the k < 0 region results in a smooth QLM profile while taking absolute value yields a "kink" at the horizon. In this case, a non-smooth QLM profile (k 0 −|k|) is clearly more natural than a smooth one (k 0 − k) and one might in general expect some non-analytic behavior of QLM around MOTSs or horizons. For multiple black holes, taking the absolute value of k yields results consistent with [43]. That is, for large spheres one recovers the ADM mass as expected while for small spheres approaching each puncture x i , an asymptotic expansion yields the ADM mass associated to each puncture (19). Numerical calculation for CES around Σ outer is also included in Fig. 5 (blue) and reveals a similar behavior as in the Schwarzschild case up to Θ + ≈ −0.2 where a turn-over happens in the multiple holes BL data (see below for discussion). Numerical calculation for CES toward each puncture, i.e. around Σ A,B , is shown in Fig. 6 and again reveals a similar behavior as in the Schwarzschild case.
We emphasize that in Brill-Lindquist data, only CES families outside Σ outer and inside each individual Σ A,B are comparable with the Schwarzschild case. The region bounded by Σ outer surrounding all punctures and Σ A,B surrounding each individual puncture does not have a direct correspondence in the Schwarzschild case. The CES family interpolates between Σ outer and the unstable common MOTS Σ inner and does not approach either of the two asymptotic ends x A,B (Fig. 7). Moreover, in this region, constant expansion or constant mean curvature (CMC) surfaces fail to foliate space and monotonicity of the QLM indeed fails here. These explain peculiar features around Σ inner in BL data seen in Fig. 5. We remark that choosing other families of foliating surfaces such as coordinate spheres yields qualitatively similar conclusions. The choice of families of CESs, which in time symmetry are constant mean curvature surfaces, allows us to get arbitrarily close to a MOTS, which is a CES itself.
Our extension confirms that QLM at the MOTS (16) 5 10 r(Σ 2 ) should be where |Σ| denotes the area of Σ and the Minkowski inequality is invoked. This is already assumed in earlier studies [2].

B. QLM in non-time-symmetric slices
As one numerically evolves the time-symmetric initial data, the slices S become non-time-symmetric and the mean curvature vector of a surface Σ ⊂ S may acquire a timelike component. The Wang-Yau quasi-local mass will then in general differ from the Brown-York mass.
This has various consequences, one being that the monotonicity of the QLM along geometric flows is not guaranteed by (31) anymore, although we numerically find monotonicity remains true, as can be seen in Fig. 8. An analytic generalization of (31) to the non-timesymmetric case is under study.

QLM at a MOTS without time symmetry
In this section, we will numerically determine the QLM on CES families near Σ outer . The goal is to justify the definition (16) of the QLM on a MOTS Σ by exploring its behavior as we approach Σ along the family from the Θ + > 0 side, where the QLM is well defined.
Although ρ (8) and QLM (7) seem well-defined even for |H| → 0, the OEE that determines τ is clearly singular at a MOTS: |H| appears as denominator in both (9) and (27). Therefore, the definition of the QLM cannot trivially be extended to the case of a MOTS. In other words, as Θ + → 0, |H| → 0, and so the assumption of having a spacelike mean curvature vector in the Wang-Yau QLM breaks down. We focus on examining this issue here numerically. A mathematical treatment of this case will be given in Sec. VI.
To numerically explore what happens as Θ + → 0, we look at the individual terms in the OEE (9) that determines τ . As can be seen in Fig. 9, ρ remains finite while τ θ = τ ′ approaches zero, so their product vanishes at a MOTS. Since τ ′ → 0, τ → const and clearly α H0 → 0. The remaining terms ψ θ and α H both remain finite in the limit. However, within numerical limits, they cancel in j θ as Θ + → 0. The net result is thus τ → 0 is a solution to the OEE at a MOTS.
This suggests that as we approach a MOTS, τ becomes constant and j a vanishes, so that ρ → |H 0 |. In terms of the QLM, this limit is shown in Fig. 10 together with the value of the QLM calculated at the MOTS using (16).

C. Time evolution of the QLM at a MOTS
Having argued that one can extend the Wang-Yau QLM to the common apparent horizon and each individual horizon by (16), we now examine its time evolution during the head-on merger of two black holes.
If one interprets the Wang-Yau QLM as a measure to separate quasi-local degrees of freedom from traveling gravitational waves, then Fig. 11 indicates the region surrounded by Σ outer loses energy/mass while subregions surrounded by each individual horizon Σ A,B gain energy/mass during the collision. Furthermore, in the longer simulation Sim2, we find an oscillation in energy/mass contained inside the apparent horizon Σ outer (Fig. 12). It is well known that for two black hole collisions, the intrinsic geometry of the apparent horizon experiences oscillations at the (lowest) quasi-normal frequency of the final black hole [64]. While the apparent horizon area monotonically increases despite the oscillation at the quasi-normal frequency, QLM of the apparent horizon fails to maintain monotonicity with the oscillation. Nevertheless, as the final black hole settles down to equilibrium, the measure of QLM and area for the apparent horizon converges.
First of all, although the Wang-Yau quasi-local mass is defined for a 2-surface and does not depend on the choice of slicing S, determining the apparent horizon through Θ + does depend on the choice of slicing S. Bearing in mind this ambiguity associated with the apparent horizon, one tends to conclude that the Wang-Yau quasi-local mass suggests that the region bounded by the outermost MOTS Σ outer could lose energy to infinity. This is a different picture from that indicated by the area of the horizons, which increase monotonically, and the standard first law of black hole thermodynamics. This difference is actually expected: equation (6.20) in [16] indicates that the Brown-York quasi-local mass for a Schwarzschild black hole satisfies a balance equation that differs from the standard black hole thermodynamics. This again emphasizes that quasi-local mass as defined by Brown and York or Wang and Yau measures a different quantity than irreducible mass. More specifically, as the horizon expands, more negative gravitational field energy is taken into account by QLM which counteracts the positive contribution from growing area and absorbed gravitational waves. This issue is crucial in understanding Wang-Yau quasi-local mass and will be carefully studied with a direct calculation of gravitational wave energy in a future study.

Investigating the hoop-conjecture
A viable notion of QLM in general relativity should lead to numerous applications. We conclude this section on numerical results by exploring one such application, namely the hoop conjecture [65,66]. This conjecture addresses the question of under what conditions a black hole forms. As we shall shortly see, several aspects of the conjecture are not precisely formulated and numerical relativity has a role to play in exploring and evaluating various possibilities; see e.g. [67] for recent work in this direction.
In the case of gravitational collapse, the intuitive pic- The smaller horizon ΣA (not shown) has a qualitatively similar behavior, though less pronounced, as ΣB. For easier comparison, the QLM has been divided by 2 to account for the fact that for Schwarzschild, the QLM of the horizon is twice the ADM mass.
ture is that of matter fields getting compressed due to their self-gravity and eventually becoming sufficiently dense to form a black hole. As originally formulated by Kip Thorne: Horizons form when and only when a mass M gets compacted into a region whose circumference situations.
If a notion of QLM is generally viable, it should be possible to use it as the appropriate mass within the hoop conjecture. For the Schwarzschild spacetime, as we have seen, the Wang-Yau QLM for the horizon is twice the ADM mass. Thus, one might expect the relevant hoop conjecture inequality should be modified to C/4πM < 0.5. For the hoops, we shall use closed geodesics lying within the constant expansion surfaces that we have already found. In our present case, we do not deal with gravitational collapse, but rather with a binary black hole merger where we always have black holes present on any time slice. We instead seek to investigate the issue of when the common horizon forms, and whether its formation can be predicted by a hoop conjecture argument using the Wang-Yau QLM. We assume further that the hoop conjecture applies to the formation of marginally trapped surfaces. In our case, the constant expansion surfaces and the marginally trapped surfaces turn out to be prolate so that the polar circumference C p is larger than the equatorial circumference. Thus we calculate the ratio C p /4πM for constant expansion surfaces on different time slices in the vicinity of the time when the common horizon is first formed.
Our results are shown in Fig. 13. We see that the ratio C p /4πM approaches 0.5 asymptotically as expected. At earlier times, this ratio is somewhat larger than the limiting value 0.5. This is not unusual -similar results were found in e.g. [67] where the ratio was about 12% larger than the limiting value predicted by the hoop conjecture. At each time just before the horizon is formed, the value of C p /4πM approaches the limiting value as the expan- sion becomes smaller. These results are suggestive but inconclusive -it is not yet clear whether this can be used as a prediction for the formation of the common horizon. This would require us to identify a suitable threshold for C p /4πM for these surfaces. In the above discussion we have considered only the constant expansion surfaces. It is plausible that there exist other surfaces which have a smaller value of C p /4πM . Therefore, while we shall not do so here, it would be more satisfactory to consider a more general class of 2-spheres and to minimize the ratio C p /4πM over these spheres. Following the spirit of the hoop conjecture, one could then look for a threshold on this minimum value of C p /4πM , and investigate whether it can predict the formation of a black hole horizon.

VI. DEFINING THE QLM ON A MOTS
The above numerical results have already used a "working definition" (16) for evaluating the Wang-Yau QLM on a MOTS. In this section, we explore the limit Θ + → 0 from a mathematical perspective to justify this definition. It seems plausible that an extension to the case with angular momentum is possible. However, this will be left to future work.
Lemma VI.1. Consider an axisymmetric collision with no angular momentum involved. Further assume isometric embedding or time function τ being axisymmetric, then j ≡ 0 when it is well-defined, i.e. when the mean curvature vector H is spacelike.
Proof. Use the definition (9). We first observe that under the above assumptions, j ϕ ≡ 0. It is easy to see that ∇ ∂ϕ τ = ∇ ϕ sinh −1 ρ∆τ |H0||H| = 0 by assumption. That connection one-forms vanish can be seen from computing Christoffel symbols, using subscripts 3, 4 for the H and J direction, and noting that the spacetime under consideration possesses a symmetry ϕ → −ϕ such that cross terms in the metric involving ϕ all vanish. Next we show j θ = 0. The most general axisymmetric metric for a 2-surface is Using that Q, P do not vanish (metric non-degenerate), sin(0) = sin(π) = 0 and that j is well-defined, one gets const = 0 and hence j θ ≡ 0 on the surface.
Remark VI.2. This may sound contradictory to the well-known fact that there are infinitely many nonvanishing, smooth, divergence-free vector fields on S 2 . We elaborate on this. Denote the volume form by ω = √ σdθ ∧ dϕ. Take divergence-free vector field j = j θ ∂ θ with j ϕ = 0. Then div(j)ω = L j ω = dι j ω = 0 .
Given that H 1 dR (S 2 ) = 0, there exist a smooth function f (θ, ϕ) such that ι j ω = df . That is, It follows that ∂f ∂θ = 0. Then noting that the RHS is only a function of ϕ while the LHS clearly has a dependence on θ, it follows that j θ = f = 0. We again reach j ≡ 0 when axisymmetry is imposed.
Remark VI.3. When angular momentum is present, the metric would have cross terms involving ϕ and the connection one-form α H would generally not vanish (α H0 = 0 since the reference spacetime is static). In fact, the quasi-local angular momentum as defined in [23] vanishes exactly when j ϕ = 0, consistent with our proof above. Recall that quasi-local angular momentum is defined as If one assumes Σ and τ both axisymmetric, ⟨K, T 0 ⟩ = 0. E(Σ, X, T 0 , K = ∂ ϕ ) = 0 if and only if j ϕ = 0.
Theorem VI.4. Assume H 0 remains spacelike as H turns into null or |H| → 0, then the solution τ to the OEE approaches a constant as the surface approaches the apparent horizon.
First look at the α H term. Use (27) and the fact that the family of Σ is constant expansion surface with ∂ a Θ + ≡ 0, and hence bounded as θ + → 0. Next look at the α H0 term.
To gain more understanding about the limit τ → const, we use another formula for j that does not invoke picking the specific frame {e H , e J } j = ρ∇τ − αě 3 + α e3 .
Remark VI.5. The above proposition shows that one can extrapolate the definition of the QLM to |H| = 0 with the optimal embedding being τ = const. In this case, i.e. τ = const and |H| = 0, Remark VI.6. As proved in [28], a solution τ to the optimal embedding equation is a local minimum of Wang-Yau quasi-local energy if where H τ is the mean curvature vector of the isometric embedding with time function τ . The condition is satisfied for constant expansion surfaces with Θ + > 0 (Fig. 14). Assume continuity, τ = const is a local minimum of Wang-Yau quasi-local energy at the apparent horizon.
The numerical calculations in Sec. V B 1 comply with the above results.

VII. CONCLUSIONS
In this work, we have studied the Wang-Yau quasilocal mass during a binary black hole merger. The Wang-Yau QLM uses an embedding of 2-surfaces in Minkowski space. We have solved the optimal embedding equation numerically and applied it to the head-on collision of two non-spinning black holes starting with Brill-Lindquist initial data. We numerically determined the QLM for surfaces close to the horizons and for families of surfaces approaching the asymptotically flat ends and study their time evolution, and also present a preliminary investigation of the hoop conjecture applied to the formation of the common horizon. Finally, we have suggested an extension of the Wang-Yau QLM to marginally trapped surfaces. For a Schwarzschild black hole, our calculations agree with the well known result of the Brown-York mass, i.e. the QLM is twice the ADM mass on the horizon. The Brown-York mass decreases monotonically as one moves outward from the horizon, and for the sphere at infinity, it yields the ADM mass. This is in sharp contrast to other quasi-local mass definitions such as Hawking and Bartnik masses. The Wang-Yau quasi-local mass also inherits this monotonic decreasing property. Such monotonicity is clearly demonstrated numerically for our family of constant expansion surfaces. Therefore, it is a crucial property of the Wang-Yau quasi-local mass that some measure of negative gravitational field energy is accounted for. In particular, for surfaces Σ in a time-symmetric slice S, an explicit expression for gravitational field energy was written down [57]. An analogous expression for the case of non-time-symmetric slices is under study.
We have extended the definition of Wang-Yau quasilocal mass for 2-surfaces of spacelike mean curvature vector to 2-surfaces of null mean curvature vector, i.e. (16). With this extension, we examined the time evolution of QLM for a black hole, defined as QLM at Σ outer , during the head-on collision of two non-spinning black holes. As is well known, the area increases monotonically throughout the evolution. At late stages, the area evolution exhibits damped oscillations which are known to be associated with the quasi-normal modes of the remnant black hole (see e.g. [68]). In contrast, QLM decreases monotonically at first and starts to lose monotonicity when oscillations take place. We see from the bottom panel of Fig. 12 that the oscillation frequency of the QLM is similar to that of the area and one might expect these to also be associated with the quasi-normal modes. Monotonicity could be included as one of the additional requirements for a QLM, and in future work we shall investigate the possibility to modifying the definition of the QLM appropriately to make it monotonic.
That QLM and area evolve differently seems to comply with a variation equation for the Brown-York mass in the Schwarzschild black hole case [16], which differs from the standard first law of black hole thermodynamics. One might again invoke negative gravitational field energy to explain this difference. Nevertheless, both of these measures, namely the area and half of QLM of the apparent horizon converge to the same value asymptotically as the final black hole settles down to its equilibrium state. We expect that employing estimates of gravitational wave energy will greatly clarify on various distinctive features of the Wang-Yau quasi-local energy.
Future work will extend this study in various directions. Further extension of quasi-local definition for mass and angular momentum to surfaces of time-like mean curvature vector is of great interest. An example is trapped surfaces inside the horizon, whose quasi-local mass might reveal important information about black holes. An attempt definition based on Brown-York mass was made in a previous study [17], where it is found the quasi-local mass either experiences an infinite slope or a cusp at the horizon while the former was preferred by those authors. A naive extension of Wang-Yau QLM presented here reveals a similar behavior. A careful study of this issue will be presented elsewhere.
Furthermore, it will be important to calculate the QLM in cases with rotation or angular momentum. In these cases, we can also calculate the quasi-local angular momentum. Quasiloal mass and angular momentum for surfaces inside Kerr's ergoregion might exhibit interesting features related to the Penrose process. The time variation of the mass and angular momentum should be related to appropriate fluxes. Moreover, the event horizon, whose determination requires the knowledge of the whole spacetime evolution, may reveal different information from the apparent horizon studied here. Finally, an extension of the QLM and angular momentum to cosmological spacetimes would be of great interest. Here it would be necessary to consider a reference configuration not in Minkowski spacetime, but in de Sitter or antide Sitter spacetime [69].