{'fraction_non_alphanumeric': 0.031634224544926655, 'fraction_numerical': 0.009795753490368351, 'mean_word_length': 4.403082799072432, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Initial Conditions for Numerical Relativity -- Introduction to numerical methods for solving elliptic PDEs

Numerical relativity became a powerful tool to investigate the dynamics of binary problems with black holes or neutron stars as well as the very structure of General Relativity. Although public numerical relativity codes are available to evolve such systems, a proper understanding of the methods involved is quite important. Here we focus on the numerical solution of elliptic partial differential equations. Such equations arise when preparing initial data for numerical relativity, but also for monitoring the evolution of black holes. Because such elliptic equations play an important role in many branches of physics, we give an overview of the topic, and show how to numerically solve them with simple examples and sample codes written in C++ and Fortran90 for beginners in numerical relativity or other fields requiring numerical expertise.


Introduction
Numerical relativity is now a mature science, the purpose of which is to investigate non-linear dynamical spacetimes. A traditional example of the application and importance of numerical relativity concerns the modelling of gravitational-wave emission and consequent detection. In order to detect gravitational waves(GWs) from black hole-neutron star (BH-NS), BH-BH or NS-NS binaries, one needs to accurately understand the waveforms from these sources in advance, because their signals are quite faint for our detectors. 1 To understand why the problem is so difficult, consider Newtonian gravity, as applied to systems like our very own Earth-Moon. In Newton's theory, binary systems can move on stable, circular or quasi-circular orbits. However, in binary systems heavy enough or moving sufficiently fast, the effects of General Relativity become important, and the notion of stable orbits is no longer valid: GWs take energy and angular momentum away from the system and energy conservation implies that the binary orbit shrinks until finally the objects merge and presumably form a final single object. Accordingly, the evolution of binary stars can be divided into an inspiral, merger and a ring-down phase. 2 In the inspiral phase, GW emission is sufficiently under control by using slowmotion, Post-Newtonian expansions because the stars are distant from each other and their gravitational forces can be described in a perturbation scheme. 3,4 The ringdown phase describes the vibrations of the final object. Because of the uniqueness theorems, 5 GWs can be computed by BH perturbation methods. 6,7 Advanced BH perturbation methods are reviewed in Ref. 8 Numerical Relativity enables us to obtain the GW form in all phases. [9][10][11] Furthermore, we note that techniques of numerical relativity are also available in a variety of contexts. For example, but by no means the only one, it became popular to investigate the nature of higher dimensional spacetimes, 12 most specially in the framework of large extra dimensions. [13][14][15][16] It was pointed out that a micro BH can be produced from high energy particle collisions at the Large Hadron Collider(LHC) and beyond, 17,18 and while some works use shock wave collisions [19][20][21] the full-blown numerical solution is clearly desirable to investigate the nature of gravity with high energy collisions in four dimensional spacetime. [22][23][24][25] If we consider spacetime dimensions higher than four in our simulations, larger computational resources will be required. However, it can be reduced to four dimensional problem with small changes assuming the symmetry of spacetimes. [26][27][28][29][30][31] As a result, to investigate the nature of higher dimensional spacetimes is in the scope of numerical relativity. [32][33][34] Fortunately, open source codes to evolve dynamical systems with numerical relativity are available. [35][36][37][38] All that one needs to do is to prepare the initial data describing the physics of the problem one is interested in. Here, we explain precisely how this is achieved, to prepare initial data for numerical relativity in this paper. Briefly, it amounts to solving an elliptic partial differential equation (PDE) and we explain how to solve the elliptic PDE from the numerical point of view of beginners in numerical studies.

ADM formalism
In numerical relativity, we regard our spacetimes as the evolution of spaces. We begin by showing how to decompose the spacetime into timelike and spacelike components in the ADM formalism. [39][40][41] Then, we derive evolution equations from Einstein's equations along the lines of York's review. 42

Decomposition
First, we introduce a family of three-dimensional spacelike hypersurfaces Σ in fourdimensional manifold V . Hypersurfaces Σ are expressed as the level surfaces of a scalar function f and are not supposed to intersect one another. We can define a one-form Ω µ = ∇ µ f which is normal to the hypersurface.
Let g µν be a metric tensor in four-dimensional manifold V . The norm of one-form Ω µ can be written by a positive function α as where α is called lapse function. We define a normalized one-form by ω µ = αΩ µ , g µν ω µ ω ν = −1.
The orthogonal vector to a hypersurface Σ is written by whose minus sign is defined to direct at the future and we note that this timelike vector satisfies n µ n µ = −1 by definition.

Induced metric
The induced metric γ µν on Σ and the projection tensor ⊥ µ ν from V to Σ are given by the four-dimensional metric g µν , where one can show n µ γ µν = 0, which yields that timelike components of γ µν vanish and only spacelike components γ ij exist. The induced covariant derivative D i on Σ is also defined in terms of the four-dimensional covariant derivative ∇ µ .
where ψ and W λ denote arbitrary scalar and vector on Σ. By a straightforward calculation, one can show that the induced covariant derivative satisfies D i γ jk = 0.

Curvature
Riemann tensor on Σ is defined using an arbitrary vector W i by where [ ] denotes the antisymmetric operator for indices and R ijkℓ denotes Riemann tensor on Σ. Ricci tensor is determined by the contraction of the induced metric and Riemann tensor on Σ. Ricci scalar is also determined by the contraction of the induced metric and Ricci tensor. We define the extrinsic curvature on Σ, which describes how the hypersurface is embedded in the manifold V . The extrinsic curvature is defined by where ( ) denotes the symmetric operator for indeces. One can also show that the extrinsic curvature is spacelike by multiplying the normal vector n µ in the same manner as γ µν . In addition, by the definition of the projection tensor, we obtain the following relation between the covariant derivative of the normal vector and their projection, where the relation n λ n λ = −1 is used in the last equation. Then, the extrinsic curvature can be rewritten by where £ n γ µν denotes the Lie derivative of the tensor γ µν along the vector n µ . The geometrical nature of the three-dimensional hypersurfaces can be determined by the induced metric and extrinsic curvature on Σ. K ij and γ ij must satisfy the following geometrical relations to embed Σ in V .

Geometrical relations
We derive geometrical relations by the projection of the four-dimensional Riemann tensor to the hypersurface Σ. First, in order to obtain the relation between the four-dimensional Riemann tensor R µνρλ and the three-dimensional Riemann tensor R ijkℓ , we rewrite the definition of R ijkℓ with R µνρλ and the extrinsic curvature, Eq. (13) is called Gauss' equation. Secondly, we project the four-dimensional Riemann tensor contracted by an orthogonal normal vector n λ .
Eq. (14) is called Codazzi's equation. Finally, we consider a Lie derivative of the extrinsic curvature to the time direction. We define a timelike vector t µ with a lapse function and a shift vector β µ which satisfies Ω µ β µ as Then, with the Lie derivative along t µ and β µ , we rewrite the four dimensional Riemann tensor contracted by two orthogonal normal vectors as which Eq. (16) is called Ricci's equation.

Decomposition of Einstein's equations
Let us now use the geometric relations to decompose Einstein's equations. Let us for convenience define the Einstein tensor G µν , where G denotes the gravitational constant and c denotes the speed of light and hereafter we set G = c = 1 for simplicity. We start by decomposing the energy momentum tensor as where ρ ≡ T µν n µ n ν , j µ ≡ − ⊥ ρ µ T ρλ n λ and S µν ≡⊥ ρ µ ⊥ λ ν T ρλ . We multiply Gauss' equation (13) by an induced metric γ ik and obtain In addition, Eq. (19) contracted by γ jℓ gives twice as much as the Einstein's tensor contracted by two orthogonal normal vectors n µ and n ν . Then, we obtain Similarly, Codazzi's equation (14) contracted by γ jk results in  (17) can also be expressed with the trace of the energy momentum tensor T ≡ g µν T µν as The projection of Einstein's equations (22) yields where S = γ ij S ij . Ricci's equation (16) is rewritten with Eq. (19) and Eq. (23) as

Propagation of Constraints
In ADM formalism, Einstein's equations are regarded as evolution equations in time with the geometrical constraints on each hypersurface. In general, it is numerically expensive to guarantee the constraints on each step because we must solve elliptic PDEs as described in Sec. 3.1. However, in principle, one does not have to solve the constraint equations if the initial data satisfy the constraints. [43][44][45] This works as follows. We first define the following quantities, where C = 0 corresponds to the Hamiltonian contstraint, C µ = 0 correspond to the momentum constraints and C µν = 0 denote the evolution equations in ADM formalism. Einstein's equations can be decomposed in terms of C, C µ and C µν as Thanks to the Bianchi identity which is a mathematical relation for the Riemann tensor, the covariant derivative of Einstein's tensor vanishes. Besides, the covariant derivative of the energy momentum tensor which appears in the right-hand-side of Einstein's equations denotes the energy conservation law, Let us project the covariant derivative of Einstein's equations to n µ direction and to the spatial direction with ⊥ ρ µ .
where D µ denotes the covariant derivative with respect to the induced metric, noting that C µ and C µν are spatial. Thus, we show the propagation of constraints along the timelike vector as Therefore, we can keep the Hamiltonian and momentum constraints satisfied during evolution, as long as we evolve the initial data satisfying the constraints C = 0 and C µ = 0 on Σ by the evolution equation C µν = 0. It should be noted that the ADM evolution equations are numerically unstable and not suitable for numerical evolutions; instead one uses hyperbolic evolution equations for time integration, for example, BSSN and Z4 evolution equations. 26, 46-49

Initial Condition for Numerical Relativity
As emphasized in Sec. 2.2, initial data cannot be freely specified, as it needs to satisfy the Hamiltonian and momentum constraints on a hypersurface Σ. In general, the problem of constructing initial data is called "Initial Value Problem". The standard method for solving an initial value problem is reviewed in Ref. 50,51 In this section, we derive the equations for the initial value problem and then introduce some examples as initial data for numerical relativity.

Initial Value Problem
There are twelve variables(γ ij , K ij ) as the metric part to be determined and four constraint equations in ADM formalism. One can obtain four variables by solving constraints after assuming eight variables by physical and numerical reasons.

York-Lichnerowicz conformal decomposition
To begin with, we introduce the conformal factor ψ as where we define detγ ij ≡ 1 and we have one degree of freedom in ψ and five degrees of freedom inγ ij . By the conformal transformation, the following relations between variables with respect to γ ij andγ ij are immediately given by whereD i ,Γ i jk ,R and△ are respectively covariant derivative, Ricci scalar, Christoffel symbol and Laplacian operator defined by△ψ =γ ijD iDj ψ with respect toγ ij .

Transverse-Traceless decomposition
As for the extrinsic curvature, we start by decomposing it into a trace and a tracefree part, where γ ij A ij = 0 and K = γ ij K ij and we have one degree of freedom in K and five degrees of freedom in A ij . Then, we also define the conformal transformation for A ij as According to the definition of derivatives with respect to γ ij andγ ij , we obtain the following relation, In addition, we decompose the conformal traceless tensorÃ ij into a divergenceless part and a "derivative of a vector" W j part. Hereafter we assume the divergenceless part vanishes for simplicity. The conformal traceless extrinsic curvature is described byÃ The covariant derivative ofÃ ij is written bỹ where we used the definition of Riemann tensor.

Constraints as initial value problem
With the above conformal transformation, the Hamiltonian and momentum constraints are rewritten as

Schwarzschild Black Hole
Let us consider an exact solution of Einstein's equations as initial data for numerical relativity. The Schwarzschild BH is the simplest BH solution in static and spherically symmetric spacetimes. 52 The line element of the Schwarzschild BH in spherical coordinates (r, θ, φ) is given by where we define f 0 with BH mass M , Let us define the coordinate transformation bȳ where r denotes the isotropic radial coordinate and we introduce a scalar function ψ 0 . Then, we solver under the coordinate transformation and obtain ψ 0 and the relation betweenr and r as After straightforward calculations, the line element is rewritten by where η ij denotes the flat metric and we define α 0 . In the isotropic coordinates, all spatial metric components remain regular, in contrast to the ones in the standard Schwarzschild coordinates. The range [2M <r < ∞] in the spherical coordinates corresponds to [ M 2 < r < ∞] in the isotropic coordinates. In addition, when we change to a new coordinater ≡ (M/2) 2 /r, we obtain the same expression as . The solution is inversion symmetric at r = M/2 and corresponds to the Einstein-Rosen bridge. 53 Obviously, the extrinsic curvature K ij of the Schwarzschild BH vanishes because the spacetime is static, and therefore the momentum constraints are trivially satisfied. Besides, the Hamiltonian constraint is also satisfied as because the Schwarzschild BH is an exact solution of Einstein's equations.

Puncture Initial Data
One can analytically solve the momentum constraints with the following conditions, The derivative operator becomes quite simple assuming conformal flatness. We also note that Eq. (42) and Eq. (43) are decoupled with K = const. condition.

Single Black Hole
Next, let us consider a BH with non-zero momentum(P i = 0), for which the momentum constraints become non-trivial. However, a solution to conditions (52) can still be found. In this case, the momentum constraints are given by We have a simple solution to satisfy Eq. (53) as where P i and S i are constant vectors corresponding to the momentum and spin of BH and n i ≡ x i /r denotes the normal vector. Then, we obtain the Bowen-York extrinsic curvature 54 by substituting Eq. (54) into Eq. (40), On the other hand, to satisfy the Hamiltonian constraint (42) we must in general solve an elliptic PDE, even if simple-looking, Let us define the function u as a correction term relative to the Schwarzschild BH, We can regularize the Hamiltonian constraint (56) when ru is regular at the origin.
A ij is at most proportional to r −3 at and ψ is proportional to r −1 , so that the divergent behavior ofÃ ijÃ ij is compensated by the ψ −7 term at the origin. Therefore, the Hamiltonian constraint yields

Multi Black Holes
We can easily prepare the initial data which contains many BHs without any momenta under the condition (52) because the Hamiltonian constraint is the same as Eq. (51) and we know that the following conformal factor satisfies the Laplace equation.
where M n and x n denote the mass and position of n-th BH, respectively. The initial data defined with ψ = ψ M andÃ ij = 0 is called Brill-Lindquist initial data. 55 As for BHs with non-zero momenta, we also use the Bowen-York extrinsic curvature and the same method for the Hamiltonian constraint as ψ ≡ ψ M + u, In principle, it is possible to construct initial data for multi BHs with any momenta and spins by solving an elliptic PDE. 56

Kerr Black Hole
It should be noted that we have the exact BH solution of a rotating BH for Einstein's equations and we can also use it as initial data. The Kerr BH is an exact solution of Einstein's equations in stationary and axisymmetric spacetime. 57 The line element of the Kerr BH in Boyer-Lindquist coordinates 58 is defined by where where M and a denote the mass and spin of BH respectively. ∆ vanishes when the radial coordinate r BL is at the radius of the inner or outer horizon r ± , which is a coordinate singularity. Let us introduce a quasi-isotropic radial coordinate in the same manner as for the Schwarzschild BH by Thus, the line element of the Kerr BH yields The spatial metric components in the quasi-isotropic coordinates also remain regular. 59,60 One can show that the extrinsic curvature of the Kerr BH in the quasiisotropic coordinates is given by which comes from the shift vector β φ .

Apparent Horizon Finder
Now we can perform long-term dynamical simulations containing BHs with numerical relativity. For the sake of convenience, we usually use the apparent horizon (AH) to define the BH and investigate the nature of BH during the evolution. In this section, we introduce the concept of AH and derive the elliptic PDE to determine the AH.

Apparent Horizon
The region of BH in an asymptotic flat spacetime is defined as the set of spacetime points from which future-pointing null geodesics cannot reach future null infinity. 61 To find the BH, one can use the event horizon(EH) which is defined as the boundary of such region. It is possible to determine the EH by the data of the numerical simulation because in principle, one can integrate the null geodesic equation for any spacetime points forward in time during the evolution, where x µ and λ denote the coordinates and the affine parameter. The numerical cost to find the EH is normally high, because we need global metric data. 62 We define a trapped surface on the hypersurface Σ as a smooth closed twodimensional surface on which the expansion of future-pointing null geodesics is negative. The AH is defined as the boundary of region containing trapped surfaces in the hypersurface and is equivalent to the marginally outer trapped surface on which the expansion of future-pointing null geodesics vanishes. 63 The EH is outside the AH if the AH exists. 61 We often use the AH to find the BH in the numerical simulation instead of the EH because the AH can be locally determined and then the numerical cost is lower compared with finding the EH.
Let us introduce the normal vector s i to the surface and define the induced two-dimensional metric as where γ µν denotes the induced metric on the three-dimensional hypersurface Σ. The null vector is described with s i and the normal vector to Σ by Then, the following equation should be satisfied on the AH by definition.
where Θ denotes the expansion of null vector and D i denotes the covariant derivative with respect to γ ij .

Apparent Horizon Finder
We can find the AH during the dynamical simulation by solving Eq. (73). [64][65][66] Let us define the radius of the AH by Thus, the normal vector s i can be described with h(θ, φ) bỹ wheres i is introduced for convenience and we raise their indeces of s i ands i by γ ij andγ ij respectively. Incidentally, the divergence of the normal vector is given by where γ denotes the determinant of γ ij . Therefore, we obtain the equation to determine the AH as the elliptic PDE consisting of first and second derivatives of h(θ, φ). Note that because the AH equation is originally a non-linear elliptic PDE, we change the AH equation to the flat Laplacian equation with non-linear source term, 64 which has the advantage of fixing the matrix with diagonal dominance mentioned in Sec. 5.2.2. Specifically, we solve the following equation.
where ζ denotes a constant to be chosen by the problem and the source term is given by the flat laplacian term and the AH equation as where ξ ij ≡γ ij − η ij .

Mass and Spin of Black Hole
The area of the AH is defined by where S denotes the surface of AH. We also compute the quantities related to the AH, the polar and equatorial circumfential length(C p , C e ).
If the BH relaxes to a stationary state during the evolution, the BH would be the Kerr BH because of no-hair theorem. The quantities related to the AH of the Kerr BH can be obtained by where M BH , a and r + denote the mass, spin and outer horizon radius defined by r + = 1 + √ 1 − a 2 and E(z) denotes an elliptic integral defined by

Numerical Methods for solving elliptic PDEs
Constructing the initial data for numerical relativity is, in general, equivalent to solving the elliptic PDEs (42) and (43) with appropriate conditions. In order to solve a binary problem with high accuracy, the spectral method should be the standard method for solving an elliptic PDEs. In fact, there are useful open source codes, for example, TwoPuncture 67, 68 and LORENE. 69 Futhermore, we have to solve another elliptic PDE to find the BH in simulations within numerical relativity as described in Sec. 4. In this case, fast methods to solve the ellitptic PDE are preferred. Because there are many elliptic PDE solvers, the method has to be chosen according to the specific purpose. In this section, we introduce some classical numerical methods for beginners. It would also be the basis for Multi-Grid method mentioned in Appendix B.

Discretization
We should discretize our physical space to the computational grid space by finite difference method because we cannot take continuum fields into account on the computer. Consider first one-dimensional problems for simplicity, and introduce the grid interval ∆x. Taylor expansion of a field Q(x) is given by Thus, the derivative of the field Q(x) can be expressed as where Q j+1 , Q j and Q j−1 denote Q(x + ∆x), Q(x) and Q(x − ∆x) respectively and both accuracies of the forward and backward difference method for derivatives are O(∆x). In addition, the central difference method whose accuracy is O(∆x 2 ) can be defined by both Taylor expansions as Similarly, the second-order derivative of Q(x) is written by One can increase accuracy of the calculation by using many points. For example, using five values Q(x + 2∆x), Q(x + ∆x), Q(x), Q(x − ∆x) and Q(x − 2∆x) around x, the fourth-order accuracy scheme are defined by noting that higher accuracy scheme can also be defined. We also note that we can discretize our space in more than two dimensions in the same way.

Relaxation Method
Hereafter, let us focus on Poisson equations (△ψ = S) with a field ψ and a source S as elliptic PDEs. These are a sufficiently general and complex class of problems that they embody all necessary elements to solve Poisson equation for constructing initial data for numerical relativity or finding an apparent horizon of BH. We explain how to solve general elliptic PDEs in Appendix A. One of the simple methods to solve elliptic PDEs, so-called relaxation method, 70-72 is described in this section. Let us introduce a virtual time τ to solve an elliptic PDE and our equation of elliptic type can be transformed to the equation of parabolic type as which denotes the original Poisson equation after ψ relaxes by iteration. We adopt Cartesian coordinates in three-dimensional spaces and discretize the Poisson equation with second-order accuracy as ψ n+1 j,k,l − ψ n j,k,l ∆τ = ψ n j+1,k,l − 2ψ n j,k,l + ψ n j−1,k,l ∆x 2 + ψ n j,k+1,l − 2ψ n j,k,l + ψ n j,k−1,l ∆y 2 + ψ n j,k,l+1 − 2ψ n j,k,l + ψ n j,k,l−1 where the superscript n denotes the label of virtual time and the subscript j, k and l denote labels of x−, y− and z−direction, respectively. Therefore, the field in the next step of the iteration is determined by ψ n+1 j,k,l = 1 − 2∆τ 1 ∆x 2 + 1 ∆y 2 + 1 ∆z 2 ψ n j,k,l − ∆τ S j,k,l + ∆τ ψ n j+1,k,l + ψ n j−1,k,l ∆x 2 + ψ n j,k+1,l + ψ n j,k−1,l ∆y 2 + ψ n j,k,l+1 + ψ n j,k,l−1 We continue to update the field ψ by the above expression until ψ relaxes and obtain the solution of the Poisson equation (△ψ = S).

Jacobi Method
In order to discuss method in practice, we consider one-dimensional Poisson equation and discretize it as We consider Eq. (98) as the equation to determine ψ j , where the superscript n denotes the label of time step and we attach the label J on the field of next time step in Jacobi's method. Thus, we repeat updating the field until ψ converges. In other words, the flowchart of Jacobi method is as follows.
1. Give a trial field ψ n . 2. We obtain a new field ψ n+1 by Eq. (99). 3. Set the obtained field to a new trial field. 4. Repeat these steps (1.-3.) until the change of ψ is within a numerical error.

Matrix expression
Discretized elliptic PDEs can be expressed by matrices and vectors. Once we describe the elliptic PDE via a matrix expression, the problem involves solving the inverse of the matrix. For example, a matrix expression for Jacobi method is given as follows. We introduce a solution vector ψ I and source vector b I . Then, Eq. (98) can be expressed as where A IJ is the coefficient matrix corresponding to the Laplacian operator and the first and last rows of A IJ denote boundary conditions to be determined by physics. It is formally expressed by A IJ ψ J = b I and the problem leads to solving the inverse of coefficient matrix as ψ J = A −1 IJ b I . There are many methods to numerically solve the inverse of matrices. In general, Jacobi method for an arbitrary coefficient matrix is expressed as We note that other elliptic PDEs may not be solvable by the Jacobi method because the iteration is not always stable; this can be shown by the von Neumann numerical stability analysis. However, the Poisson equation is fortunately stable, which is equivalent to that the matrix is diagonally dominant.

Gauss-Seidel Method
In iterative methods to solve Poisson equations as Jacobi method, it depends on a trial field ψ n how fast we obtain solutions. We usually expect that it becomes better solution as iteration step goes forward. In order to obtain a closer trial field to the solution, we should actively use updated values. Thus, by Gauss-Seidel method, we can determine a next trial field as where GS denotes that the field is determined by Gauss-Seidel method. In addition, the matrix expression for Gauss-Seidel method is given by

SOR Method
It turns out that the Poisson equation is faster to solve with Gauss-Seidel method than with the Jacobi method, as shown later. Although the speed with which the numerical solution converges depends on the trial field in iterative methods, the Gauss-Seidel method gives a "better" field than Jacobi's in that respect. Thus, it is possible to accelerate convergence by specifying trial guess of the field more aggressively. This method is called Successive Over-Relaxation(SOR) method, which is defined by where the superscript S denotes the label of SOR method and ω denotes an acceleration parameter whose range is defined in 1 ≤ ω < 2 by the stability analysis. When we set the acceleration parameter as unity, SOR method is identical to the Gauss-Seidel method by definition.

Results
In this section, we introduce sample codes to solve Poisson equations using different methods. These codes are sufficiently general that they can be applied to other problems in physics, provided one slightly changes the source term and boundary conditions.

Code Tests
As tests for our codes, we use the following analytical solutions. We show numerical results of Poisson equations with simple linear source and sufficiently non-linear source. In addition, the code to find the AH of the Kerr BH is also shown as an example. Some sample codes parallelized with OpenMP are also available Appendix C.

Linear source
Let us consider simple source term for the Poisson equation as In numerical computation, we set the range as 0 ≤ x ≤ 1 and boundary conditions by dψ dx x=0 = 0 , Neumann B.C., Then, we obtain the analytical solution by integrating Eq. (106) twice with boundary conditions (107), Arbitrary initial guess for the solution can be given and we set ψ(x) = 1 at initial for those boundary conditions. We set the resolution of the computational grid as ∆x = 1/100. Fig. 1 (a) shows the numerical solution by Jacobi method as compared to the analytical solution. We note that the accuracy of the numerical result depends on the computational resolution and how the accuracy increase with resolution depends on the scheme of discretization; Fig. (b) is compatible with second-order accuracy. We compare Poisson solvers in Fig. 1 (c) by the time steps needed to obtain the solution. Curves show the difference of methods and we choose the Jacobi method, Gauss-Seidel method, SOR methods with ω = 1.5 and ω = 1.9. The SOR method gives the solution about 10-100 times faster than Jacobi method, and depends on the acceleration parameter ω.

Non-linear source
Next, let us consider a weak gravitational field, namely Newtonian gravitational source. A gravitational potential Φ can be determined by the Poisson equation, where we omit the Newton constant by using G = 1 units. Suppose gravitational sources are distributed with spherical symmetry as where ρ 0 is a constant. Corresponding Poisson equation is rewritten by Thus, we obtain the analytical solution of the source (110) by solving the equations separately as the region (r > 1) with the boundary condition Φ → 0 at infinity and (r ≤ 1) with the regularity condition at the origin.
We consider this analytical solution to test our code. The Poisson equation with spherical symmetry can be regarded as one dimentional Poisson equation with the non-linear source in our method, whose range to be considered as 0 ≤ x ≤ 10 and boundary conditions are set as dΦ dr r=0 = 0 , Neumann B.C., where the Robin boundary condition is chosen because we expect Φ → r −1 at large distance. Fig. 2 (a)

Apparent Horizon of Kerr Black Hole
Let us apply our code for Poisson solver to solving the AH of Kerr BH. The AH equation (79) should be reduced to simpler equation with axisymmetry ∂ φ = 0. The normal vector s i is defined with axisymmetry and the normalization C is determined by the Kerr metric as To be concrete, we note that the non-trivial part of the AH equation with axisymmetry in isotropic coordinates can be written by A ,r = 4 r 2 BL + a 2 r BL dr BL dr − ∆ ,r a 2 sin 2 θ, A ,θ = −2∆a 2 sin θ cos θ, On the other hand, we can also solve the AH in Boyer-Lindquist coordinates, only to change the following part.
In Fig. 3 (a), we show the surface of AH of the Schwarzschild BH in isotropic coordinates with the code "sor AHF SBH ISO.f90" and show the three dimensional shape of the AH in 1/8 spaces of computational grid. Fig. 3 (b) shows the difference of the shape on x-z two dimensional plane among different spin parameters with the code "sor AHF KBH ISO.f90". The AH radius shrinks as the spin of BH increses.

Kerr Black Hole and Single Puncture Black Hole
As the last example, let us compare Kerr BH to single puncture BH with a spin as initial data for numerical relativity. A Kerr BH in quasi-isotropic coordinates can be used as the initial data discussed in Sec. 3.4. A single puncture BH is obtained by solving the Hamiltonian constraint (58) without any momenta P i = 0 and with a spin S z in the Bowen-York extrinsic curvature (55).
In order to check whether our AHF for this comparison works well, in Fig. 4 (a) we show the relation between AH area of the Kerr BH and AH radius in isotropic coordinates as a function of spin parameter. The blue line denotes the analytical AH area and red crosses denote numerical results by solving AH equation for Kerr BH. The green circles show the coordinate radii where the AHs with different are located. Much larger computational resources should be required to obtain the solution with a high BH spin because high resolution in the coordinate radius is required in this regime.
We perform numerical relativity simulations with the initial data of single puncture BH and Kerr BH in Fig. 4 (b). The BSSN evolution equations which give stable dynamical evolution [46][47][48] are adopted in these simulations. The color difference shows the difference among spins and the type of lines denotes the difference between Kerr BH and single puncture BH. The spins of single puncture BHs settle down at late time, which shows BHs relax to almost the stationary state and one can compare results of Kerr BHs at late time. The single puncture BH with the higher spin does not reach at the spin which we expect. This is because we assume the conformal flatness for constructing puncture BH but Kerr BH should not be expressed by the conformal flat metric. However, it should be noted that the puncture BH can represent the small spin BH well and it is actually powerful to construct the initial data for multi BHs system.

Conclusions
In these notes, we showed how to prepare the initial data for numerical relativity and how to obtain the apparent horizon of BHs, which are reduced to solving elliptic PDEs in general. We presented several BH solutions as initial data for numerical relativity and described several numerical methods to solve elliptic PDEs. In particular, sample codes to solve Poisson equations with linear and non-linear sources are available online to public users. It is worth noting that these simple, "classical" methods are still powerful enough to be of use for current problems. In addition, we note that modern methods (e.g. Multi-Grid method) can help to eventually upgrade these classical methods in terms of numerical costs and consuming-time. Of course, one should carefully choose which method to use to solve elliptic PDEs, according to the problem at hand.

Acknowledgments
The author would like to thank Vitor Cardoso for giving the opportunity to lecture on this school, and to the Organizers and Editors of the NR/HEP2: Spring School at Instituto Superior Técnico in Lisbon. The author would also thank an anonymous referee for a careful reading of the manuscript and many useful suggestions. The author is thankful to Ana Sousa who helps to improve English on this notes, Sérgio Almeida who maintains the cluster "Baltasar-Sete-Sóis" and Takashi Hiramatsu who maintains the "venus" cluster. Numerical computations in this work were carried out on the cluster of " Let us decompose a matrix A IJ into the lower and upper triangular matrices defined as L IJ and U IJ respectively, where D KK denotes the diagonal matrix. Then, the solution vector can be solved ψ J step by step as Thus, it is easy to obtain the solution as the following precedures. First, we obtain an auxiliary vector ξ K as solving Eq. (A.3) from ξ 1 to ξ N , Therefore, the solution vector ψ J is written by similarly solving Eq. (A.4) from ψ N to ψ 1 , As the last of this section, we note how we compute the lower and upper matrices from our matrix A IJ , which is the time-consuming part. The matrix A IJ is written with the diagonal, lower and upper triangular matrices by lower (I > J).
Thus, the components of marices are obtained in turn by Although LDU decomposition allows us to numerically solve general elliptic PDEs, the large numerical costs will be required in many cases.

Appendix B. Multi-Grid method
Multi-Grid method is proposed by R. Fedorenko and N. Bakhvalov and developed by A. Brandt. [73][74][75][76] The SOR method as mentioned in Sec. 5.2.4 has the advantage of reducing the high frequency components of residual between the exact solution and numerical solution, because the values near the grid point to be updated are used for next trial guess during the iteration. On the other hand, it would take much time to reduce the low frequency modes of redisual with this iteration method. When we consider different resolution grids, however, the low frequency modes on the finer grid can be the high frequency modes on the coarser grid. The low frequency modes of residual on the finer grid can efficiently be reduced on the coarser grid. The Multi-Grid method is based on the concept of reducing different frequency modes of residual with different resolution grids. In fact, it was implemented by some groups. [77][78][79] Appendix B.1. Multi-Grid structure Suppose we have different resolution grids and the level of different grids is labeled by k, which the larger k denotes the finer grid. One can solve the Poisson equation on the level k by any iterative methods described in Sec. 5 and obtain the numerical solution, where φ (k) is the numerical solution on the level k. We define the residual on the level k between φ (k) and the exact solution by In general, the communication of the quantities such as the residual with different grid levels is needed. Now, we just use the Lagrange interpolation to communicate with each other level defined by where F, x j , x and N denote the quantity to be interpolated, the coordinate on the level, the location to be interpolated, and the number of grid points to be used by the interpolation, respectively.

Restriction operator
After we obtain the solution on the finer grid k, we transfer the information of the solution from the finer grid k to the coarser grid k − 1. Now we use the second-order discretization scheme and choose the third-order Lagrange interpolation. We define the modified source term on the coarser level k − 1 with the information of the solution on the finer grid k by where R k−1 k denotes the restriction operator to the coarser grid k − 1 and φ (k−1) c denotes the smoothing solution by the restriction operator. Roughly speaking, the modified source term S (k−1) consists of that on the level k with smoothing operation and the correction by the difference of Laplacian operator between two levels. Then, we obtain the numerical solution φ (k−1) on the level k − 1 to solve the Poisson equation with the modified source term.

Prolongation operator
The solution with the modified source term on the coarser level k − 1 is to be brought back to the finer level k. Now the communication is also done by third-order Lagrange interpolation.
where P k k−1 denotes the prolongation operator and φ There are some ways of deciding the order of the level to compute. Fig. (5) shows the difference of such order between the methods of V-cycle and W-cycle as examples. Now we choose V-cycle because it is easier to implement to the code. We use the restriction operator before computing on the coarser level and the prolongation operator before computing on the finer level. This cycle is repeated until we obtain the expected error of the Poisson equation.   We note that the boundary of the finer grid is given by the interpolation.

Appendix C. List of Sample codes
We have some sample codes for the lecture on NR/HEP2: Spring School at Instituto Superior Técnico in Lisbon and they are available online. In this section, we show the simplest code to solve an elliptic PDE and the sample code which is parallelized with OpenMP. One can see what is the parallel computing in Ref. 80 Here is the list of sample codes which are available in http://blackholes.ist.utl.pt/nrhep2/?page=material,