Revisiting the apparent horizon finding problem with multigrid methods
Hon-Ka Hui Lap-Ming Lin



















Department of Physics, The Chinese University of Hong Kong, Hong Kong SAR, China
hkhui@link.cuhk.edu.hk, lmlin@cuhk.edu.hk







     Apparent horizon plays an important role in numerical relativity as it provides a 
     tool to characterize the existence and properties of black holes on three-dimensional spatial  slices in 3+1 numerical spacetimes. Apparent horizon finders based on different techniques have been developed. 
     In this paper, we revisit the apparent horizon finding problem in numerical relativity using multigrid-based algorithms. 
 We formulate the nonlinear elliptic apparent horizon equation as a linear Poisson-type equation with a nonlinear source, and solve it using a multigrid algorithm with Gauss-Seidel line relaxation.
 A fourth order compact finite difference scheme in spherical coordinates is derived and employed to reduce the complexity of the line relaxation operator to a tri-diagonal matrix inversion. 
    The multigrid-based apparent horizon finder developed in this work is capable of locating apparent horizons in generic spatial hypersurfaces without any symmetries.
    The finder is tested with both analytic data, such as Brill-Lindquist multiple black hole data, and numerical data, including off-centered Kerr-Schild data and dynamical inspiraling binary black hole data. 
    The obtained results are compared with those generated by the current fastest finder AHFinderDirect (Thornburg, Class. Quantum Grav. 21, 743, 2003), which is the default finder in the open source code Einstein Toolkit. 
Our finder performs comparatively in terms of accuracy, and starts to outperform AHFinderDirect
at high angular resolutions (∼ 1^∘) in terms of speed. Our finder is also more flexible to initial guess, as opposed to the Newton's method used in AHFinderDirect. This suggests that the multigrid algorithm provides an alternative option for studying apparent horizons, especially when high resolutions are needed. 



Keywords: apparent horizon, multigrid, compact finite differences, Poisson equations






§ INTRODUCTION


    Apparent horizon has been an important concept for the numerical studies of spacetime that involve black holes. It serves as a real-time indicator of the existence of a black hole in numerical simulations, as its presence implies the presence of a surrounding event horizon <cit.>. The determination of apparent horizon during numerical simulations is made possible because of its local nature. This is in contrast to the teleological nature of event horizon, whose determination requires the full knowledge of the future evolution of a spacetime. Knowing the presence of an apparent horizon during simulation is important not only because it allows the estimation of the properties of the black hole, such as its mass and spin, but it also marks the region within which singularity may appear at a later time <cit.>, such region sometimes is necessary to be avoided through excision to maintain stability of the simulation <cit.> (also see, e.g., <cit.> for using singularity-avoiding gauge instead of excision).

    The dynamical evolution of the apparent horizon is also of interest for understanding black holes and gravitational wave physics. The notion of dynamical horizons (e.g., <cit.>, also see <cit.> for a review), which is loosely speaking 
    a world tube formed by the foliation of the apparent horizon at different time slices, can be used instead of event horizons to better study the interaction and merging of black holes <cit.>. Recent studies have also suggested correlation between the dynamics of (common) apparent horizon in black hole binaries with the gravitational waves emitted <cit.>, which could hint a possible imprint of event horizon properties on the signal.
    
    
    
    

    Numerous efforts have been dedicated to developing an efficient apparent horizon finder, as discussed in a comprehensive review <cit.>. 
    Currently, the fastest algorithm for locating the apparent horizon is based on Newton's method (<cit.>), along with some variants of the matrix inversion scheme for the associated Jacobi equation. In general, such a method suffers from the poor scaling of the computational time complexity and the reliance on a good initial guess. On the other hand, while multigrid methods in general exhibit much better scaling, they have not been applied to the apparent horizon searching problem to our best knowledge. Multigrid methods, owing to their efficiency in solving elliptic partial differential equations (PDEs), have been applied to different problems in numerical relativity, including the initial data construction of binary stellar objects (see, e.g., <cit.>) and more recently to a relativistic hydrodynamic code under the conformally flat assumption <cit.>. It should also be pointed out that a review article introducing the multigrid method for numerical relativists was already written about 40 years ago <cit.>. With the expectation that multigrid methods can solve elliptic equations efficiently, we aim to implement an apparent horizon finder in a multigrid approach and examine whether this could be a better algorithm in terms of speed and robustness, especially for cases when one needs to resolve the apparent horizon with high enough resolution to increase the precision of the inferred properties of the black hole in simulations. 
    

    

    
    This paper is organized as follows. In section <ref>, we first describe the equation to solve for the apparent horizon, and present the implementation of the multigrid algorithm in our apparent horizon finder. Section <ref> includes the test results of our finder in terms of accuracy and efficiency in several benchmark data. Finally, we summarize our work and give future prospects in section <ref>.



§ METHOD AND IMPLEMENTATION


    Locating the apparent horizon requires solving a non-linear elliptic PDE—the expansion equation. In this section, we first introduce the variant ansatz of the expansion equation we employed. We then give a brief overview of the multigrid methods in general and the specific scheme we used in this work.

    

 §.§ The Expansion Equation

        The apparent horizon is a 2-dimensional surface which has a vanishing expansion Θ in its embedding spatial hypersurface. We briefly describe how the expansion function Θ depends on the spacetime geometric variables, and the ansatz we adopted to solve for an apparent horizon.

        

  §.§.§ Notations and definitions

            Let n^a denote the timelike future-pointing unit normal vector to a spacelike hypersurface Σ_t, and s^a denote the spacelike outward-pointing unit normal vector to an arbitrary 2-dimensional smooth closed surface 𝒮 embedded in Σ_t. The spatial 3-metric γ_ab induced on Σ_t and the 2-metric m_ab induced on 𝒮 are given by
            
            
    γ_ab   = g_ab + n_a n_b,
                    
    
                    m_ab   = γ_ab - s_a s_b = g_ab + n_a n_b - s_a s_b ,

            
            where g_ab is the spacetime metric.
            By denoting the unit vector along the future-pointing outgoing null geodesics by k^a = ( s^a + n^a )/√(2), the expansion function can be expressed as (see, e.g., <cit.>),
            
    Θ≡ m^ab∇_a k_b = (D_i s^i + K_ij s^i s^j - K )/√(2),

            where ∇_a and D_a are the covariant derivatives associated with the metric tensors g_ab and γ_ab, respectively; K_ij is the extrinsic curvature, and K denotes its 
            trace (i.e., K ≡γ^ij K_ij).
            𝒮 is called a marginally-trapped surface if Θ=0 everywhere on 𝒮. 
            The apparent horizon is defined to be the outermost of such surfaces.
            
            By assuming that the apparent horizon is topologically equivalent to a 2-sphere and is a star-shaped surface around an interior local coordinate origin[We refer the reader to 
            <cit.> for the removal of the star-shaped assumption in the axisymmetric case.], we use standard spherical coordinates (r, θ, ϕ) to parameterize it. The apparent horizon surface is represented by a horizon function h(θ, ϕ) which measures the radial coordinate distance of the surface from the local origin. The aforementioned outward-pointing normal s^i can be constructed using the gradient of a level-set function F= r- h(θ, ϕ), namely
            
    s_i = λ m_i ≡λ∂_i F = λ (1, -∂_θ h, -∂_ϕ h),

            where λ is the normalization factor such that s_i s^i = 1. Using such 
            parameterization, (<ref>) can be cast into an elliptic PDE in terms of h.

        

  §.§.§ The ansatz adopted

            Following <cit.>, we separate out a linear elliptic operator from the non-linear expansion equation Θ=0, such that it is in a suitable form to be solved by standard relaxation method. We first consider
            
    Δ_θϕ h - 2 h = Δ_θϕ h - 2 h + ρΘ.

            where ρ = ρ(h, θ, ϕ) is a scalar function to be determined and Δ_θϕ is the flat-space Laplacian on a 2-sphere which is given by
            
    Δ_θϕ h = ∂^2_θ h + θ ∂_θ h + 1/sin^2θ∂^2_ϕ h.

            We briefly discuss how the Laplacian shows up in the expansion equation and how the scalar function ρ should be chosen.
            
            A 3-metric tensor γ_ij (r, θ, ϕ) on a spacelike hypersurface Σ_t can be expressed by its conformally related metric γ̅_ij and be further decomposed into the flat metric η_ij along with a tensor field h_ij, i.e.,
            
    γ^ij = ψ^-4γ̅^ij = ψ^-4 (η^ij + h^ij),

            where the conformal factor ψ is given by 
            
    ψ = ( γ_ij/η_ij)^1/12.

            By introducing the covariant derivative D_i associated with the flat metric tensor η_ij and the corresponding Levi-Civita connections Γ^k_ij, the divergent term in the expansion function (<ref>) is
            
            
            
    D_i s^i    = m^ij D_i (λ m_j) = λ m^ij D_i m_j
                    
       = λψ^-4η^ijD_i m_j + λ( ψ^-4 h^ij - s^i s^j ) ∂_i m_j
                    - λ (γ^ij - s^i s^j) Γ^k_ij + λψ^-4η^ijΓ^k_ij m_k.

            
            One can see that the Laplacian, i.e., the linear elliptic part we are looking for, appears in the first term. In spherical coordinates, it reads
            
    η^ijD_i m_j
                    = - 1/h^2( Δ_θϕ h - 2h ).

            The expansion function now becomes
            
            
            
    √(2)Θ   = D_i s^i - (γ^ij - s^i s^j) K_ij
    √(2)Θ   = - λψ^-4/h^2( Δ_θϕ h - 2h ) + λ( ψ^-4 h^ij - s^i s^j ) ∂_i m_j
                    
       = - (γ^ij - s^i s^j) ( λΓ^k_ij + K_ij) + λψ^-4η^ijΓ^k_ij m_k.

            
            Comparing with (<ref>), the combination Δ_θϕh - 2h on the right hand side cancels out if we choose ρ = √(2) h^2 ψ^4 / λ, yielding a PDE in suitable form to be solved by relaxation schemes. We make an additional remark that from 
            (<ref>), the linear elliptic operator in the expansion equation is indeed
            
    h^2 η^ij( D_i m_j + Γ^k_ij m_k ) = - ( ∂^2_θ h + 1/sin^2θ∂^2_ϕ h ),

            but it turns out that such form is not suitable for relaxation methods to work properly.
            
            Finally, the ansatz we adopt in this work in explicit form is given by
            
    Δ_θϕh - (2 - η) h = S ( h, ∂_i h, ∂^2 h; γ_ij, ∂_k γ_ij, K_ij; η),

            with the source term S being
            
            
            
    S    = Δ_θϕ h - (2 - η) h + h^2 ψ^4/λ (√(2)Θ )
                    
       = h^2 ψ^4 [ ( ψ^-4 h^ij - s^i s^j ) ∂_i m_j 
                    - ( γ^ij - s^i s^j ) ( Γ^k_ij m_k + K_ij/λ) ]
                    
       = + h^2 η^ijΓ^k_ij m_k + η h,

            
            where the extra parameter η can speed up convergence in the relaxation method if appropriately chosen <cit.>.


    

 §.§ Implementation of Multigrid Apparent Horizon Finder

    
        

  §.§.§ Multigrid method

            Iterative schemes have been a standard tool for solving general elliptic PDEs numerically. For iterative schemes that use a finite difference discretization setup, there is sometimes a trade-off between choosing a grid of higher resolution for more accurate results and a grid of relatively lower resolution for faster convergence. The persistent low-frequency components in the numerical error is one of the factors that drag the convergent rate, but they can be eliminated rather efficiently at lower resolutions.
            
            A multigrid scheme combines multiple levels of grid with different resolutions to tackle low-frequency components in the error in the coarser grids while retaining good accuracy of the solution in the finest grid. The major components of a multigrid scheme are as follows. 
            Smoothers (e.g., standard relaxations) are applied on each but the coarsest level to eliminate error modes at different spectral range, and a solver is employed at the coarsest level at which an exact/approximate solution can be found more efficiently. With multiple levels of grid, the inter-grid data transfer operators—restrictions and prolongations—are there to bridge successive grid levels by transferring the current solutions and/or error corrections between them. Finally, a cycling scheme is chosen to specify the exact scheduling, when to jump between different levels of grid, of a multigrid solver. We refer the reader to <cit.> for detailed discussion on the multigrid methods.
            
        

  §.§.§ Linear multigrid algorithm

            Owing to the fact that the principle Laplacian in (<ref>) is linear, we choose to use a linear multigrid algorithm, which will be briefly outlined in the following.
            
            Suppose we solve a linear elliptic equation ℒ(u)=f on a uniform grid of 
            spacing k, and we write
            
    ℒ_k (u_k) = f_k,

            where ℒ is the elliptic operator, f is the source term and u is the exact solution. Let u_k denote the intermediate approximate solution, and e_k = u_k - u_k denote the corresponding error. Since ℒ is linear, (<ref>) becomes
            
    ℒ_k (e_k) = f_k - ℒ_k (u_k) ≡ r_k,

            where r_k is called the residual. This residual equation is easier to solve on a coarser grid of spacing 2k if we consider an approximation version of it, i.e.,
            
    ℒ_2k (e_2k) = r_2k,

            where the coarser grid residual r_2k is found by restricting the finer grid residual r_k using the restriction operator ℛ, i.e., r_2k = ℛ (r_k). The solution to (<ref>) can be thought as the correction to the approximate solution u_k we have for the original problem (<ref>). Now denote the approximate solution to (<ref>) by e_2k, we can interpolate it to the finer grid by the prolongation operator 𝒫 such that the approximate solution u_k is updated by
            
    u_k^ new = u_k + 𝒫 (e_2k).

            
            The above procedure is an illustration on a 2-grid structure and can be easily generalized for
            an n-grid solver. In particular, coarser grids can be recursively constructed when solving for the correction in (<ref>).

        

  §.§.§ Grid discretization

            We discretize the apparent horizon surface h(θ, ϕ) uniformly using a vertex-centered grid with N_θ points along the polar direction and 
            N_ϕ points along the azimuthal direction, respectively. The grid points sit at 
            
            
    θ_i = i Δθ = i ( π/N_θ - 1),         i=0,…, N_θ -1;
                    
    ϕ_j = j Δϕ = j ( 2π/N_ϕ - 1),         j=0,…, N_ϕ-1.

            

            We further discretize the apparent horizon equation (<ref>) using a 4th order central finite difference representation for the Laplacian such that it is in a suitable form for conventional relaxation methods.
        
            Special treatments are carried out to the boundary points of this vertex-centered spherical grid. Both the Laplacian Δ_θϕ and the expansion function Θ exhibit singularity at the poles. The polar singularities are avoided by not updating the horizon function h with (<ref>) there. Instead, the values of h at the poles are interpolated by a cubic polynomial using neighboring points. Periodic boundary conditions are applied in the azimuthal direction, i.e., h(θ_i, 2π) = h(θ_i, 0).

        

  §.§.§ Smoothers and inter-grid transfer operators

            In the multigrid algorithm, we need to specify a way to obtain approximate solutions to (<ref>) and (<ref>). We use the Gauss–Seidel line relaxation method <cit.> in the ϕ-direction as a smoother (see section <ref> for more details) at all grid levels including the coarsest level. The smoothing operator applied before restrictions and after prolongations consists of 1 and 2 sweeps of relaxation, respectively. At the coarsest level, the smoothing operator consists of 100 sweeps of relaxation.
        
            For the inter-grid transfer operators, we choose the bi-linear prolongation operator 𝒫 and the full-weight restriction operator ℛ. These operators can be represented in the stencil notations (see, e.g., <cit.>) by
            
    𝒫 = 1/4[ 
                    1     2     1 
                    2     4     2 
                    1     2     1
                    ],     ℛ = 1/16[ 
                    1     2     1 
                    2     4     2 
                    1     2     1
                    ].


        

  §.§.§ Cycling algorithm and the solving procedure

            Among different multigrid cycling algorithms, we choose to use the V-cycle (see Figure <ref>) owing to its simplicity. Successive V-cycles are carried out until some prescribed tolerance is satisfied to locate the apparent horizon. 
           
            
            
            Assuming that the geometric objects γ_ij and K_ij are given on a hypersurface, an initial trial surface h^(0) is specified for the calculation of the source term S according to (<ref>). Note that we do not need the value of S at the poles since these points are excluded from the relaxation domain. One V-cycle is then performed to obtain a new guess surface h^(1) while keeping the source term fixed. The new surface h^(1) is used to update the source term S before entering the next V-cycle. This process is repeated until the maximum change in h between successive trial surfaces, denoted by δ h_∞, 
            and/or the maximum value of the expansion on the final guess surface, denoted by Θ_∞, are less than some tolerance ϵ_h and ϵ_Θ, respectively, within the relaxation domain (i.e., excluding the poles).

        

  §.§.§ Handling of numerical spacetimes

            The ability to process numerical spacetime is important for an apparent horizon finder to be applicable in numerical simulations, where the geometric objects γ_ij and K_ij are given only at discrete locations. To update the source term (<ref>) between each V-cycle with numerical spacetime data, we need the values (in spherical basis) of γ_ij, K_ij and also ∂_k γ_ij on the current trial surface h(θ,ϕ). The current implementation of the interpolator in our finder takes input of γ_ij and K_ij on a uniform Cartesian grid, interpolates γ_ij, ∂_k γ_ij and K_ij in Cartesian basis onto the trial surface using tricubic Hermite spline, and finally transforms the components from Cartesian basis to spherical polar basis. In approximating the first to third order derivatives of γ_ij and K_ij in Cartesian basis which are required for tricubic Hermite interpolation, a standard 4th order finite differencing scheme is used to maintain an overall third order accuracy in the interpolator.  

    

 §.§ Line Relaxation Smoothing Operator

    
        As mentioned in the previous section, we have chosen the Gauss-Seidel line relaxation in the azimuthal direction as the smoothing operator in our multigrid code. Due to the inherent asymmetry 
        between the polar and azimuthal directions of the spherical Laplacian Δ_θϕ, the line relaxation has better convergence, especially near the poles, than a simple point-wise relaxation with or without red-black ordering updates <cit.>. Because of the periodic boundary condition in the azimuthal direction, one step of the ϕ-line relaxation requires solving a cyclic banded diagonal system <cit.>. In particular, a cyclic penta-diagonal system is needed to be solved for each relaxation step if we use a standard 4th order 5-point finite difference formula. This obviously requires much more work than solving a tri-diagonal system, which is the case with a 2nd order 3-point finite difference scheme. In order to take advantage of the efficiency in solving tri-diagonal system while achieving 4th order accuracy, we aim for a compact finite difference scheme for (<ref>) (see e.g., <cit.> and <cit.> respectively for such treatment in Cartesian Poisson equation and polar Helmholtz equation).

        The compact finite difference scheme, also known as the Mehrstellenverfahren discretization (e.g., <cit.>), aims at increasing the accuracy of a finite difference scheme without involving more grid points that are too far away from the current pivoting point. As an example, for the Poisson equation (∂_x^2 + ∂_y^2)u = f on a square Cartesian grid of step size k, the traditional 4th-order finite difference scheme involves the second next neighboring grids, as represented by the stencil notation
        
    1/12k^2[        -1              16       
                -1   16   -60   16   -1        16              -1       ] u_ij
                =
                f_ij,

        while its 4th order compact finite difference representation is <cit.>
        
    1/6k^2[ 
                1     4     1 
                4     -20     4 
                1     4     1
                ] u_ij
                =
                1/12[ 
                0     1     0 
                1     8     1 
                0     1     0
                ] f_ij.

        Notice the smaller stencil size (width) in the compact finite difference scheme. This compact representation can be viewed as the incorporation of a correction to the 2nd order finite difference formula to eliminate the leading error in its Taylor's expansion (see, e.g., <cit.> for a variety of PDEs that such method has been applied to). The compact form for the spherical Poisson equation, to which the apparent horizon equation (<ref>) belongs, can be derived in a similar manner. 
        We refer the reader to <cit.> for a summary of compact scheme for Poisson equation on Cartesian grids.
        
        Since our ultimate goal is to create a stencil that is suitable for reducing the ϕ-line relaxation to solving a tri-diagonal system, we only need to “squeeze” the usual 5-point-wide stencil in the ϕ direction, which we refer to as semi-compact. 
        We derive such semi-compact stencil for a general 2-dimensional spherical Poisson equation in the following[For completeness, we also derive the fully-compact representation in <ref> for a general 2-dimensional Poisson equation on the sphere. We do not use the fully-compact form because it requires the value of the source term f at the poles, but the apparent horizon equation (<ref>) contains the expansion function Θ which is undefined at the poles in spherical coordinates (see (<ref>)).]. 
        Suppose we solve
        
    ∂^2_θ u + θ ∂_θ u + 1/sin^2θ∂^2_ϕ u = f(θ, ϕ).

        First observe that using a 2nd order finite difference method for the third term, we have
        
    1/sin^2θ[2]uϕ 
                = ( 1/sin^2θ[2]uϕ)_ FD2 - k^2/12sin^2θ[4]uϕ + 𝒪(k^4),

        where the subscript FDn indicates that the corresponding term is approximated by an n-th order finite difference scheme, and k is the step size in the ϕ direction. By differentiating (<ref>) twice w.r.t ϕ and substituting the result to (<ref>), we obtain
        
    1/sin^2θ[2]uϕ 
                = ( 1/sin^2θ[2]uϕ)_ FD2
                -k^2/12[ [2]fϕ - ∂^4u/∂θ^2∂ϕ^2 - θ∂^3u/∂θ∂ϕ^2] + 𝒪(k^4).

        With a factor of k^2 in front, the terms inside the square bracket only need to be evaluated at 2nd order accuracy to maintain an overall 4th order accuracy. Finally, assuming that the step size in θ direction is also k for simplicity, the semi-compact 4th order finite difference representation of (<ref>) can be given by
        
    ( [2]uθ + θuθ)_ FD4 + ( 1/sin^2θ[2]uϕ)_ FD2   
    
                -k^2/12[ [2]fϕ - ∂^4u/∂θ^2∂ϕ^2 - θ∂^3u/∂θ∂ϕ^2]_ FD2 + 𝒪(k^4)
                   = f.

        The adaption of this form to our apparent horizon equation (<ref>) is straight-forward, despite the fact that the source term now also depends on the solution. Notice that the stencil is only 3-point-wide in the ϕ direction such that a ϕ-line relaxation can be done by a more efficient tri-diagonal algorithm. 
        
        
        
        We make some final remarks to conclude this section. First, a simple point-wise Gauss-Seidel relaxation using the usual 5-point-wide finite difference scheme can certainly serve as the smoothing operator, but the efficiency of the solver will be much degraded (See Table <ref> and the discussion in Section <ref>). Second, the asymmetry in the spherical Laplacian may change direction from the pole to the equator depending on the grid sizes in the θ- and ϕ-directions <cit.>. In the case when the number of points in the θ-direction is much larger than that in the ϕ-direction, the ϕ-line relaxation may not provide good convergence, especially near the equator, or it may even fail. A remedy is to use a combined relaxation <cit.>, sometimes referred to as the segment relaxation scheme, which applies the θ- and ϕ-line relaxations to different domains according to the direction of asymmetry in the Laplacian. Although this could be relevant for further increasing the efficiency of our finder when applied to spacetimes with axisymmetry, we will continue to use the ϕ-line relaxation here.




§ RESULTS


    We have tested our multigrid apparent horizon finder with several benchmark spacetimes and the results are reported in this section. All tests are run on a 3.6GHz processor and the timings are reported in user CPU time.
    
    The geometric objects γ_ij and K_ij are given analytically in spherical coordinates for the calculation of the source term S (<ref>), except for numerical data where they are interpolated onto the surface. Unless otherwise specified, the initial guess surface is a sphere of coordinate radius 1.5M, where M is the mass of the black hole. The tolerance to stop the iteration process is set to ϵ_h = ϵ_Θ = 10^-8, and the convergence parameter is set to η=0.8. After locating the apparent horizon, its proper area is calculated numerically by <cit.>
    
    A =    ∫_0^2πϕ∫_0^πθ [ ( γ_rr h_,θ^2 + 2γ_rθ h_,θ + γ_θθ) ( γ_rr h_,ϕ^2 + 2γ_rϕ h_,ϕ + γ_ϕϕ) 
            
       - ( γ_rr h_,θ h_,ϕ + γ_rϕ h_,θ + γ_rθ h_,ϕ + γ_θϕ)^2 ]^1/2.


    

 §.§ Kerr-Schild spacetime

    
        For describing a Kerr black hole of mass M=1 and spin parameter a (0≤ a<1), one can use the Kerr-Schild coordinates <cit.>, of which the spatial metric γ_ij and the extrinsic curvature K_ij are given in Cartesian coordinates by
        
        
    γ_ij   = η_ij + 2H ℓ_i ℓ_j;
                
    
                K_ij   = 2α H ℓ^k ∂_k (H ℓ_i ℓ_j) + α [ ∂_i (H ℓ_j) + ∂_j (Hℓ_i) ],

        
        where η_ab =  diag(-1,1,1,1) is the flat metric and α≡√(1+2H) is the lapse function. 
        The function H and the auxiliary components ℓ^i are defined by
            
    H ≡MR^3/R^4+a^2 z^2;
                        ℓ_i = ℓ^i ≡( Rx+ay/R^2+a^2, Ry-ax/R^2+a^2, z/R),
 
            with the parameter R being the Boyer-Lindquist radial coordinate which satisfies
            
    x^2+y^2/R^2+a^2 + z^2/R^2 = 1.

        The analytic forms of γ_ij, ∂_k γ_ij and K_ij, after a coordinate transformation, serve as the input for the calculation of the source term S in (<ref>). There are two horizons at R=r_±≡ M(1±√(1-a^2)) and the apparent horizon is given by the outer one, which is, in spherical coordinates (r, θ, ϕ),
        
    r^2 = r_+^2 ( r_+^2 + a^2)/r_+^2 + a^2 cos^2θ,

        with the corresponding proper area given by 
        
    A_ Kerr = 4π (r_+^2 + a^2).


        Before considering the performance of our apparent horizon finder, we first use the analytic 
        Kerr-Schild spacetime as a test case to show the advantage of using a line relaxation smoother over a point-wise one in a multigrid code as noted in Section <ref>. 
        As shown in Table <ref>, the number of V-cycles required to reach the solution is less when a line relaxation is used. The difference becomes more significant with higher resolution and more grid levels. Combined with the enhancement given by tuning the convergence parameter η, the use of line relaxation over point relaxation leads to an order-of-magnitude  
        overall boost in terms of speed in certain test cases.
        
        We use several resolutions to locate the apparent horizon in Kerr-Schild spacetime with spin parameter a=0.6 and 0.9. The results are reported in Table <ref>. With the same number of grid levels, increasing the resolution slightly increases the number of iteration needed to reach convergence. This, however, can be compensated by introducing more grid levels. Overall, we find that the number of iterations (V-cycles) to reach convergence is more or less the same for a fixed a, independent of the resolutions. The run time reported here is only determined by the searching algorithm, excluding the interpolation routines required in numerical spacetimes (see Section <ref>). We find that the run time is proportional to the total number of grid points on the finest grid level.

        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        

        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        

        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        

    

 §.§ Brill-Lindquist spacetime

        
        The second test for our finder is to locate the apparent horizon in Brill-Lindquist multiple black hole data <cit.>. Under time symmetry and the assumption that the spatial slice is conformally flat, the 3-metric and the extrinsic curvature for N-black-hole spacetime in isotropic coordinates (, θ, ϕ) are, respectively,
        
        
    γ_ij   = ( 1+ ∑_α=1^N M_α/2 |*-*_α|)^4  diag(1, ^2, ^2 sin^2θ ), 
                
    
                K_ij   =0,

        
        where M_α and *_α are the mass at infinite separation and the coordinate position of the α-th black hole, respectively. We consider the case of N=2 and N=3. We use very high resolution of =0mu N_θ× N_ϕ=385×769 (angular spacing <0.5) with 6 levels of grid in these tests without assuming any symmetry of the apparent horizon. We relax the convergence criteria in this section to ϵ_h=e-8 and ϵ_Θ=e-6.


        

  §.§.§ Brill-Lindquist 2-black-hole spacetime

            Let us first consider the Brill-Lindquist data for two black holes. 
            When the two black holes are far away from each other, each black hole possesses its own apparent horizon. With decreasing separation distance d, there exists a critical separation d=d_ crit at which the two black holes start to have a common apparent horizon. One of the benchmark tests is to determine the value of d_ crit. The center of the two black holes are on the z-axis at (0,0,± d/2) in the following.
            
            
            
            
            We first consider two equal-mass black holes with M_1=M_2=1. Our finder reports the critical separation to be d_ crit=1.5323948, which agrees well with the generally accepted values in the literature (e.g., 1.5323949 <cit.>, 1.532 <cit.>). The critical area is found to be A_ crit=196.40795 (agrees with <cit.> to 5 decimal places). At the critical separation, the maximum expansion on the apparent horizon is =0mu Θ_∞ = 2.3e-8. We also report the results for finding the apparent horizon near the critical separation at d=1.532≈ d_ crit for lower resolutions in Table <ref>.
            

            
            
            
            We then consider unequal-mass black holes with M_1=0.2 and M_2=0.8. The initial guess surface for this case is changed to a sphere of coordinate radius =1.5. 
            
            Our finder reports the critical separation to be d_ crit=0.6987161 with the proper area being A=49.688602. The maximum expansion on the surface is =0mu Θ_∞ = 2.1e-7. The value of the critical separation agrees with that reported recently in <cit.> by using an axisymmetric apparent horizon finder. 


        

  §.§.§ Brill-Lindquist 3-black-hole spacetime

            
            
            
            We now consider a system of three black holes which are placed to form an equilateral triangle on the x̅-z̅ plane, each at a coordinate distance R from the origin, and one of them sits on the +z̅-axis. Each black holes has unity mass. We find that the critical parameter for the three black holes to have a common apparent horizon is R_ crit = 1.1954995, which agrees with the result in <cit.> to seven decimal places. 
            The critical area and the maximum expansion on the surface are A=444.75623 and =0mu Θ_∞ = 2.3e-8, respectively. 
            The left panel of Figure <ref> shows the cross section (red dashed line) of the common horizon on the x̅-z̅ plane for the case R=1.195 ≈ R_ crit. 
            The right panel shows the side view of the full surface of the common horizon. 
            
            

    

 §.§ Numerical spacetimes

    
        We prepare numerical spacetimes to further test our finder using the open source code Einstein Toolkit <cit.>. The simulation variables γ_ij and K_ij are post-processed to produce data on a uniform rectangular Cartesian grid before importing to our finder. In the case of dynamical binary black hole simulation where adaptive mesh refinement is turned on, we use the python package Kuibit <cit.> for the post-processing. We also activate the inherent apparent horizon finder in the code,  <cit.>, for comparison. 

        

  §.§.§ Numerical Kerr-Schild spacetime

            We first test our finder with an off-centered Kerr-Schild data. The numerical spacetime data spans the Cartesian grid space of x,y,z∈ [-2.5,2.5] with a grid spacing of Δ x = Δ y = Δ z=0.05. We use a sphere of coordinate radius 1.5 centered at the coordinate origin as the initial guess surface. We use the ILUCG matrix routines in all cases for . 

            

            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            

            We consider a black hole of unity mass and spin a that is located at (0.2,0.2,0.2) in Cartesian coordinates. We tabulate the results for a=0.6 and a=0.8 in Table <ref>. The time to locate the solution surface mainly depends on two factors—the overall convergence of the algorithm and the computational cost within each iteration. The convergence order of the algorithm is important in the sense that as many as 30 interpolations for the geometric variables are required at every grid point at each (outer) iteration. Since we adopt a pure multigrid algorithm, which has a slower convergence than the quadratic convergence in Newton's method, one can see that our finder needs more iterations (number of V-cycles) for convergence, so more time is spent on the interpolation routines when compared to . Notice that the number of iterations is insensitive to the resolution used for fixed a. Had we used the pointwise Gauss-Seidel relaxation as the smoothing operator, the number of iterations would grow linearly with the number of grid points, and the speed of the code would become very slow due to the large amount of interpolations required.


            At high resolutions, however, we find that our finder generally locates the apparent horizon faster than . This is somewhat expected because, at a rough estimation, the operation cost of 
            performing the matrix LU decomposition of the Jacobian equation in  is 𝒪(N^3) but that of one multigrid cycle is 𝒪(N), where N is the total number of grid points used to describe the horizon surface and N=N_θN_ϕ in our case. The total search time is a competition between the two factors mentioned above. The computational cost of each iteration becomes more important with increasing resolutions. Although the overall convergence rate of our multigrid algorithm is not as fast as Newton's method, there is a critical resolution at which our multigrid algorithm starts to work faster. 
            In particular, our finder can be more than two times faster when the 
            angular separation between two neighboring grids is about 1 as shown in Table <ref>. 
            We have performed the same test for additional values of the spin parameter a at 1 resolution, and the runtimes are plotted in Figure <ref>. Although the runtime of our finder grows with a, it takes lesser time to converge in all cases. Recall that the performance of our finder has extra dependence on the convergence parameter η. At the extreme case of a=0.99, the optimal value of η shifts to around 1.45, and the runtime is reduced by half comparing to the case η=0.8.

            
            


            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            


            We further test on the case where the spin parameter is a=0.99, but with the black hole relocated to (0.2,0.2,0). The initial guesses are spherical surfaces centered at the coordinate origin as before, with radius r_0 ranging from 1.5 to 1.6 with increments of 0.01.  
            We found that  is more volatile to the initial guess surface, and is able to find the apparent horizon 
            only at r_0=1.54, 1.57, 1.58 and 1.6 (it fails to converge otherwise). On the other hand, our finder is able to find the correct solution for all these cases. 


        

  §.§.§ Dynamical binary black hole spacetime

            We simulate the merger of an unequal-mass black hole binary
            using Einstein Toolkit <cit.>. The initial mass ratio of the two non-spinning black holes is set to be q≃ 0.7. We use standard evolution schemes and gauge conditions, with multiple layers of mesh refinements to perform the simulation. 
            The finest grids of grid spacing 
            Δ x= Δ y= Δ z=0.03125 cover the two black holes at all times. We export the numerical data to an uniform grid of the finest grid spacing by Kuibit <cit.> to search for the horizons using our finder. For the two individual horizons, the local coordinate origin is chosen by the position of singularity provided by the PunctureTracker thorn in the simulation code, while that of the common apparent horizon is simply set to be the origin of the simulation grid. 

            

            
            
            We report the results of horizon finding at the simulation time when the common apparent horizon first appears in Table <ref> and the corresponding snapshot of the three horizon surfaces in Figure <ref>. The horizon surfaces are modeled by 5 angular resolution 
            (N_θ×N_ϕ = 37×73). In Figure <ref>, the common apparent
            horizon is represented by the purple surface, inside of which are the two individual horizons represented by the green and blue surfaces.                      
            In terms of the proper areas, our finder agrees with the default  in the simulation code to within 0.2%. The difference is the biggest for the horizon of the smaller black hole, as the spherical grids become overly stacked at that radius.





§ CONCLUSIONS


    We have presented a new apparent horizon searching algorithm using a multigrid method in this work. We have evaluated the performance of our finder on analytic spacetimes, including the Brill-Lindquist data with two or three black holes. We have recovered the results of the critical parameter for a common apparent horizon to appear in these spacetimes. We have also tested our finder on numerical spacetimes, namely the off-centered Kerr-Schild data and the inspiraling binary black hole spacetime. We have demonstrated that, in terms of speed, the multigrid searching algorithm begins to outperform the currently fastest algorithm (Thornburg's  <cit.>) at high resolutions. We have also shown that our finder is capable of capturing the first appearance of the common apparent horizon in merging black hole binary simulation.
    
    From a numerical perspective of viewing our finder as a general multigrid solver for Poisson equations in spherical coordinates, there are some notable aspects. We have verified that the line relaxation scheme is preferable (in terms of convergence) to a pointwise relaxation scheme, even when dealing with a solution-dependent nonlinear source term. 

    We have derived the (semi-)compact 4th order finite difference scheme for solving the spherical Poisson equation. The derivation follows a similar 
    procedure used for Cartesian grids (see <cit.> and references therein). 
    Our study demonstrates the applicability of this scheme in reducing the complexity of a 
    line relaxation technique on spherical coordinates, from a penta-diagonal matrix inversion to a more efficient tri-diagonal matrix inversion. 
    We believe the formula derived should be directly applicable to other 2-dimensional Poisson equations on the sphere.
    
    We make some final remarks on how our multigrid algorithm can be further improved. 
    First, the segment relaxation scheme (See Section <ref>) could be used to improve the convergence. 

    We only use ϕ-line relaxation as the smoother in our multigrid solver, but the line relaxation direction is preferred to be aligned with the anisotropy direction in the differential operator as it could give better convergence <cit.>. 
    Second, using a variant multigrid algorithm might also improve the convergence. We solve the non-linear apparent horizon equation by turning it into a linear elliptic equation with a non-linear source term. In this way, the nonlinearity of the principle equation is avoided.
    However, by not addressing it directly, the rate of convergence could be compromised. 
    The common multigrid methods to tackle nonlinearity directly are the full approximation scheme
    <cit.> or Newton-multigrid <cit.> algorithm. We have not adopted these methods in the current work since they in principle bring up technical implementation issues in the case of the apparent horizon equation. Exploring these alternative multigrid schemes could be one future investigation direction for increasing the efficiency of our multigrid code. 





§ DATA AVAILABILITY STATEMENT

No new data were created or analyzed in this study.



§ ACKNOWLEDGMENTS

We thank Kenneth Chen and Yun-Kau Lau for useful discussions. LML is also grateful for the hospitality of the Morningside Center of Mathematics, Chinese Academy of Sciences, where the idea of this project was initiated. 



§ ORCID IDS


Hon-Ka Hui 0000-0002-2815-527X <https://orcid.org/0000-0002-2815-527X>

Lap-Ming Lin 0000-0002-4638-5044 <https://orcid.org/0000-0002-4638-5044>



§ REFERENCES



This preprint is prepared with the IOP Publishing  preprint class.






§ FULLY-COMPACT 4TH-ORDER FINITE DIFFERENCE SCHEME FOR SPHERICAL POISSON EQUATIONS


    Suppose we solve a general Poisson equation in spherical coordinates,
    
    ( [2]θ + θθ + 1/sin^2θ[2]ϕ) u(θ, ϕ) = f(θ, ϕ),

    on a spherical grid of grid size[We follow the standard notation to use h for the grid spacing, which should not be confused with the horizon function h(θ, ϕ) used elsewhere in the paper.] h=Δθ and k=Δϕ. Using a 2nd-order finite difference approximation, the first two terms on the left side are
    
    [2]uθ + θuθ
            = ( [2]uθ + θuθ)_ FD2
            - h^2/12( [4]uθ + 2θ[3]uθ) + 𝒪(h^4).

    We then seek an expression to replace the third and fourth order derivatives in the second parenthesis. Differentiating (<ref>) w.r.t θ gives
    
    [3]uθ + θ[2]uθ - 1/sin^2θuθ
            -2cosθ/sin^3θ[2]uϕ + 1/sin^2θ∂^3u/∂θ∂ϕ^2
            = fθ.

    The second derivative reads
    
    [4]uθ + θ[3]uθ - 2/sin^2θ[2]uθ + 2cosθ/sin^3θuθ 
            + 6-4
            sin^2θ/sin^4θ[2]uϕ
    
            -2cosθ/sin^3θ∂^3u/∂θ∂ϕ^2
            -2cosθ/sin^3θ∂^3u/∂θ∂ϕ^2 + 1/sin^2θ∂^4u/∂θ^2∂ϕ^2
            = [2]fθ.

    Combining the last two equations, we have
    
    [4]uθ + 2θ[3]uθ 
            =   ( θfθ + [2]fθ)
            + ( 1+1/sin^2θ) [2]uθ 
            - cosθ/sin^3θuθ
       + 2sin^2θ-4/sin^4θ[2]uϕ 
            + 3cosθ/sin^3θ∂^3u/∂θ∂ϕ^2 
            - 1/sin^2θ∂^4u/∂θ^2∂ϕ^2.

    The right hand side of this equation can be evaluated by 2nd order finite difference approximation and be substituted back to (<ref>) while the overall error is maintained at 𝒪(h^4). We can then combine the result with (<ref>) to obtain the fully-compact 4th order finite difference approximation for the spherical Poisson equation (<ref>). The explicit expression is
    
    ( [2]uθ + θuθ)_ FD2 + ( 1/sin^2θ[2]uϕ)_ FD2
            -k^2/12[ [2]fϕ - ∂^4u/∂θ^2∂ϕ^2 - θ∂^3u/∂θ∂ϕ^2]_ FD2
    
            -h^2/12[
                θfθ + [2]fθ 
                + ( 1+1/sin^2θ) [2]uθ
                - cosθ/sin^3θuθ
                + 2sin^2θ-4/sin^4θ[2]uϕ
    
                + 3cosθ/sin^3θ∂^3u/∂θ∂ϕ^2 
                - 1/sin^2θ∂^4u/∂θ^2∂ϕ^2]_ FD2
            + 𝒪(h^4, h^2k^2, k^4)    = f.




