{'fraction_non_alphanumeric': 0.033524345165455154, 'fraction_numerical': 0.017827770998036823, 'mean_word_length': 4.222509582967718, 'pattern_counts': {'":': 0, '<': 2, '<?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': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The initial value problem as it relates to numerical relativity

Spacetime is foliated by spatial hypersurfaces in the 3+1 split of General Relativity. The initial value problem then consists of specifying initial data for all relevant fields on one such a spatial hypersurface. These fields are the 3-metric and extrinsic curvature together with matter fields such as fluid velocity, energy density and rest mass density. There is a lot of freedom in choosing such initial data. This freedom corresponds to the physical state of the system at the initial time. At the same time the initial data have to satisfy the Hamiltonian and momentum constraint equations of General Relativity and can thus not be chosen completely freely. We discuss the conformal transverse traceless and conformal thin sandwich decompositions that are commonly used in the construction of constraint satisfying initial data. These decompositions allow us to specify certain free data that describe the physical nature of the system. The remaining metric fields are then determined by solving elliptic equations derived from the constraint equations. We describe initial data for single black holes and single neutron stars, and how we can use conformal decompositions to construct initial data for binaries made up of black holes or neutron stars. Orbiting binaries will emit gravitational radiation and thus lose energy. Since the emitted radiation tends to circularize the orbits over time, one can thus expect that the objects in a typical binary move on almost circular orbits with slowly shrinking radii. This leads us to the concept of quasi-equilibrium which essentially assumes that time derivatives are negligible in corotating coordinates, for binaries on almost circular orbits. We review how quasi-equilibrium assumptions can be used to make physically well motivated approximations that simplify the elliptic equations we have to solve.

11 Numerical codes used to construct initial data for binaries 51 12 Conclusions 52 Appendix A Lie derivatives 53 Appendix B Elliptic Equations 55

Introduction
After the first direct detection of a gravitational wave signal emitted by a binary black holes system [1,2], the problem of faithfully simulating the evolution of binary systems of compact objects has become increasingly important. The detectors (LIGO [3,4], VIRGO [5,6], GEO600 [7]) use laser interferometry to measure the strains associated with passing gravitational waves [8], and offer much higher sensitivity than previous experiments aimed at direct detection of these waves. Additional detectors such as KAGRA [9] or the space borne eLISA/NGO [10] and DECIGO [11] are in planning and construction stages. In all these detectors, the measured strains are tiny and contaminated by noise. For detection and especially parameter estimation, it is thus necessary to compare the observed signals to theoretical gravitational wave templates that describe different systems. One of the most promising sources for gravitational wave detectors is the inspiral and merger of compact objects such as black holes or neutron stars. Due to the emission of gravitational waves, the binary loses energy and the orbit tightens until finally the two objects merge. As long as the two compact objects are far apart, post-Newtonian calculations [12] can give highly accurate approximations for the orbital motion and the gravitational waves emitted by the binary. When the two objects get closer, the post-Newtonian expansion becomes more and more inaccurate and eventually breaks down. Thus the full non-linear equations of General Relativity have to be solved and computer simulations are used to obtain numerical answers [13]. From each of these computer simulations we can then extract a gravitational wave template for the simulated system.
In order to start such numerical simulations we need to specify initial data that describe the initial state of the binary system. These initial data should be chosen as accurately and as realistically as possible. Otherwise the subsequent numerical evolution may not simulate the kind of system the detectors observe.
In this paper we will concentrate on how one can construct reliable initial data. In Sec. 2 we describe the basic equations that are used in numerical relativity. This is followed in Sec. 3 by a discussion of several conformal decompositions that are used to derive the equations commonly used in the construction of initial data. Sections 4 and 5 describe various methods to construct initial data for single and binary black holes, and Sec. 6 presents possible improvements. We proceed in Sec. 7 by introducing equations that are needed to describe neutron star matter. In Secs. 8 and 9 we discuss single and binary neutron star data, while Sec. 10 deals with mixed binaries. Section 11 provides a short overview of the most common computer codes that are used to compute initial data. We provide some conclusions in Sec. 12. Throughout we will use geometric units where G = c = 1.

Numerical relativity
In this section we introduce the basic equations of numerical relativity. We start by discussing some equations of General Relativity and explain how they are reformulated for use in numerical simulations.

General Relativity
We recall here some of the basic equations of General Relativity with the purpose of illustrating our notation.
In Relativity theory 3-dimensional space and time are unified into 4-dimensional spacetime. In the simplest case of a flat spacetime distances in this 4-dimensional spacetime are measured as where x 1 , x 2 , x 3 are the usual Cartesian coordinates of Euclidean space, and x 0 is time.
Usually this distance is written as where the indices µ and ν run from 0 to 3, and we use the Einstein summation convention and sum over repeated indices. Here g µν are the components of a symmetric tensor called the spacetime metric or 4-metric. For flat spacetime in Cartesian coordinates we have g µν = diag(−1, 1, 1, 1), but if we change coordinates g µν will have different values. The inverse of the 4-metric will be denoted as g µν . Throughout we will use Greek letters to denote 4-dimensional indices that run from 0 to 3. Later when we refer to the spatial components of tensors we will use letters from the middle of the Latin alphabet. An example is the spatial 3-metric γ ij = diag(1, 1, 1) of flat space, where both i and j run from 1 to 3. We denote the 4-dimensional covariant derivative operator by ∇ α . It is defined by in order to be compatible with the 4-metric. The covariant derivative of any tensor can be given in terms of partial derivatives and the Christoffel symbols For example, the covariant derivative of a 4-vector V µ can be computed from If matter or energy is present, spacetime will become curved and the 4-metric has to be determined from the Einstein equations Here T µν is the stress-energy tensor of matter. The Ricci tensor 4 R µν = 4 R ρ µρν and the Ricci Scalar 4 R = 4 R µν g µν are both related to the Riemann tensor 4 R µ ναβ , defined by [14] 4 R ν αβµ ω ν : where ω µ is any one-form. Indices on all tensors can be raised and lowered with the spacetime metric g µν . We will follow the sign conventions of Misner, Thorne and Wheeler [15], so that the Riemann tensor when computed from the Christoffel symbols in Eq. (4) is given by As we have seen any tensor of any rank can be written using the index notation introduced above. However, sometimes it can be advantageous to use index free notation. We use u to denote a 4-vector u µ and ω to denote a one-form ω µ . The contraction of the two is then written as u · ω = u µ ω µ .

The 3+1 split of spacetime
Since both the 4-metric and the Ricci tensor are symmetric, it is clear that the Einstein equations in (6) are ten independent equations for the ten independent components of the 4-metric. From Eq. (8) we see that these ten equations are linear in the second derivatives and quadratic in the first derivatives of the 4-metric. However, only six of these ten equations contain second time derivatives. These six second order in time equations represent evolution equations. The other four are not evolution equations, but rather represent constraints that the 4-metric has to satisfy.
If we can solve the Einstein equations in (6) we obtain the 4-metric everywhere in space and time. While a solution for all of spacetime is certainly very useful, in practice it is possible to find such direct solutions of Eq. (6) only in rare cases where one assumes symmetries such as time independence and axisymmetry. In addition, in many astrophysical problems of interest it may not be necessary to find a solution for all of time. One usually would like to start from some given initial conditions at some initial time and then to calculate the state of the system at some later time. For example, for a particle in classical mechanics one starts by giving the initial position and velocity, and then calculates position and velocity at a later time. We would like to do the same for General Relativity and specify initial data (i.e. some fields and their time derivatives) at some initial time and then calculate these fields at a later time. In order to do this we use the 3+1 split of spacetime [16]: We foliate spacetime by spacelike 3-dimensional hypersurfaces or slices, i.e. surfaces with a normal vector n µ that is timelike. As time coordinate x 0 we use a function t which is constant on each hypersurface, but increases as we go from one hypersurface to the next. We then define the lapse by The normal vector thus must satisfy in order for it to be normalized such that g µν n µ n ν = −1.
The minus sign in Eq. (10) is chosen such that n µ points in the future direction of increasing t. Notice that while n µ is orthogonal to each hypersurface, n µ is not in general tangent to lines of constant spatial coordinates x i = const, since we are still free to choose arbitrary spatial coordinates on each hypersurface. Indeed the tangent vector t µ to lines of constant spatial coordinates is given by where we have normalized t µ such that The latter also implies from which it follows that where β µ is an arbitrary spatial vector, i.e. any vector that satisfies β µ n µ = 0. Once the foliation is given the lapse α and the normal vector n µ are fixed. Yet the direction of the vector t µ still depends on the choice of spatial coordinates on each hypersurface. This freedom to choose spatial coordinates is encapsulated in the spatial vector β µ = (0, β i ), which is called the shift vector.

3-metric, spatial covariant derivative and extrinsic curvature
We now define the 3-metric as which is the projection of the 4-metric onto the hypersurface with normal vector n µ . This is a proper 4-tensor whose indices can be raised and lowered with the 4-metric g µν . Since n i = 0 we find γ ij = g ij and γ i j = g i j = δ i j . From Eq. (16) we also find γ 0ν = γ µν n µ /(−α) = 0 and thus γ ik γ kj = γ iν γ νj = γ iν g νj = γ i j = δ i j . From Eqs. (12) and (15) we find that and where β j = g µj β µ = γ ij β i . Using the latter two equations, the line element in Eq. (2) can now be written as The 3-metric defined in Eq. (16) can be used to project any tensor T α 1 α 2 ... β 1 β 2 ... onto the spatial hypersurface. The resulting tensor S α 1 α 2 ...
is called a spatial tensor since it is orthogonal to n α in each index. Note that the indices on any spatial tensor can be raised or lowered by using either the 3-metric γ ij or the 4-metric g µν . The result will be the same since γ µν n ν = 0. We define the spatial covariant derivative of any spatial tensor S α 1 α 2 ...
where we have projected all free indices onto the spatial hypersurface. From Eq. (16) it immediately follows that so that the spatial covariant derivative operator D µ is compatible with the 3-metric γ αβ . This also implies that the spatial covariant derivative D i S k 1 k 2 ... l 1 l 2 ... of any spatial tensor S k 1 k 2 ... l 1 l 2 ... can be computed with the help of 3-dimensional Christoffel symbols, which can be computed from the 3-metric γ ij .
In order to describe the curvature of a t = const hypersurface one introduces the extrinsic curvature It can be shown that K µν is symmetric and also given by where £ n is the Lie derivative (see Appendix A) along the vector n µ . We can thus express K ij also as (see Appendix A.4) which tells us how time derivatives of the 3-metric are related to the extrinsic curvature. From its definition in Eq. (22) it is clear that the extrinsic curvature tells us how the normal vector n µ changes across a t = const hypersurface. It thus measures how this hypersurface is curved with respect to the 4-dimensional space it is embedded in. However, the 3-dimensional space within each hypersurface can be curved as well. This so called intrinsic curvature is described by the 3-dimensional Riemann tensor R i jkl that can be computed from the 3-metric γ ij and its resulting 3-dimensional Christoffel symbols using Eqs. (4) and (8) for γ ij instead of g µν . We will see below how the curvature of spacetime described by the 4-dimensional Riemann tensor 4 R µ ναβ is related to the extrinsic curvature K ij and the intrinsic curvature R i jkl . Also notice that in three dimensions, the Riemann tensor has the same number of independent components and carries the same information as the 3-dimensional Ricci tensor R jl = R i jil .

Einstein's equations and the 3+1 split
Going back to our example of a particle in classical mechanics, we see that the 3-metric γ ij is analogous to the particle position, while the extrinsic curvature K ij is analogous to its velocity. Of course, in order to compute the time evolution of a particle we also need the equation of motion that tells us the time derivative of the velocity. In our case we will use Einsteins equations to find additional equations that tell us how K ij evolves in time. This is done by relating the 4-dimensional Riemann tensor 4 R αβµν to the 3-dimensional Riemann tensor R αβµν with the help of the Gauss-Codazzi and Codazzi-Mainardi relations given by (see e.g. [17] or [18]) and One can also show that By looking at the different projections of the Einstein equations (6) onto n µ and γ µν and using Eqs. (25), (26), (27) and (7) we find that Einstein's equations split into the evolution equations and the so called Hamiltonian and momentum constraint equations Here R ij and R are the Ricci tensor and scalar computed from γ ij , D i is the derivative operator compatible with γ ij and all indices here are raised and lowered with the 3-metric γ ij . The source terms ρ, j i , S ij and S = γ ij S ij are projections of the stress-energy tensor T µν given by and correspond to the energy density, flux and stress-tensor seen by an observer moving with 4-velocity n µ .

Constraints, gauge freedom and physical degrees of freedom
The evolution equations in Eqs. (28) and (29) tell us how we can obtain γ ij and K ij at any time, if we are given γ ij and K ij at some initial time. As one can see, the constraint equations (30) and (31) do not contain any time derivatives and thus are equations that need to be satisfied on each spatial hypersurface. Thus when we construct initial data, we have to ensure that both γ ij and K ij satisfy the constraint equations at the initial time. It can then be shown that the evolution equations will preserve the constraints. For this reason constructing initial data is a non-trivial task in the sense that we cannot freely choose the 12 fields γ ij and K ij . Rather these fields are subject to the 4 constraints in Eqs. (30) and (31). This leaves us with only 8 freely specifiable fields. However, in General Relativity we are free to choose any coordinates we like. Since there are 4 coordinates that we can freely choose, the 8 freely specifiable fields in γ ij and K ij can really contain only 4 physical fields that are independent of our coordinate choice. Since we specify both γ ij and its time derivative K ij (see Eq. (24)) 2 of these 4 fields are merely time derivatives of the other 2 fields. Thus we are dealing with only 2 physical degrees of freedom, as expected for gravity. The freedom to choose coordinates is usually referred to as gauge freedom and connected to the choice of lapse α shift β i . If we only change the 3 spatial coordinates, all components of γ ij and K ij as well as the shift will change in general, but the spatial hypersurfaces themselves are unchanged, so that the normal vector n µ and the lapse α will not change. If we change the time coordinate t, the hypersurfaces themselves will also change so that now n µ and the lapse α will change as well [19].

Conformal decompositions and initial data construction
As already mentioned, the initial data γ ij and K ij cannot be freely chosen, since they must satisfy the constraint equations (30) and (31). However, since the constraint equations alone cannot determine the initial data either, Eqs. (30) and (31) do not have unique solutions and are thus not directly useful in constructing initial data. For this reason so called conformal decompositions have been developed. Using these decompositions it is possible to start by specifying arbitrary γ ij and K ij that may not satisfy the constraints. In a second step these γ ij and K ij are then modified in such a way that they satisfy the constraints. This second step involves solving certain elliptic equations for 4 auxiliary quantities that are used to compute the final constraint satisfying γ ij and K ij . These elliptic equations have the great advantage that they generally have unique solutions once appropriate boundary conditions are specified (see Appendix B). Except for simple cases with many symmetries these elliptic equations are usually solved numerically [20].
Following Lichnerowicz [21] and York [22,23] we start by decomposing the 3-metric γ ij into a conformal factor ψ and a conformal metricγ ij such that Using this decomposition the Ricci scalar can be written as HereD k is the derivative operator compatible with the conformal metricγ ij andR the Ricci scalar computed fromγ ij . It is also possible to verify that for any symmetric and tracefree tensor M ij we have [24] Let us also introduce a symmetric tracefree differential operator L. It is defined by its action on any vector W i by It is sometimes referred to as conformal Killing form since it is related to the Lie derivative of the tensor density γ −1/3 γ ij of weight −2/3, where γ = det(γ ij ). Using Eq. (A.12) we find So if W i is a Killing vector of the metric γ −1/3 γ ij with £ W (γ −1/3 γ ij ) = 0, we will have (LW ) ij = 0. Notice that in general γ −1/3 γ ij is not the same as the conformal metricγ ij . One can show that ijD k W k is computed using the conformal metric. Next, the extrinsic curvature is split into its trace K = γ ij K ij and its tracefree part A ij by writing it as The tracefree part is rescaled as Note that the factor of ψ −10 in Eq. (41) has been picked so that Eq. (35) applies to A ij .
Inserting Eqs. (34), (40) and (41) into the Hamiltonian constraint (30) yields where indices on the barred quantities are raised and lowered withγ ij . We will use this equation for the Hamiltonian constraint to compute the conformal factor in the different conformal decompositions discussed below. Note, however, that these decompositions differ in howĀ ij is split further.

The conformal transverse traceless decomposition
We consider first the conformal transverse traceless (CTT) decomposition, where we split [25]Ā in two pieces. HereM ij is tracefree but otherwise arbitrary andσ is some positive weighting factor that we can choose. The vector W i will be computed below.
The Hamiltonian constraint is written as in Eq. (42) withĀ ij defined as in Eq. (43). When we insert Eqs. (40), (41) and (43) into the momentum constraint (31) we find Both Eqs. (42) and (44) are elliptic equations. The weighting factor is often simply set toσ = 1. In this case we obtain the standard CTT decomposition often called the York-Lichnerowicz decomposition [21,22,23]. However, other choices are possible. For example, forσ = ψ −6 we obtain the so called physical CTT decomposition [20]. Let us now discuss how Eqs. (42) and (44) are used in practice. We start with some reasonable guess or approximation for the physical situation we want to describe. For example, if we want to find initial data for two black holes we could simply take the superposition of the 3-metrics and extrinsic curvatures for two single black hole solutions ‡. Since the superposition principle does not hold in a non-linear theory like General Relativity, this superposition will at best be an approximate solution that we denote here by γ approx ij and K ij approx . Hence γ approx ij and K ij approx will not satisfy the Hamiltonian and momentum constraints. However, using the above CTT decomposition we can now construct a γ ij and K ij that will satisfy the constraints if we set and then solve Eqs. (42) and (44) for ψ and W i . Once we have these solutions we can use Eq. (33) and Eqs. (40), (41) and (43) to obtain a γ ij and K ij that are now guaranteed to satisfy the Hamiltonian and momentum constraints.

The conformal thin sandwich approach
A form of the conformal thin sandwich formalism (called the Wilson-Mathews approach [26,27]) was first introduced under the assumption that the conformal metric is flat, i.e.γ ij = δ ij . Here we relax this assumption and present the more general conformal thin sandwich (CTS) decomposition proposed by York [28]. Like the CTT decomposition, it starts again with the conformal metric decomposition in Eq. (33). However, in order to have a direct handle on time derivatives one also introduces and demands that i.e. that the time evolution of conformal factor be chosen such that the determinant of the conformal metricγ is instantaneously constant. The latter condition yields where γ = det(γ ij ). Together with ∂ t ln γ = γ ij ∂ t γ ij , Eqs. (24) and (40) we obtain Using Eq. (39) for W i = β i , Eq. (49) can be rewritten as where we have defined Inserting Eqs. (34), (40) and (50) into the Hamiltonian constraint (30) yields again Eq. (42), but this time withĀ ij defined according Eq. (50). When we insert Eqs. (40) and (50) into the momentum constraint (31) we find From Eqs. (50) and (52) we can see that the CTS equations can be obtained from the CTT equations by settinḡ σ = 2ᾱ, Up to this point the CTS approach may not seem to differ very much from the conformal transverse traceless approach. Once we specifyγ ij ,ū ij , K andᾱ, we can solve Eqs. (42) and (52) for ψ and β i and then use Eq. (33) and Eqs. (40) and (50) to obtain a γ ij and K ij that are now guaranteed to satisfy the Hamiltonian and momentum constraints. Notice, however, that we do not just obtain initial data for the 3-metric and extrinsic curvature in this way. We also obtain a preferred lapse α = ψ 6ᾱ and shift β i .
So farᾱ was an arbitrarily chosen function. It is possible to relate it to the time derivative of K. From Eqs. (29) and (30) we find that Thus the lapse α is related to ∂ t K. Using Eqs. (51) and (42), Eq. (54) can be rewritten as an elliptic equation forᾱ. We find With the latter it is now possible to specify ∂ t K instead ofᾱ. The fact that we can directly specify the time derivativesū ij and ∂ t K is usually exploited when one wants to construct initial data in quasi-equilibrium situations, where a coordinate system exists in which the time derivatives of 3-metric and extrinsic curvature either vanish or are small in some way, as we will discuss below.

Boundary conditions at spatial infinity
The conformal decompositions discussed above result in elliptic equations that require boundary conditions. In the case of the CTT decomposition we have to solve Eqs. (42) and (44) for ψ and W i . This is usually done in asymptotically inertial coordinates, i.e. coordinates such that the 4-metric approaches the Minkowski metric diag(−1, 1, 1, 1) at spatial infinity which we denote by r → ∞. In this case we obtain Usually we also choose so that Eq. (33) implies To ensure lim r→∞ K ij = 0 we usually choose which together with Eqs. (40), (41) and (43) yields Once Eqs. (42) and (44) are supplemented by the boundary conditions (59) and (61) at infinity we will obtain unique solutions for ψ and W i . In the case of the CTS decomposition we have to solve Eqs. (42), (52) and (55) for ψ, β i andᾱ. If we again work in asymptotically inertial coordinates which result in Eq. (57) at spatial infinity, we find again Eq. (59), as well as lim r→∞ᾱ = 1 (62) and However, sometimes it is convenient to work in a corotating frame, which is obtained by changing the spatial coordinates such that they rotate with a constant angular velocity Ω with respect to the asymptotically inertial coordinates at infinity that are used in Eq. (57). When we change coordinates to this corotating frame the 4-metric components change such that now Here Φ i is an asymptotic rotational Killing vector at spatial infinity. Hence the boundary conditions for ψ andᾱ are unchanged, while the boundary condition for the shift becomes Together with these boundary conditions the elliptic equations (42), (52) and (55) have unique solutions, provided no other boundaries are present. Notice that the shift condition (65) can be problematic for some numerical codes since β i β i → Ω 2 r 2 as r → ∞, so that the shift becomes infinite at spatial infinity. This problem is usually avoided by replacing the shift β i by in equations such as (50) and (52). These equations contain the shift only in the form (Lβ) ij , and since whenγ ij is flat, Eqs. (50) and (52) still have the same form, only with β i replaced by B i . The boundary condition for B i then is simply

ADM quantities at spatial infinity and the Komar integral
The ADM mass M ADM ∞ and angular momentum J ADM ∞ at spatial infinity (r → ∞) are given by and Here Φ i is again the asymptotic rotational Killing vector at spatial infinity. Note that these definitions are valid only in coordinates where the 4-metric approaches the Minkowski metric at spatial infinity. In this case lower and upper spatial indices do not need to be distinguished, so that summation is implied over any repeated spatial indices. The integrals here are over a closed sphere (denoted by S ∞ ) at spatial infinity.
In spherical coordinates dS j is given by 4πr 2 (sin θ cos φ, sin θ sin φ, cos θ)dθdφ. M ADM If the spacetime possess a Killing vector ξ µ , i.e. a vector satisfying then the Komar integral [29,30,14] integrated over any closed 3-surface S containing the sources yields a value that is independent of S. For a time-like Killing vector that is normalized to ξ µ ξ µ = −1 at spatial infinity I K (ξ, S) = M ADM ∞ , while for an axial Killing vector that approaches the asymptotic rotational Killing vector Φ µ at spatial infinity I K (ξ, S) = −2J ADM ∞ [30]. We now consider a helical Killing vector ξ µ that is given by at spatial infinity, where T µ and Φ µ are the asymptotic time-translation and rotational Killing vectors at spatial infinity and Ω is a constant. If ξ a is normalized such that lim r→∞ ξ a n a = −1, then the Komar integral over a sphere at r → ∞ becomes [30,31] (78) § Analogously to the ADM momentum P ADM i ∞ , it is also possible to compute the components of ADM angular momentum from J ADM

Quasi-equilibrium assumptions for the metric variables of binary systems
We now make some additional simplifying assumptions to better deal with quasiequilibrium situations. For concreteness, let us discuss a binary system made up of two orbiting objects that could be stars or black holes. Such a system will radiate gravitational waves that carry away energy. Thus over time the two objects will spiral towards each other. According to post-Newtonian predictions [32,33] the orbits circularize on a timescale which is much shorter than the inspiral timescale due to the emission of gravitational waves. The same circularization effect has also been shown in the extreme mass ratio limit in Kerr spacetime [34,35,36]. Hence when we set up initial data we often assume that the binary is in an approximately circular orbit. The circular orbital motion then can be transformed away by changing to a corotating coordinate system, i.e. a coordinate system that rotates with the orbital angular velocity of the binary with respect to asymptotically inertial coordinates. In this corotating system the two objects do not move, which means ∂ t g µν vanishes at least approximately. This implies that the time evolution vector t µ is an approximate symmetry of the spacetime, which in turn implies the existence of an approximate helical Killing vector ξ µ (see e.g. [31,37]) with Here £ ξ is the Lie derivative (see Appendix A.3) with respect to the vector ξ µ . Corotating coordinates are understood as the coordinates chosen such that t µ lies along ξ µ , such that Since we are considering orbiting binaries, this Killing vector would have to be a helical Killing vector, which means that at spatial infinity (r → ∞) where T µ and Φ µ are the asymptotic time-translation and rotational Killing vectors at spatial infinity and Ω is the angular velocity with which the binary rotates. In asymptotically inertial coordinates we have , where we have chosen the rotational Killing vector to correspond to a rotation about the a line through the center of mass x i CM parallel to x 3 -axis. An approximate helical Killing vector with £ ξ g µν ≈ 0 implies that the Lie derivatives with respect to ξ µ of all 3+1 quantities similarly vanish approximately. For the initial data construction using the CTS decomposition we will use In a corotating coordinate system where the time evolution vector t µ lies along ξ µ , the time derivatives of these metric variables are then equal to zero. From ∂ tγij = 0 it follows thatū ij in Eq. (50) vanishes. The assumption ∂ t K = 0 can be used in Eq. (55) to findᾱ.
As we have seen above, the approximate helical Killing vector has the form ξ µ = αT µ + ΩΦ µ at spatial infinity, where Ω is the orbital angular velocity of the binary. It is clear that circular orbits are only possible for one particular value of Ω, for a given distance between the two components of the binary. We will now discuss a method that can be used to find this Ω. If we insert ξ a as given in Eq. (80) into Eq. (75) we obtain Using Eq. (81) we find that the second term integrated at r → ∞ is Therefore combining Eqs. (78), (83) and (84) the condition must hold if the helical Killing vector of Eq. (80) exists. This condition can then be used to fix the value of Ω. The term I K (αn, S ∞ ) is often called the Komar mass. For asymptotically flat stationary spacetimes the equality of Komar and ADM mass has already been shown in [38]. It is closely related to the virial theorem in general relativity [39] and was first used in [40,41] to determine the orbital angular velocity Ω for binary black hole initial data.

Apparent horizons of black holes
As is well known, black holes are enclosed in so called event horizons. An event horizon is defined as the boundary of a spacetime region from where nothing can escape to spatial infinity. Thus in order to find the event horizon we need to know the entire future spacetime. Hence this definition is not very useful when we want to construct initial data on a single spatial slice. For this reason the concept of apparent horizons has been introduced. One considers outgoing photons and considers the spatial 2-surface S that is given by the location of these outgoing photons at time t. In flat spacetime this surface will expand, so that its area increases (if the photons are outgoing). An apparent horizon is defined as the outermost closed surface where this expansion is zero. This concept captures the idea of no escape from a black hole, but this time the definition involves only a single spatial slice. Let s µ be the outward pointing spatial normal vector (with s µ n µ = 0) of the surface S, normalized to The induced metric on S is then given by If we denote its determinant by h, the surface area of S is given by integrating √ h over S. We now define the outgoing lightlike vector The expansion can then be defined as Using The apparent horizon is now defined as a surface on which Θ = 0. In order to find this surface it is convenient to define it as the location where a function F (x i ) = 0. The normal vector then becomes Inserting into Eq. (90) we obtain as the equation that determines the location of the apparent horizon [17]. A review on algorithms for finding apparent horizons can be found in [42]. Given a spacetime with an event horizon and given a particular slicing, we can try to find the apparent horizon on each slice by solving Eq. (92). When they exist, the apparent horizons at different times t form a world tube through spacetime, that will either be inside or coincide with the event horizon [43].

Quasi-equilibrium boundary conditions at black hole horizons
As we will see below (when we discuss puncture initial data) it is possible to construct black hole initial data without any inner boundaries at or inside the black hole horizons. However, in many cases such as e.g. Misner or Bowen-York, or superimposed Kerr-Schild initial data one needs a boundary condition at each black hole as well. Here we will discuss quasi-equilibrium boundary conditions introduced by Cook et al. [44,45,46] that can be imposed at the apparent horizon of each black hole within the CTS approach. As clarified in [47], the assumptions required for these quasi-equilibrium boundary conditions are essentially the same as those required of an isolated horizon [48,49,50]. One starts by picking some surface S (e.g. a coordinate sphere) that one wishes to make into an apparent horizon in equilibrium. To ensure that S is indeed an apparent horizon one imposes Since this horizon is supposed to be in equilibrium, one further demands that the shear of outgoing null rays on S should be [50,47] σ µν | S = 0.
Finally, coordinates adapted to equilibrium should be chosen such that the apparent horizon does not move as we evolve from the initial slice at t to the next one at t + dt. A photon moving along the outgoing light ray given by l µ of Eq. (88) moves along the curve dx µ dλ = l µ , where λ is the parameter along the curve. This curve intersects the slice at t + dt for dx 0 = l 0 dλ = dt, and thus at n 0 dλ = dt, which yields dλ = αdt. So in the time interval dt a photon moves in the direction αl µ while a point with constant spatial coordinates moves in the direction t µ . We could now choose our coordinates such that t µ − αl µ = 0, which would ensure that all photons that make up the surface S stay at the same spatial coordinate. However, for the surface S to remain unchanged it is sufficient when each photon only moves along the surface. So a less restrictive condition is to demand that which implies that this difference has no component perpendicular to S. It is equivalent to and is interpreted as a boundary condition on the component of the shift perpendicular to S. Following [46] we now introduce a conformally rescaled normal vector and surface metric by defining Combining Eqs. (90), (93) and the CTS definitions for K ij in Eqs. (40), (41) and (50) we arrive at [44,45,51] This is a Robin-type boundary condition on the conformal factor, derived under the assumption (93) of vanishing expansion, that defines the apparent horizon. So far we have boundary conditions only for the conformal factor and the component of the shift perpendicular to S. If we define the parallel component by it can be shown [46] that the vanishing shear assumption (94) leads to whered i is the derivative operator compatible with the 2-metrich ij . So ifū ij = 0, β i is a conformal Killing vector of the metrich ij on S. In practice one uses as boundary condition for the shift. Here χ i is a rotational conformal Killing vector on h ij with affine length of 2π and with χ i s i = 0. The constant Ω r is can be chosen to change the spin of the black hole.
Finally, we also need a boundary condition forᾱ. Several conditions will work [46]. The simplest choice is to use a Dirichlet boundary condition where one prescribes a value for αψ on S.

Black hole initial data for single black holes
We will now discuss how one can construct initial data for single black holes. We will present several analytic black hole solutions to the Einstein equations. These solutions represent single stationary black holes. Since they are obtained by solving the Einstein equations they automatically satisfy the constraints.

Non-spinning black holes
A non-spinning spherically symmetric black hole of mass m can be described by the Schwarzschild metric where 0 ≤ r s ≤ ∞. The event horizon is located at r s = 2m. The region r s ≤ 2m is inside the black and at r s = 0 there is a physical singularity where curvature invariants are infinite. However, this metric is rarely used in numerical relativity because it has a coordinate singularity at r s = 2m. We can change the radial coordinate by setting where we have defined From this we obtain the metric in isotropic coordinates where 0 ≤ r ≤ ∞. As we can see the 3-metric is now conformally flat. The event horizon is at r = m/2. Notice however, that the new r coordinate does not cover the inside of the black hole since r s ≥ 2m for any r. For r → ∞ the metric becomes flat. Indeed, there is another asymptotically flat region at r = 0. This can be seen by performing the coordinate transformation which maps the region with r < m/2 into r > m/2. As one can easily verify, the 3-metric in Eq. (106) remains conformally flat under this transformation. Its new conformal factor is ψ = 1 + m 2r , i.e. it is invariant under this transformation. The transformation in Eq. (107) is thus an isometry. This shows that at r = 0 (i.e. r → ∞) we obtain another asymptotically flat region.
One can also switch to Cartesian spatial coordinates, which yields where now r = x 2 + y 2 + z 2 . If we apply the 3+1 split we find which is often used as initial data for a single non-spinning black hole.

Black holes with spin
In order to describe a spinning black hole one can use the Kerr metric that is also an exact solution to Einstein's equations. The 4-metric in so called Kerr-Schild coordinates (see e.g. [15]) is given by where we the x µ stand for (t, x, y, z), η µν is the Minkowski metric, and Here m is the mass and a is the spin parameter which takes values in the range 0 ≤ a ≤ m. The function r here depends on (x, y, z) and is given implicitly by The event horizon is located at r = m + √ m 2 − a 2 . The physical curvature singularity occurs on a ring given by x 2 +y 2 = a 2 and z = 0, which is inside the event horizon. If a is non-zero, no spatial slicing exists in which the 3-metric can be written in a conformally flat way [52]. The Kerr-Schild metric in 3+1 form is given by [53] while the extrinsic curvature can be obtained by using the derivative operator D i compatible with the 3-metric. The Kerr-Schild coordinates described here are horizon penetrating in the sense that there is no coordinate singularity in γ ij and K ij at the horizon, and that they cover both the inside and outside of the black hole.

Initial data for black hole binaries
The binary black hole problem is more complicated and no analytic solutions to the Einstein equations are known. Here we will discuss several approaches to set up initial for a black hole binary. We do not attempt to discuss every possible method, rather we will focus on some of the most widely used approaches.

Black holes at rest
It is possible to find analytic initial data for N black holes at rest. We start by observing that satisfies the momentum constraint (31) for the case of vacuum where ρ = 0 = j i . We further choose the 3-metric to be conformally flat, i.e. we set for the conformal metric of Eq. (33). Thus the Hamiltonian constraint (42) simplifies tō whereD i = ∂ i is the flat space derivative operator. A solution that satisfies the boundary condition (59) is given by Here A labels the black holes, so that m A is the mass of black hole A, and is the conformal distance from the black hole center located at (x, y, z) = (c x A , c y A , c z A ). These initial data are known as Brill-Lindquist initial data [54,55] and result in N black holes at rest that will fall toward each other if evolved forward in time. They are a straightforward generalization of the Schwarzschild solution in isotropic coordinates presented above. Notice that the conformal factor in Eq.
there are singular points that have to be removed from the manifold. The total ADM mass of these data is given by The main difference between Brill-Lindquist initial data and the Schwarzschild metric in isotropic coordinates is that we now have N + 1 different asymptotically flat regions. One is located at r = ∞, the others are at r A = 0 for A = 1, ..., N . An observer sitting near r = ∞ will observe N black holes, while an observer sitting in one of the other asymptotically flat regions will see only one black horizon. Due to this asymmetry the asymptotically flat regions are not isometric to each other as was the case for a single Schwarzschild black hole. However, an approach due to Misner [56] allows us to construct another solution to Eq. (118) where the 3-metric has two isometric asymptotically flat hypersurfaces connected by N black holes. The case of two black holes that satisfies this isometry condition is known as Misner initial data. This solution is written in terms of an infinite series expansion. As in the case of Brill-Lindquist data there are singular points, but this time there are an infinite number of singular points for each black hole.

Puncture initial data
It is possible to generalize Brill-Lindquist data to an arbitrary number of moving black holes that have non-vanishing K ij . We present here so called puncture initial data which is one of the main initial data types used in numerical relativity. We start again with the conformally flat metric in Eq. (117), but this time we only assume Furthermore, we use the CTT decomposition withσ = 1 and choosē in Eq. (43). Then Eq. (44) simplifies tō A solution of this equation is [57] Since the momentum constraint in Eqs. (31) and (44) is a linear equation, linear superpositions of several such solutions will again be a solution. We can thus generalize Brill-Lindquist data by writinḡ where again A labels the black holes, and c i A , P i A and S i A are constant vectors. One can show that the ADM momentum and angular momentum for this extrinsic curvature are given by at r = ∞. Thus the vectors P i A and S i A can be interpreted as momentum and spin of black hole A.
The next step is to solve the Hamiltonian constraint (42). As we have seen in Eq. (119), Brill-Lindquist data are singular at each black hole. Inserting the ansatz into Eq. (42) we obtain where we have introduced As we can see from Eq. (119), Therefore we expect u to be regular at r A = 0. Indeed, Brandt and Brügmann [58] have shown that Eq. (129) can be solved on R 3 without any points removed. Then u will be at least C 2 at r A = 0. This method is known as the puncture approach. It constitutes a significant simplification since it obviates the need for any boundary conditions at r A = 0, that would need to be specified if we were to solve on This approach is widely used by the numerical relativity community. It allows us to set up N black holes with both momentum and spin and like Brill-Lindquist data it results in a solution with N + 1 asymptotically flat regions.
Note that the ADM mass depends on u and thus is not equal to the sum of the mass parameters m A , which therefore are referred to as bare masses that do not have a direct physical meaning.

Bowen-York initial data
Puncture initial data do not have an isometry. But Bowen and York [57] have shown that the solution to the momentum constraint in Eq. (124) can be generalized tō TheĀ ij ± represent two inversion symmetric solutions, i.e. solutions that are isometric under the transformation r = b 2 /r . Here r = b is the radius that constitutes the fixed-point set of the isometry. Since the momentum constraint is linear we can again consider superpositions of thisĀ ij ± to generate multi hole solutions. The process of constructing momentum constraint solutions with isometry fromĀ ij ± is rather complex and again results in an infinite-series sum. However, it is possible to evaluate this sum numerically [59]. With the solution forĀ ij at hand, the Hamiltonian constraint (42) becomesD which is an elliptic equation that has to be solved numerically. The boundary condition at infinity is again Eq. (59). This time, however, we also need a boundary condition at each black hole. It is derived by demanding isometry at each black hole and reads It needs to be imposed at the surfaces given by r A = b A which constitutes the fixedpoint set of the isometry at each hole. Here N i A is the normal to the surface r A = b A . The presence of these inner boundaries complicates the numerical method, when we compare with the puncture approach. For this reason these initial data are rarely used in practice. Furthermore, in many applications one is interested only in the spacetime outside the black holes, so that it is not important whether the initial data indeed possess an isometry.

Superimposed Kerr-Schild initial data
As we have seen above, both puncture initial data as well as Bowen-York initial data are generalizations of the Schwarzschild solution in isotropic coordinates. Both use a conformally flat 3-metric and make use of analytic solutions to the momentum constraints, the so called Bowen-York extrinsic curvature given in Eqs. (124) and (131). Let us consider a single black hole with spin only, i.e. the extrinsic curvature is given by Eq. (124) with c i = 0 = P i and a non-zero S i . Once we solve the Hamiltonian constraint for the conformal factor ψ, we obtain valid initial data. These data, however, cannot be a spatial slice of the stationary Kerr spacetime, since no spatial slicing of Kerr exists with a conformally flat 3-metric [52,60,61]. This means that these initial data do not describe a stationary situation. Rather, we obtain a black hole that is surrounded by gravitational waves, of which some will move out and carry away some of the angular momentum. In fact, if we increase the magnitude of S i , the ADM mass also increases so that the dimensionless spin ratio |S i |/M 2 ADM does not increase above about 0.928 [62,63]. Hence we cannot come close to the extreme Kerr limit of |S i |/M 2 ADM = 1. Analogous effects occur if we consider a single moving black hole with S i = 0 and a non-zero P i . Again the data contain gravitational waves that themselves will contain some of the momentum.
For this reason another approach has been developed [53,64] where one starts from the single spinning Kerr solution given in Eq. (110). In order to obtain a moving black hole, the 4-metric of Eq. (110) can be boosted. Due to its special structure the 4-metric in Eq. (110) is invariant under boosts. After the boost it is given by where H and k µ can be obtained by transforming the expressions in Eqs. (111) and (112) like a scalar and a vector respectively. Since Hk µ k ν in the 4-metric drops off away from the black hole, a superposition of two boosted black holes of the form will be an approximate solution to Einstein's equations, that does not satisfy the constraints exactly. This approximation can then be used, as explained after Eq. (45), to generate constraint satisfying initial data. This idea has first been used within the context of the CTT decomposition [64]. More recently it has been also used within the CTS approach to generate very highly spinning black hole binaries [63]. Since all time derivatives are assumed to be zero in corotating coordinates, one has to specify only the conformal 3-metric and the trace of the extrinsic curvature. They are then superposed according toγ where γ A ij and K A are the the conformal 3-metric and the trace of the extrinsic curvature of the boosted spinning Kerr-Schild black hole A. Here e −r 2 A /w 2 A is a Gaussian weighting factor that depends on a length scale w A and the conformal distance r A from black hole A. This factor is similar to the attenuation factors in [64,65]. It ensures that near each black hole we have a metric that is very close to a single black hole solution. The length scale w A is chosen to be larger than the size of black hole A, but smaller than the distance to the nearest neighboring black hole.
As boundary conditions [63] uses Eqs. (59), (62) and (65) at infinity. The apparent horizon S A of black hole A is chosen to be the coordinate location of the event horizon of the boosted spinning black hole A in the superposition. For the conformal factor and the shift Eqs. (99) and (102) are used as boundary conditions on S A . The lapse boundary condition there is given by With this approach it is possible to construct initial data with very highly spinning black holes. In [63] a dimensionless spin of 0.9997 is achieved. Further improvements for such data constructed within the CTS approach are presented in [66].

Conformally flat CTS binary black hole data
The CTS approach can also be used with conformally flat initial data [45,46,51,67]. In this case we assume all time derivatives are zero and simply set As boundary conditions at infinity one uses Eqs. (59), (62) and (65). The fact that we are dealing with black holes is incorporated into the initial data by specifying two conformal coordinate spheres S A where one imposes the apparent horizon boundary conditions (99) and (102) for the conformal factor and the shift. Usually one uses dαψ for the lapse boundary condition at the apparent horizons. Such conformally flat data are easier to construct than superimposed Kerr-Schild type initial data. But they suffer from the same problems as puncture initial data, in that one cannot generate highly spinning black holes.

Excision versus puncture approaches
When we construct puncture initial data, we have no inner boundary conditions and thus obtain data everywhere. However, if we use the apparent horizon boundary conditions (99) and (102) to construct black holes, in either the superimposed Kerr-Schild or the conformally flat CTS approach, we obviously obtain initial data only outside the apparent horizon. The lack of data inside the black holes leads to complications when we evolve such data numerically. First of all, one needs to use a more complicated topology for the numerical domain, where two spheres have to be excised. This approach is called black hole excision. Once this is done, one might think that it should be possible to evolve only the outside of the black holes, since nothing physical can come out of a black hole horizon. This means we should impose no boundary condition on any physical field or mode at the horizon, since no physical mode can enter the numerical domain outside the horizon. This approach works well for evolution systems such as the generalized harmonic evolution system [68,69] that have no superluminal modes, when using their standard gauges. We note however, that the often used BSSNOK system [70] as well as the Z4c system [71,72] possess superluminal gauge modes that can cross the horizons, when the usual moving puncture formalism [73,74] is used. Thus one needs boundary conditions for some modes, which further complicates the numerics. For this reason, one sometimes fills the black hole interiors [75,76,77,78] with approximate data when using the BSSNOK system. Both black hole excision and black hole filling techniques can work. Nevertheless puncture initial data are easier to evolve with the BSSNOK system using the moving puncture formalism, since neither excision nor filling are then required. For evolutions with the generalized harmonic evolution system, black hole excision is required even for puncture initial data, since the system in its standard form cannot handle the divergences in the conformal factor and the Bowen-York extrinsic curvature. Notice, however, that with a recent modification of the generalized harmonic system, where the 3-metric is conformally rescaled as in the BSSNOK system, it has been possible to evolve Brill-Lindquist initial data without excision [79].

Quasi-circular orbits
Above we have discussed the examples of puncture and superimposed Kerr-Schild initial data. In both cases we can freely choose the momentum of both black holes. As we have discussed already in Sec. 3.5, we expect most binaries to be in circular orbits with a slowly shrinking radius. The question thus arises how one should pick the black hole momenta for such a quasi-circular configurations.
For punctures an answer has been given in [31,37]. We construct puncture initial data (using the standard CTT decomposition) for two black holes. We choose the center of mass to be at rest and let P i 1 and P i 2 be perpendicular to the line connecting the two black holes. The momenta are then characterized by the parameter P = |P i 1 | = |P i 2 |. It is assumed that an approximate helical Killing vector ξ µ satisfying Eq. (79) exists for the correct choice of P . We choose coordinates, i.e. a lapse and shift such that the time evolution vector t µ = ξ µ . In these coordinates all time derivatives should approximately be zero. We then use ∂ t K = 0 to turn the evolution equation (54) into an elliptic equation for the lapse. We find Eq. (56) with K = 0 =R and ρ = 0 = S for punctures in vacuum. This elliptic equation can be solved for α, with the puncture ansatz where v is a function that is regular at r A = 0. Then Eq. (56) becomes The value of the lapse needs to be specified at r → ∞ and in the other two asymptotically flat regions at r 1 = 0 and r 2 = 0. For r → ∞ we use α = 1. The other two values α| r 1 =0 = −c 1 and α| r 2 =0 = −c 2 are adjusted together with P until the equality of Komar integral and ADM mass of Eq. (85) is satisfied in all three asymptotically flat regions. This yields a unique P . For initial data constructed within the CTS decomposition it is somewhat easier to find the black hole momenta corresponding to the quasi-circular orbits, since we can directly set certain time derivatives to zero, as required when the time evolution vector t µ is chosen to lie along an approximate helical Killing vector ξ µ . To then solve the CTS equations for say the superimposed Kerr-Schild initial data, we need to specify the value of Ω in the shift boundary condition (65). We can then simply adjust Ω (as well as the corresponding boost velocities of the individual Kerr-Schild holes) until Eq. (85) is satisfied.
When we use Eq. (85) to choose the black hole momenta, we obtain quasi-circular initial data. These initial data can then be evolved to study the orbits of the black hole binaries. For both punctures and the superimposed Kerr-Schild initial data constructed within the CTS approach, we find that the orbits are reasonably circular inspiral orbits. The solid line in Fig. 1 shows the evolution of the coordinate distance between the two black holes when we start from puncture initial data and choose the initial momenta perpendicular to the line connecting the two black holes, and use the equality of Komar integral and ADM mass (discussed above) to determine the magnitude of the initial black hole momentum parameter. As we can see the distance between the two black holes is not monotonically decreasing during the inspiral due to some residual eccentricity. This can be explained by the fact that in this case the initial data are constructed in an approximation that assumes that the orbits are exactly circular, i.e. both black holes have only tangential momentum. Yet, in a real inspiral the black holes also have a small inward momentum component. This means the initial momenta are not the ones needed for a realistic inspiral situation. The situation can be improved if we also allow . This plot shows the coordinate distance r for two different binary black hole simulations versus time. In both cases the mass ratio is 3 and the spins are pointing along the orbital angular momentum with dimensionless spin magnitude of 0.4 for the smaller and 0.6 for the larger black hole. Both simulations are starting from puncture initial data, but with slightly different initial separations and black hole momenta. The solid line corresponds to the quasi-circular setup where we have only tangential momentum, chosen by the equality of Komar integral and ADM mass. The broken line is for data where we have added a small radial momentum and tuned the tangential momentum to obtain less eccentric orbits [80].
for a radial momentum component and adjust the tangential momentum magnitude. When this is done according to the prescription described in [80] we obtain an inspiral with much less eccentricity, as shown by the broken line in Fig. 1. By using such methods [81,82,83,84,85,80] it is possible to achieve eccentricities on the order of 10 −4 for puncture initial data. As already mentioned, for initial data constructed within the CTS approach we also find eccentric orbits comparable to the solid line in Fig. 1, if we assume a helical Killing vector and construct orbits with only tangential velocity. The reason is again that a true inspiral does not have exactly a helical symmetry. However, as described in [51,86] one can replace this helical symmetry by an inspiral symmetry, i.e. by assuming that the approximate Killing vector has the form where we have added a radial component to the helical vector and assume that orbital angular velocity is along the x 3 -direction. Here v r is the radial velocity, x i CM the center of mass position and r 12 the distance between the two black holes. We can now adjust Ω and v r to mimic true inspiral orbits. In comoving coordinates we still have so that at spatial infinity the shift now must have the form Within the CTS equations we can now use Eq. (146) as a boundary condition on the shift at r → ∞. Using an iterative method [51,86] it is possible to tune the values of Ω and v r to achieve very low eccentricities on the order of 5 × 10 −5 [86] for conformally flat CTS initial data. Similar results can also obtained using superimposed Kerr-Schild CTS initial data [63]. In fact, this kind of eccentricity reduction works better for CTS based initial data than for CTT puncture initial data. The CTS decomposition gives us a direct handle on the shift (via Eq. (146)) and allows us to directly set certain time derivatives to zero. It produces a preferred lapse and shift that when used in the evolution, result in coordinates that are well adapted to quasi-equilibrium. For CTT based initial data one does not have such a preferred lapse and shift. Hence oscillations in the coordinate distance r are not due to real eccentricity alone, but also due to the fact that the coordinates are still evolving as well. This problem is quite visible after eccentricity reduction in Fig. 1. The broken line shows a dip during the first 200M of the evolution. This non-monotonic behavior makes it harder to define or measure eccentricity, so that the reduction procedure works less well.

Toward more realistic binary black hole initial data
As discussed above, two orbiting black holes will emit gravitational radiation. From post-Newtonian calculations we expect to see a so called chirp signal during the inspiral phase, i.e. a gravitational wave with slowly increasing amplitude and frequency (see e.g. broken line in Fig. 2). The initial data should reflect this and contain such gravitational waves. However, when we construct conformally flat initial data using either the CTT or CTS approach, the resulting conformal factor will be monotonically decreasing as we move away from the black holes. Hence the 3-metric does not contain any wiggles. The same is true for the extrinsic curvature. This means that the initial slice of our spacetime contains no gravitational radiation at all, even though we have two orbiting black holes. This is clearly unrealistic, even though the data satisfy the constraints by construction. When such data are evolved in time the system starts to radiate and over time fills the spatial slice with gravitational waves. The results of a evolution starting from conformally flat puncture initial data is shown in Fig. 2. The solid line depicts the dominant mode of the Weyl scalar Ψ 4 which encodes the emitted gravitational waves. This mode has been extracted at a distance of about 90M from the center of mass. As we can see, there is no gravitational wave signal until disturbances from the black holes reach the extraction radius at a time of about 90M . At this time we see a strong transient high frequency signal. This signal is usually referred to as junk radiation, since it is an artifact of the unrealistic initial data. The broken line shows the same mode computed from a corresponding post-Newtonian calculation, which contains no such junk radiation. Similar problems also occur for superimposed Kerr-Schild type data, because each individual black hole in the superposition is a stationary black hole solution of Einstein's equations that does not contain any gravitational waves.

Post-Newtonian based initial data
We know that post-Newtonian calculations can be highly accurate during the inspiral phase before the black holes merge. This is also the regime in which we would like to construct initial data. It thus seems natural that we should try to incorporate post-Newtonian information into the initial data construction. Notice, however, that post-Newtonian theory is an approximation that assumes low velocities (v/c 1) and weak gravity (GM/r 1). While black holes may be moving slowly enough when they are still well separated, their gravitational fields are always strong near each black hole. For this reason post-Newtonian theory alone can give us reliable initial data only away from each black hole. In [87,88] post-Newtonian initial data in ADMTT gauge [89] has been investigated. In this gauge one can directly obtain the 3-metric and the extrinsic curvature as post-Newtonian expansions. It is then possible to resum these expansion so that they take the form and with Here E A are the energies of the point particles used in the post-Newtonian theory, h T T ij and k ij T T are higher order post-Newtonian terms that we do not write out here, and A ij BY is the Bowen-York extrinsic curvature given in Eq. (125) with S i A = 0. As we can see the conformal factor is of Brill-Lindquist form as in Eq. (119). Thus these approximate initial data look very similar to puncture initial data. The main difference is the appearance of the extra terms h T T ij and k ij T T , so that the 3-metric is no longer conformally flat. It is thus clear that these data do indeed contain black holes, even though the post-Newtonian expressions were derived using point particles. Following the approach in [87,88] we can now construct constraint satisfying initial data by using the CTT decomposition with the free datā and using the puncture ansatz for the conformal factor. Inserting these expressions in Eqs. (42) and (44) we obtain elliptic equations for u and W i . These equations can be solved if we impose the boundary conditions lim r→∞ u = 0 (154) and Eq. (61). We then obtain initial data that satisfy Hamiltonian and momentum constraints. This was first done in [87] using a near zone approximation for the term h T T ij in Eq. (147), but the full h T T ij was later calculated in [88], and used to solve the constraints in [90]. As we can see from Eq. (153) the conformal factor gets modified when we solve for the constraint Eqs. (42) and (44). This leads in general to an increase in the ADM mass. Such changes are common when we use one of the conformal decompositions to find constraint satisfying data from approximate data. In [87] and [90] it was found that for the case of the post-Newtonian data in ADMTT gauge, this change in the mass can be prevented by also modifying the post-Newtonian conformal factor to ψ P N → ψ P N − q m 1 m 2 2r 12 This modification (parameterized by the number q) counters the effect of the term u in Eq. (153) on the ADM mass. The parameter q should be chosen such that the resulting initial data after solving the constraints in Eqs. (42) and (44)  such that the binding energy after solving the constraints is as close as possible to the post-Newtonian binding energy. In Fig. (3) we show the dominant l = m = 2 mode of the gravitational waveform for different initial data. As we can see, the ADMTT based initial data leads to better waveforms with reduced junk radiation. Similar improvements in the waveform even exist if we do not solve the constraint equations and just directly use the approximate post-Newtonian data [91,78].

Approximate initial data from matching
As we have seen, post-Newtonian theory can give very useful input when we want to construct initial data with gravitational waves already built into the initial slice. Yet, as discussed in the previous subsection, post-Newtonian theory is not valid near a black hole. In fact one had to judiciously resum the post-Newtonian expressions to even obtain black holes. These resummations were guided by the aim to make the post-Newtonian data similar to puncture initial data, because the latter contains well understood black holes. Of course this resummation is ad hoc, and it is not clear if we really obtain a good approximation near the black holes and whether the black holes are in equilibrium. In fact, since there is still some amount of junk radiation, we do know that they are not really equilibrium. For this reason an alternative method based on asymptotic matching has been proposed [92,93,94,95]. In this method one matches a post-Newtonian 4metric to two different perturbed black hole metrics (one for each black hole) to obtain an approximate solution everywhere. Then no ad hoc resummations are necessary, and we obtain an approximation that is valid everywhere. The most sophisticated data of this kind to date has been presented in [95], where a post-Newtonian 4-metric for circular orbits containing outgoing gravitational waves has been matched to two tidally perturbed Schwarzschild metrics. When we evolve initial data coming from this approximate 4-metric without solving the constraint equations, we find that the data have waves built in from the start and that there is less junk radiation than for puncture initial data [78]. However, somewhat surprisingly the junk radiation in the dominant l = m = 2 mode has still about the same size as the one coming from the resummed ADMTT gauge post-Newtonian data discussed in the previous subsection, at least for a black hole separation of 10M . However, further studies are necessary as it is likely that the situation will improve once we move the black holes further apart, since the then all the approximations used in the matching will improve.
In [96] the two tidally deformed black hole 4-metrics from [95] (but not the post-Newtonian 4-metric) have been used to construct a superposition as in Eqs. (136) and (137). This superposition thus contains black holes with the correct tidal deformations, but without any realistic gravitational waves, since the matched post-Newtonian 4metric is not included. The constraint equations are then solved for this superposition following the methods described in Sec.5.4 within the CTS approach. The only difference is that as boundary condition at each black hole, [96] uses a Dirichlet condition for the shift, since tidally deformed black holes do not have a vanishing shear, so that Eq. (102) may not be appropriate. When these data are evolved, junk radiation is reduced when compared to conformally flat or superimposed Kerr-Schild data.

Matter equations for perfect fluids
So far we have been discussing how one can create initial data for vacuum solutions such as black holes. When we need initial data for situations where matter is involved, we have to consider additional equations that describe this matter. Here we will concentrate on matter that can be described as a perfect fluid. This is sufficient for e.g. neutron stars.

Perfect fluids
For a perfect fluid the stress-energy tensor is given by Here ρ 0 is the mass density (which is proportional the number density of baryons), P is the pressure, is the internal energy density divided by ρ 0 and u µ is the 4-velocity of the fluid. The matter variables in Eq.(32) are then where u 0 = −n µ u µ /α. From ∇ ν T µν = 0 we obtain the relativistic Euler equation which together with the continuity equation governs the fluid. The five fluid Eqs. (158) and (159) are not enough to determine the the six quantities ρ 0 , , P and u i . We thus also need an equation of state of the form to close the system of fluid equations. In many cases it is also useful to introduce the specific enthalpy because then the Euler equation can be simplified to Using u µ u µ = −1, the latter can also be written as u · dp = 0.
Here u is the fluid velocity 4-vector and the one-form p has the components The dot in Eq. (163) indicates contraction of the vector u with the first index of the two-form dp .

Expansion, shear and rotation of a fluid
Using the projector covariant derivatives of the fluid 4-velocity we can be split into where the expansion, shear, rotation and acceleration of the fluid are defined as and If the 4-velocity is of the form where f is any scalar function, it immediately follows from P µν u ν = 0 that From the latter is is clear that if the velocity is derived from a potential φ, i.e. if hu = p = dφ, we immediately obtain ω µν = 0, which characterizes an irrotational fluid.

Non-spinning neutron stars
Let us consider a static spherically symmetric star. In this case the Einstein equations can be solved exactly. The metric outside the star is given by the Schwarzschild metric of Eq. (103), but with m replaced by m(r s * ), where r s * is the star radius (in standard Schwarzschild coordinates) and is the gravitational mass inside radius r s and Inside the star the metric is given by where m(r s ), Φ(r s ) and also P (r s ) are found by integrating the Tolman-Oppenheimer-Volkoff (TOV) equations [97,98] dm together with some equation of state that allows us to obtain ρ E from P . The TOV equations are ordinary differential equations that have to be integrated out from r s = 0 to the radius where P = 0, which is the location of the star surface (r s = r s * ). At r s = 0 we start with m = 0 and some value of P = P c which determines the core pressure and thus the total mass of the star. We also have to start with some particular Φ at r S = 0. We can set this value to 1 at first. The final Φ(r) can be obtained by adding a constant to it such that 2Φ(r s * ) = ln(1 − 2m(r s * )/r s * ). This shift in Φ ensures that the metric is continuous at the star surface. It is possible, and often convenient, to transform the TOV metric to isotropic coordinates by using the transformation given in Eqs. (104) and (105). From the TOV metric it is then easy to obtain the 3-metric and extrinsic curvature.

Spinning neutron stars
The metric for a spinning neutron star is no longer spherically symmetric. Instead one assumes a stationary and axisymmetric metric of the form where A, B, C and ω are functions of r and θ only. If matter is again treated as a perfect fluid and an equation of state is provided, one obtains a system of three field equations and one equation expressible as a line integral [99], both inside the star and in the vacuum region on the outside. This system of equations is rather complicated and we will not go through its derivation here. It can only be solved numerically. Several groups have worked on this problem in the past. A freely available public code (called RNS) that can solve the system of equations has been developed by Stergioulas and Friedman [100]. However, the most accurate results by far have been achieved by Ansorg et al. [99,101]. It is also worth mentioning that the metric outside a spinning star is not given by the Kerr metric that is valid for a rotating black hole. A review about a single spinning neutron stars can be found in [102].

Binary neutron star initial data
Binary neutron stars, like binary black holes, will be on approximately circular orbits if they have been inspiraling already for a long time. We can thus use many of the same methods as for black holes when we construct binary neutron star initial data. Thus again we can assume that an approximate helical or inspiral Killing vector ξ µ as in Eq. (144) should exist. So again we can choose coordinates where Eq. (80) holds, and in which time derivatives of metric variables can be approximated by zero. As mentioned previously, such a situation can be captured best by using the CTS decomposition, which we will use almost exclusively for the binary neutron star initial data described below. We will also restrict ourselves to the case of a conformally flat metric as in Eq. (117). Since gravity in neutron stars is weaker and rotational velocities are lower than in black holes, conformal flatness is generally a good approximation for neutron stars. Notice however, that a very interesting non-conformally flat method, named the waveless formulation, has been developed [103,104,105,106]. In this formulation one obtains the same fluid equations as described below. The equations for the metric variables are essentially the CTS equations, plus one additional elliptic equation for the conformal metric. This extra equation is derived from the evolution equation for the extrinsic curvature by using a particular gauge to rewrite R ij in Eq. (29) as an elliptic operator acting on the conformal metric. Thus the conformal metric does not need to be assumed to be of a particular form, rather it is computed from the additional elliptic equation. Instead one needs to specify as free data the time derivative of a conformally rescaled tracefree extrinsic curvature, in addition to the time derivative of the conformal metric.
The main difference between binary neutron star initial data and binary black hole initial data is that we now also have to find quasi-equilibrium configurations for the matter, when we solve the Euler and continuity equations in Eqs. (158) and (159). Notice, however, that the stars have only a finite extent so that we will use the same boundary conditions for the metric variables as in Sec. 3.3 at infinity.
As in the black hole case, we will not discuss every possible method to construct binary neutron star initial data here. Rather we will concentrate on conceptual issues, as well as on describing the most widely used methods.

Corotating binary neutron stars
In order to construct initial data for two neutron stars in a corotating configuration, we assume again a helical Killing vector ξ µ . The fact that the two stars are corotating means that the orbital period and both star spin periods all have the same value. Thus each fluid element is at rest in corotating coordinates. This idea is captured by assuming that the fluid in each star flows along the Killing vector ξ µ , i.e. that the fluid 4-velocity is In this case it is easy to show that the continuity equation (159) is identically satisfied. Furthermore one can show that the Euler equation leads to (see e.g. problem 16.17 in [107]) where ξ µ is the assumed helical Killing vector. With the help of the first law of thermodynamics (d[ρ 0 (1 + )] = [ρ 0 (1 + ) + P ]dρ 0 /ρ 0 ) this equation can be integrated to yield where C is a constant of integration that is different for each star. Using Eqs. (182) and (161) we arrive at which tells us the value of the specific enthalpy for each point inside the star. If we assume a polytropic equation of state we can express the mass density, the pressure and the internal energy in terms of h. We obtain The constant n here is known as the polytropic index, and κ is a constant. In order to construct binary neutron star initial data we have to solve the five elliptic equations in Eqs. (42), (52) and (56) with the matter terms given by Eqs. (157), (185) and (187). This is done through an iterative procedure [108] as described in subsection 9.4.
Examples of such corotating initial data can be found in [109,110,111,112,108]. The problem with such data is that it is very unlikely for two neutron stars to have spin periods that are both synchronized with the orbital period. As pointed out by Bildsten and Cutler [113], the two neutron stars cannot be tidally locked, because the viscosity of neutron star matter is too low. Hence barring other effects, like magnetic dipole radiation, the spin of each star remains approximately constant. This means that initial data sequences of corotating configurations for different separations cannot be used to approximate the inspiral of two neutron stars.

Irrotational binary neutron stars
Since corotating configurations are unrealistic, many groups use initial data for irrotational stars. The advantage is that, sequences of irrotational configurations can be used to approximate the inspiral of two neutron stars without spin. Of course, for such configurations the fluid velocity will not be along the helical Killing vector, which we will still assume to exist approximately. Thus the fluid velocity is now written as where u 0 = −u µ n µ /α, and V µ is a purely spatial vector with V µ n µ = 0 that describes the fluid velocity relative to the Killing vector ξ µ . In order to find the enthalpy for a quasi-equilibrium situation we insert Eq. (188) into Euler Eq. (163) to obtain ξ · dp + V · dp = 0. (189) Note that for any vector v and form ω the Cartan identity holds, where £ v ω is the Lie derivative of ω along v. The dot in Eq. (190) indicates contraction of a vector with the first index of a form. Using the Cartan identity to replace ξ · dp , Eq. (189) simplifies to The first term in this equation term will be dropped since ξ is an approximate Killing vector, i.e.
where we again have used corotating coordinates as in Eq. (80). From Eq. (163) we can see that setting will automatically solve the Euler equation. The function φ is a velocity potential from which we can then derive p µ = hu µ . Notice that a fluid with a velocity coming from a potential, according to hu µ = ∇ µ φ, has zero rotation (see Eqs. (171) and (174)), hence the name irrotational configuration for the initial data constructed in this section. If we insert Eq. (193) together with assumption (192) into Eq. (191) we finally obtain which implies that where C is a constant that may be different in each star. Inserting ξ µ = u µ /u 0 − V µ from Eq. (188) and using u µ u µ = −1 we arrive at As in the case of corotating configurations this last equation can be used to find the specific enthalpy. Note that u 0 can be obtained from u µ u µ = −1.
Using Eqs. (188) and (193) one can show that the continuity equation (159) becomes [114,115,116] where we have again dropped some Lie derivatives with respect to the Killing vector ξ.
More details about this derivation will be provided in the next section, where we consider the more general case of spinning stars. Compared to the corotating configuration, the irrotational case yields the additional elliptic Eq. (197), that we need to solve for the velocity potential φ together we the other elliptic equations for the metric variables in Eqs. (42), (52), (56). Such irrotational binary neutron stars have been produced by many groups [117,118,119,120,121,122,123,124,125,126].

Binary neutron stars with arbitrary spins
We will now consider the general case of neutron stars with spin. First we assume that the binary is in an approximately circular orbit and that the spin of each star remains approximately constant. As in the case of binary black holes (see e.g. [31,37]) this implies the existence of an approximate helical Killing vector ξ µ with £ ξ g µν ≈ 0. In order to clarify the meaning of the approximate sign we now briefly discuss two cases. If both star spins are parallel to the orbital angular momentum we have £ ξ g µν = O(P o /T ins ), where the inspiral timescale T ins is much longer than the orbital timescale P o . I.e. in a corotating coordinate system all metric time derivatives are of order O(P o /T ins ) and thus small. For arbitrary spins the situation becomes more complicated. We can again use corotating coordinates, but in this coordinate system the spin vectors will be precessing on an orbital timescale P o . This means there are matter currents that change on a timescale P o , and thus £ ξ u µ cannot be zero [116] even approximately! On the other hand, the matter distribution itself only changes on the much longer inspiral timescale T ins . In this case it is useful to consider gravity to be made up of gravitoelectric and gravitomagnetic fields [127,15]. The gravitoelectric parts of the metric are sourced by the matter distribution and thus change only on the timescale T ins , while the gravitomagnetic parts of the metric are sourced by matter currents and thus change on the shorter timescale P o . However, the gravitomagnetic parts are smaller than the gravitoelectric parts by O(v/c) [127,15]. Thus we now have £ ξ g µν = O(v/c) ≈ 0, where we assume that the orbital velocity v is smaller than the speed of light.
We start again by writing the fluid velocity as in Eq. (188) and use again the Cartan identity to bring the Euler equation into the form in Eq. (191). However, this time we cannot assume that £ ξ p is zero. Also, since both stars are spinning, the fluid flow is no longer purely irrotational, so that now we write where dφ and w denote the irrotational and rotational pieces of the velocity. Here w is chosen to be purely spatial, i.e. such that n · w = 0. It is clear that dφ describes the orbital or center of mass motion of the star, while w is due to the star spin. Using Eq. (198) together with the Cartan identity (190), the last term in Eq. (191) becomes so that Eq. (191) is equivalent to Notice that according to Eqs. (188) and (198) so that Eq. (200) finally becomes Following [116], we will now argue that all three terms on the right hand side of Eq. (202) are small. Since the center of mass motion described by dφ should be unchanged along the approximate symmetry vector ξ, we will assume that As we have explained, the spin and thus w changes along ξ. However, if we go along the center of mass direction given by ∇φ, we expect the spin to be unchanged. Thus we assume Furthermore we assume that which tells us that the spinning motion of the star does not change if we go along the direction of this spinning motion. The foregoing argument is also illustrated in Fig. 4, which schematically depicts the star over a quarter orbit. As the star moves, its spin stays approximately constant and its center moves along the Killing vector ξ. If we start at point P and follow ξ we end up at point P . If instead we follow the center of mass motion parallel to ∇φ we arrive at Q. As we can see the orbital fluid speed (excluding the spin) at P and P are the same, which is the origin of assumption 203. However, the fluid motion due to the spin at P and P are not the same. Rather we expect the spin motion to be approximately the same at point P and Q, which supports assumption 204.
Using the three assumptions in Eqs. (203), (204) and (205), we can integrate Eq. (202) and obtain Schematic picture of a star and its spin S: The star starts on the bottom and then moves along a circular arc. Its center follows the Killing vector ξ. The star spin S is approximately unchanged over an orbital timescale. The point P becomes point P if we follow the orbit of the Killing vector, while it becomes point Q if we follow the center of mass motion parallel to ∇φ. The orbital fluid speed (excluding the spin) at P and P are the same. Yet, the fluid motion due to the spin at P is not the same as at P . Rather it should be approximately the same as at P and Q.
where C is the constant of integration. This last equation is again used to find h. It reduces to Eq. (185) in the corotating case where V = 0, and has exactly the same form as Eq. (196) for the irrotational case. In order to to find u 0 we use We now multiply by h 2 , use Eq. (198) and solve for u 0 . We find If we insert Eq. (208) into the integrated Euler equation (206) we can solve for h. This yields where we use the abbreviations We now rewrite the continuity Eq. (159) as Using Eq. (27) and collecting terms we obtain Since ξ is an approximate Killing vector we can drop the last term. We also assume that If we project this equation onto the spatial slice using γ µ ν we find When the latter is inserted into Eq. (213) we obtain Notice that this expression reduces to Eq. (197) in the irrotational case where w i = 0. Equation (217) is only valid inside the star and we have to provide a boundary condition at the star surface where ρ 0 = 0. This condition can be obtained directly from Eq. (217) in the limit of ρ 0 → 0 but D i ρ 0 = 0. We find [116] Notice that the rotational velocity piece w i can be freely chosen. In [128] it is demonstrated that is a good choice, that results in almost rigidly rotating fluid configurations with low expansion and shear. Here x k C * denotes the location of the star center and ω j is an arbitrarily chosen vector that determines the star spin. Summation over the repeated indices j and k is implied To close the system of equations we need again an equation of state. The simplest choice is the polytrope given in Eq. (186), which again results in Eqs. (187). However, other choices are possible. For example, [129] uses piecewise polytropic equations of state to approximate several realistic equations of state at zero temperature. Figure 5 shows velocity vectors of fluid elements in the orbital plane for two binaries. In the irrotational case (on the left) the fluid velocity is in the orbital direction. For the case where both stars are spinning (right plot) the fluid velocity has an additional rotational contribution due to w i in Eq. (198).
The same set of equations discussed above has also been used in [130,131]. Notice, however, that [130] erroneously assumes £ ξ p = 0 = £ V w . But as shown in the appendix of [116], £ ξ p cannot be approximately zero for generic spins. Fortunately the two errors incurred by the latter two assumptions cancel, so that the right hand side of Eq. (200) is again zero, and all other equations follow as shown above. Spinning binary neutron initial data constructed with the method described above are evolved in [131], and the star spins are monitored using a newly defined quasi-local neutron star spin. It is found that the star spins follow post-Newtonian predictions during the inspiral of the two stars. This indicates that the above formalism together with the choice for w i in Eq. (219) works and results in realistic star spins.

Iteration procedures required to construct binary neutron star initial data
As we have discussed above, the construction of binary neutron star initial data requires us to solve the elliptic equations of the CTS approach in Eqs. (42), (52) and (56)  Eq. (185) are algebraic equations that specify the specific enthalpy h given the other fields that obey elliptic equations. The situation is complicated by the fact that the star surfaces are at the location where h = 1, so that one cannot know where the star surfaces will be before solving all coupled equations. Thus we do not know in which regions we indeed have to include matter terms in Eqs. (42), (52) and (56). Notice also that the values of the matter terms are undefined if we blindly use the equation of state in regions where h < 1 (see e.g. Eqs. (187)).
For these reasons the Eqs. (42), (52), (56), (217) and (209) are solved using an iterative procedure [128], that can be outlined as follows: (i) We first come up with an initial guess for h in each star (in practice we often choose Tolman-Oppenheimer-Volkoff solutions). (ii) Next, we solve the 6 coupled elliptic equations in Eqs. (42), (52), (56) and (217) for this given h. (iii) Then we use Eq. (209) to update h in each star. After updating h we go back to step (ii) and iterate until all equations are satisfied up to a given tolerance.
There is however one major ingredient missing in the above outline. Our entire method assumes a symmetry vector ξ of the form (144). Thus we need particular values for Ω and x 1 CM . Now, the star surfaces move between iteration steps. If we simply specify constant values for Ω and x 1 CM , it turns out that the stars drift away and that the iteration becomes unstable. We thus have to find a way to keep the stars from drifting. This is usually done by adjusting Ω and x 1 CM such that the star center does not move when we update h using Eq. (209). Notice that Eq. (209) depends on Ω and x 1 CM through ξ i which is contained in b from Eq. (211). The star centers x 1 C * 1/2 are defined by the maxima of h along the x 1 -axis, given by ∂ 1 h| x 1 C * 1/2 = 0. In [128,129] it is shown that the latter can be expressed as We then update Ω and x 1 CM by solving Eq. (220) for Ω and x 1 CM so that the star centers x 1 C * 1/2 remain in the same location, when we update h according to Eq. (209). Equation (220) is often called the force balance equation [118]. One noteworthy caveat is that we evaluate the derivative of ln Γ in Eq. (220) for the Ω and x 1 CM before the update.
Finally, the complete iteration procedure involves the following main steps [128]: The constant C has a different value for each star. We adjust the value for each star such that it keeps a prescribed rest mass. After updating h we go back to step (2) and iterate until all equations are satisfied up to a given tolerance.
Notice however, that sometimes other conditions are used instead of the force balance Eq. (220). As first found in [125], for massive stars in can be helpful to determine x 1 CM by demanding that the ADM linear momentum be zero [132,133,129,126].
We only need to specify x 1 CM because we can always rotate our coordinate system such that the line connecting the two star centers is along the x 1 -axis so that x 2 CM = 0 = x 3 CM .

General orbits and low eccentricity orbits
If we evolve the quasi-circular (i.e. with v r = 0 in Eq. (144)) binary neutron star initial data described in the previous sections, we find that the orbits are roughly circular inspiral orbits. However, as in the case of black holes (see e.g. Fig. 1), the distance between the stars exhibits small oscillations, which are a sign of residual eccentricity caused by our assumption of exactly circular orbits. As for black holes, we would like to adjust Ω and v r in Eq. (144), to obtain less eccentric orbits when we evolve such initial data. There is, however, one immediate problem. The value of Ω is determined by the force balance Eq. (220), and thus cannot be readily changed.
In [134] a method has been developed that lets us construct orbits with arbitrary eccentricity. The extra parameter e describing this eccentricity can then be adjusted along with v r to find low eccentricity orbits when the data are evolved. We now describe how the symmetry vector ξ µ can be modified to account for non-circular orbits. We make the following two assumptions: (i) ξ µ is along the motion of the star center. (ii) Without inspiral, each star center moves along a segment of an elliptic orbit at apoapsis. Since we only need a small segment of an elliptic orbit near apoapsis, we will approximate this segment by the circle inscribed into the orbit there. Then the radii of the inscribed circles are where d 1 and d 2 are the distances of the particles from the center of mass at apoapsis and e is the eccentricity parameter for the elliptic orbit [134]. These two inscribed circles are centered on the points where we have used d 1/2 = |x 1 C * 1/2 − x 1 CM | and assumed that apoapsis occurs on the x 1 -axis. (The minus and plus sign in the middle term of Eq. (223) corresponds to the subscripts 1 and 2, respectively.) Assumptions (i) and (ii) then tell us that the symmetry vector for elliptic orbits with inspiral must have the form The eccentricity parameter e that appears in x 1 c1/2 and the radial velocity v r can now be freely chosen to obtain any orbit we want. Using this new symmetry vector ξ 1/2 , we can still solve the initial data equations with the same methods as described above. In order to obtain a true inspiral orbit with low eccentricity, we have to adjust both e and v r . A procedure (similar to the one for black holes in [51,86]) to minimize the eccentricity has been developed in [129]. A similar approach has been also been applied in [135] and [126] for purely irrotational stars. Notice, however, that in [135] and [126] only quasi-circular orbits are considered and Ω and v r are directly adjusted. Similarly [131] considers quasi-circular orbits using the methods for spinning stars discussed above and also directly adjusts Ω and v r .
As we have seen, it is possible to construct binary neutron star initial data for spinning stars on arbitrary orbits. Nevertheless, due the many technical complications discussed above, other simpler approaches are sometimes still used. For example, one can use simple superpositions of single stars without solving the constraints and the matter equations [136,137], or possibly with only the constraint equations solved [138], or one can modify the velocity field of irrotational stars to introduce spin [139]. Of course such simple methods introduce additional errors that are not well understood.

Initial data for black hole -neutron star binaries
Binaries made up of one black hole and one neutron star can be constructed by combining the methods for black holes and neutron stars discussed above. As we have seen, there are a considerable number of choices on how one can construct both black holes and neutron stars, thus there are many possible combinations of these choices when it comes to mixed binaries.

Simple approaches
One of the simplest approaches is used in [140], where the black hole is assumed to be much more massive than the neutron star. The black hole is described by using the Kerr-Schild metric. The matter equations and constraints are then only solved in a small patch around the neutron star, which is assumed to be irrotational. The rest of the spacetime is approximated to be given by the Schwarzschild metric.
As with neutron star binaries it is possible to use superpositions of single neutron star and single black hole solutions. In [141,142] a boosted TOV star solution is superimposed with a boosted black hole metric to obtain highly eccentric orbits. This superposition solves the constraint equations only approximately, and the fluid is not in true equilibrium. Yet if the one starts from large separations the errors due to a simple superposition are small. This approach has also been extended [138] to using this superposition as free data in the CTS decomposition and then solving the constraint Eqs. (42) and (52), but without solving the equilibrium fluid equations. This method has recently also been generalized to spinning neutron stars by using a solution for a single spinning star in the superposition [143].

Non-spinning black hole -neutron stars binaries
Several other more comprehensive approaches exist for general mass ratios. In [144] the 3-metric is assumed to be conformally flat. Both black hole and neutron star are assumed to be non-spinning. Then the CTS Eqs. (42), (52), (56), together with the matter Eqs. (196) and (197) can be solved to set up an irrotational neutron star. The black hole is introduced by imposing the quasi-equilibrium apparent horizon boundary conditions (99) and (102) on a coordinate sphere. The Ω r in Eq. (102) is adjusted such that the black hole spin is zero. Again an iteration procedure, as described above for neutron star binaries, is used. During the iterations, the center of mass x i CM is fixed by demanding that the ADM momentum P ADM i ∞ = 0, while Ω is computed from the force balance Eq. (220) at the center of the star. Similar methods are also used in [145]. The main difference is that Ω r in Eq. (102) has the same magnitude as Ω so that the black hole has only approximately zero spin. This approximation is later improved in [146]. In [147] the same methods as in [145] are used, except that the conformal metric is set to be the Kerr-Schild metric.
All the approaches that use the quasi-equilibrium apparent horizon boundary conditions result in initial data only outside the black hole. For this reason a puncture based approach has been developed in [148,149]. One assumes a conformally flat metric as in Eq. (117) and sets K = 0. Then Eq. (42) can again be used with the puncture ansatz in Eqs. (128) and (119) with the black hole bare mass m 1 and m 2 = 0, to find the conformal factor ψ. Similarly, we solve Eq. (56) for the lapse with the ansatz in Eq. (142) with m 2 = 0. The extrinsic curvature is written as whereĀ ij BY is the Bowen-York extrinsic curvature given in Eq. (124) containing the black hole momentum (the spin is set to zero). The latter equation corresponds to the CTT decomposition withσ = 1 andM ij =Ā ij BY . We then use Eq. (44) to find W i , where we note thatD jM ij = 0 = K. As noted in [148,149], the fluid equations contain the shift, which we thus need to find. The shift is computed by rewriting the momentum constraint (52) of the CTS decomposition as All terms in this equation are regular at the black hole center (r 1 = 0), as are all terms in the other elliptic equations. We can thus solve all elliptic equations for the metric variables (u from Eq. (128), v from Eq. (142), W i and β i ) everywhere and do not need to excise any regions. This approach has been used together with the fluid Eq. (185) for corotating stars [148,149] and also with Eqs. (196) and (197) for irrotational stars [150]. As in other approaches, Ω is determined from the force balance Eq. (220). The black hole momentum parameter P i inĀ ij BY is picked such that P ADM i ∞ = 0. In order to find the center of mass x i CM , [150] chooses J ADM ∞ such that it corresponds to the post-Newtonian angular momentum.

Spinning black hole -neutron stars binaries
In order to deal with highly spinning black holes it is useful to again consider a conformal metric that is of Kerr-Schild form in the vicinity of the black hole. This has been done [151,132] within the CTS approach. Following [151,132], we again assume that all time derivatives suchū ij and ∂ t K can be set to zero. We then have to specify the conformal 3-metric and the trace of the extrinsic curvature. They are set tō where γ 1 ij and K 1 are the the conformal 3-metric and the trace of the extrinsic curvature of a boosted spinning Kerr-Schild black hole. The boost velocity is given by where c i BH is coordinate location of the black hole center with respect to the center of mass. The weighting factor w is chosen such that the conformal metric is flat at the neutron star. The exponential damping of the form e −r 4 1 /w 4 used here is necessary [151] to avoid the deviations from quasi-equilibrium reported in [145]. When then solve the CTS Eqs. (42), (52) and (56) together with e.g. the fluid Eqs. (196) and (197) for an irrotational neutron star.
We assume that the approximate symmetry vector is given by where we have added the term v i GT corresponding to a Galilean transformation of the helical inspiral vector, in order take into account the fact that the system may move since it radiates linear momentum. Again we assume that orbital angular velocity is along the x 3 -direction. Here v r is the radial velocity, x i CM the center of mass position and r 12 the distance between the black hole and the neutron star. In comoving coordinates we still have ξ µ sym = t µ = αn µ + β µ , so that at spatial infinity the shift now must have the form As boundary conditions for ψ and α we again use Eqs. (59) and (62) at infinity. The apparent horizon S of the black hole is chosen to be the coordinate location of the event horizon of the boosted black hole in the superposition. For the conformal factor and the shift Eqs. (99) and (102) are used as boundary conditions on S. The lapse boundary condition there is given by where α 1 is the lapse of the boosted black hole. This set of equations again has to be solved iteratively. When low eccentricity orbits are desired, Ω and v r are computed by an eccentricity reduction procedure as in [51,86]. We can thus not use the force balance equation (220) to center the star. Rather we fix the positions of both star and black hole centers (say at given locations on the x 1 -axis), and then translate the specific enthalpy h computed in each iteration so that its maximum occurs at the fixed star center. Finally, P ADM i ∞ = 0 is enforced by adjusting v i GT . This is called boost control in [133] and is particularly useful for very compact stars. The details of this procedure can be found in [151,132,133].
With this approach it is possible to construct initial data for neutron stars that are orbiting highly spinning black holes. So far this has only been done for irrotational stars, but the approach could easily be generalized to spinning stars using the methods explained in Sec. 9.3.

Numerical codes used to construct initial data for binaries
Many different computer codes have been used over the years to compute initial data. Here we only mention the ones that are currently most widely used.
One of the most widely used codes is the publicly available Lorene library [152].
Lorene provides a domain decomposition into several spherical-like shells. In each shell, spectral methods in spherical coordinates are employed to solve the relevant elliptic equations. Spectral expansions in terms of spherical harmonics and Chebyshev polynomials are used in the angular coordinates and the radial coordinate, respectively. In order to impose boundary conditions at infinity the outermost shell is compactified by using the inverse of the radius as a coordinate. The shells can be deformed so that a neutron star surface can be at a domain boundary, which allows for easy imposition of surface boundary conditions. For a binary, two sets of domains are used, that are centered around each compact object. These domains are embedded in a larger domain that extends to infinity. The Lorene library has been used in several codes. Among them are the codes that were used to construct binary black hole data [41], binary neutron star data [117,118,122,123,125] and also mixed binaries [144,147,146]. Most numerical simulations starting from binary black hole initial data of puncture type, use the spectral puncture solver described in [153]. This method uses a coordinate transformation that maps all of space into a finite cube. After mapping, the punctures are located along two parallel edges of one of the cube faces, while spatial infinity is mapped to the opposing face. This means that in the interior of the domain ψ will be C ∞ . Equation (129) is then solved by means of a pseudo-spectral method, where expansions in Chebyshev and Fourier polynomials are used. This solver is publicly available as the TwoPunctures thorn of the Einstein Toolkit [154].
Many initial data sets for binary black holes [46,51,63,67] mixed binaries [132,133] and neutron star binaries [131,126] have been computed using the spells framework [155]. In this framework space is covered by computational domains that are made up of hexahedra, spherical shells, and cylindrical shells. Within each domain appropriate spectral expansions are used. For example, Chebyshev polynomials are used in all three directions in hexahedra, spherical harmonics are used in the angular directions of shells and Chebyshev polynomials the radial direction, Fourier polynomials are used in the angular direction of cylindrical shells. The domains can be distorted by coordinate transformations, so that usually grids with domains are used that touch but do not overlap. For neutron stars, the star surface is put at a domain boundary.
The first evolutions of binary neutron stars with realistic spins [156,129] were performed using initial data computed with the SGRID code [108,128,129]. This code uses coordinate transformations [157] that compactify all of the space outside the neutron stars into two cubical domains. In each of these domains the star surface is one of the cube faces and spatial infinity is mapped into an edge of the opposing face. Two more domains are introduced for each star that cover the inside of each star. The outer one of these is again cubical with the star surface located on one of the faces. The coordinate transformations contain adjustable functions so that arbitrary star surfaces are possible. As long as the matter distribution is C ∞ inside each domain, the solutions of elliptic equations such as Eqs. (42), (52), (56), and (217) are expected to be C ∞ as well. In SGRID these equations are thus solved by means of a pseudo-spectral method, which converges exponentially if the solutions are C ∞ . This code has also been used to construct the post-Newtonian based binary black initial data studied in [90].
The codes mentioned up to here use spectral methods which tend to be more efficient than finite difference methods due to their faster rate of convergence when the computed solutions are C ∞ . However, the simulations of binary neutron stars and mixed binaries on highly eccentric orbits in [143] have used a finite differencing code [138] to construct the initial data. This code uses a full approximation storage implementation of the multigrid algorithm with adaptive mesh refinement as described in [158].
Recently the so called COCAL code has been developed [159,160,161,130]. It is different from the codes mentioned before and is based on the Komatsu-Eriguchi-Hachisu method [162] for computing equilibrium spinning neutron star models. In this method elliptic equations are solved by using Green's functions for the Laplace operator. Any additional or non-linear terms in the equations are added to the source terms, and the equations are then solved iteratively. The COCAL code uses several overlapping spherical shells. The Green's functions include terms for different boundary conditions. The integrals over the Green's functions is performed numerically by expanding into spherical harmonics. A comparison of COCAL and Lorene in [130] shows that COCAL gives the same results up to a difference of 0.05%.

Conclusions
The construction of initial data in General Relativity is a well established field of study these days. As we have seen, such initial data have to satisfy both the Hamiltonian and momentum constraints in Eqs. (30) and (31). In order to solve these constraints, conformal decompositions have been developed, that result in elliptic equations which can be solved numerically. The two most common approaches are the CTT and CTS decompositions discussed in Secs. 3.1 and 3.2. In both cases one has to specify certain free data such as the conformal metric and boundary conditions. Once specified one can then construct constraint satisfying initial data. As we have seen both decompositions can be and have been used to construct various types of initial data. The answer to the question of which one is preferable depends on what one tries to achieve. If we have a way to supply approximate initial data (e.g. using post-Newtonian input) that almost satisfy the constraints one can use either the CTT or CTS decomposition to construct constraint satisfying initial data. In this case it does not matter much which decomposition one chooses, except maybe that the CTT equations are slightly simpler. On the other hand, if one tries to construct binary initial data solely from quasiequilibrium conditions (e.g. for two objects on quasi-circular orbits), it is better to use the CTS decomposition because it allows us to directly specify certain time derivatives (which would be zero for quasi-equilibrium). This is why this method is used for most types of binary neutron star initial data, and also for many binary black hole initial data.
The main ideas for constructing binary black hole or binary neutron star data are very similar. However, neutron stars contain matter and cannot simply be modeled by certain boundary conditions on the black hole horizon. This fact adds extra complications that make the actual construction of initial data containing neutron stars more challenging. More equations have to be solved and the numerical iteration procedures become more complex because the star surfaces change during the iterations. Nevertheless, the quality (in terms of astrophysical realism) of current binary initial data containing one or two neutron stars is just as good as the quality of binary black hole initial data. One possible exception is the addition of gravitational waves to the initial data, which thus far has only been carried out for binary black hole initial data. This addition, however, could in principle be done for neutron stars with the same methods as for black holes, e.g. by using post-Newtonian input.

Appendix B. Elliptic Equations
We mention here some properties of elliptic equations that are important for the construction of initial data.
Appendix B.1. Classification of second order partial differential equations Consider a second order partial differential equation for the function u = u(x l ) of the form a ij ∂ i ∂ j u + F (x l , u, ∂ l u) = 0, (B.1) where F is a function of x l , as well as u and its first derivatives, and the matrix a ij is a function of x l and possibly u. When a ij is independent of u and F is linear in u the partial differential equation is called linear. If a ij also depends on u but F is linear in u the equation is sometimes called quasi-linear. If F is not linear in u the equation is called non-linear. The term a ij ∂ i ∂ j u with the highest derivatives is referred to as the principal part of the equation, and determines the type of the equation. If the eigenvalues of a ij are all positive or all negative for all x l , the equation is called elliptic. If one or more eigenvalues are zero, the equation is called parabolic. If no zero eigenvalues exist, but some eigenvalues are positive while others are negative the equation is called hyperbolic.
The definition for elliptic equations can be extended to systems of equations of the form where u C now stands for a vector of functions. The system of equations is called elliptic if [163] a ij AB ξ i ξ j ≥ c|ξ i | 2 , (B .3) for all ξ i , all x l and some constant c > 0.

Appendix B.2. Finding solutions of elliptic equations
Many useful theorems about the existence and uniqueness of solutions, as well as how one can construct solutions, are known for linear elliptic equations (see e.g. [163,164]). Unfortunately, much less is known about strongly coupled non-linear systems of elliptic equations, such as the the CTS system given by Eqs. (42), (52) and (55). For example, in [165] it is demonstrated that there are cases where this system does not admit one unique solution, but instead has several branches of solutions. It is beyond the scope of this article to detail the conditions under which the CTT or CTS systems are known to have unique solutions, and under which conditions not much is known. Note, however, that some of these conditions are listed in [166], and that important results can be found in [167,168,169].
Here we follow a more practical approach and simply note that in most cases of interest to us, we can numerically construct solutions to both the CTT and CTS systems. In order to do so we use the same methods as for linear elliptic equations, i.e. we treat the elliptic equations as boundary value problems. This means that in addition to the elliptic equations, we also have to supply boundary conditions at the boundaries of the domain in which we want to solve the elliptic equations. Only then can we obtain unique solutions. These conditions are usually obtained from physical considerations and are imposed at spatial infinity, at black hole horizons, or at star surfaces. For computational reasons one sometimes further divides the domain of interest into several adjacent computational domains. In this case one also has to specify boundary conditions at these domain boundaries. Since these boundaries are artificial, one simply demands continuity of the functions u C and their normal derivatives across such computational domain boundaries.