Martinez-Troncoso-Zanelli Philippe Kerr

Grandclément

Jordan Nicoules
email: jordan.nicoules@observatoiredeparis.psl.eu

Boundary conditions for stationary black holes ; Application to

This work proposes a set of equations that can be used to numerically compute spacetimes containing a stationary black hole. The formalism is based on the 3+1 decomposition of General Relativity with maximal slicing and spatial harmonic gauge. The presence of the black hole is enforced using the notion of apparent horizon in equilibrium. This setting leads to the main result of this paper: a set of boundary conditions describing the horizon and that must be used when solving the 3+1 equations. Those conditions lead to a choice of coordinates that is regular even on the horizon itself. The whole procedure is validated with three different examples chosen to illustrate the great versatility of the method. First, the single rotating black holes are recovered up to very high values of the Kerr parameter. Second, non-rotating black holes coupled to a real scalar field, in the presence of a negative cosmological constant (the so-called MTZ black holes), are obtained. Last, black holes with complex scalar hairs are computed. Eventually, prospects for future work, in particular in contexts where stationarity is only approximate, are discussed.

I.

INTRODUCTION

Black holes are objects so compact that nothing, light included, can escape their intense gravitational field. If this concept was first considered in the eighteenth century, it is with the advent of General Relativity that the mathematical description of those objects was made possible. The first solution of a single non-rotating black hole was obtained by Schwarzschild in 1916. In 1963, Kerr extended the solution to include rotation and the metric which took his name [START_REF] Kerr | Gravitational field of a spinning mass as an example of algebraically special metrics[END_REF]. Not only does this metric describe a rotating black hole, but it has been proven that, under some assumptions, it is the only possible choice. This result arises from the uniqueness theorems (see [START_REF] Chrusciel | Stationary black holes: uniqueness and beyond[END_REF] for a review).

If it is long known that massive stars should end up their life producing a black hole, direct proofs of the existence of astrophysical black holes are now available with more and more confidence. First there is the detection of the gravitational waves emitted by the coalescence of two black holes. When two such objects are orbiting each other, they lose energy by deforming the spacetime and eventually merge into a single black hole. The emitted gravitational waves can be detected by laser interferometry on Earth [3,4]. Since the first detection in 2015 [START_REF] Abbott | Observation of gravitational waves from a binary black hole merger[END_REF], several tens of such binaries have been detected. The observed waveforms are in total agreement with the prediction of General Relativity [START_REF] Abbott | Tests of general relativity with the binary black hole signals from the LIGO-Virgo catalog GWTC-1[END_REF][START_REF] Abbott | Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog[END_REF].

On the other hand, observations with very high angular resolutions enabled to confirm the existence of supermassive black holes at the center of galaxies, mainly in two cases. First, at the center of our galaxy where the orbits of stars indicate the existence of a small dark object of mass ≈ 4 • 10 6 M . Observations are consistent with this object being a classical black hole described by General Relativity [START_REF] Abuter | Detection of the Schwarzschild precession in the orbit of the star S2 near the Galactic centre massive black hole[END_REF]. Using very-long baseline interferometry it was also possible to directly image an accretion disk around the center of the galaxy M87 [START_REF] Akiyama | Event Horizon Telescope), First M87 Event Horizon Telescope results. I. The shadow of the supermassive black hole[END_REF]. The obtained image exhibits a shadow consistent with the fact that the disk is indeed orbiting a black hole.

It is expected that all those types of detectors will reach higher and higher sensitivity in the coming years. New detectors, like the LISA space interferometer [START_REF] Lisamission | Lisa consortium website[END_REF][START_REF] Amaro-Seoane | Astrophysics with the Laser Interferometer Space Antenna[END_REF] or the Einstein Telescope [START_REF] Maggiore | Science case for the Einstein telescope[END_REF], will also come online . This will enable precise tests of the black hole paradigm. If so far observations are consistent with the compact objects being classical Kerr black holes, there is possibility that they are indeed more complicated ones. An example of such an alternative model is the black hole with scalar hairs studied in Sec. V. To maximize the scientific impact of the future observations, theoretical studies of the various models of black holes are needed. A lot of those studies rely on analytic choices, in particular concerning the coordinates used. This can be a limitation when those coordinates are singular or when fast rotation is included (see [START_REF] Van Aelst | Hairy rotating black holes in cubic Galileon theory[END_REF] for instance).

In this paper, a framework is proposed that enables the numerical description of stationary black holes in a rather general context. The formalism should lead to a choice of coordinates that is regular everywhere, in particular across the black hole horizon. The formalism relies on the 3+1 decomposition of spacetime. The choice of coordinates is based on maximal slicing for the time coordinate and on the spatial harmonic gauge for the spatial ones. The presence of the black hole is enforced by demanding that spacetime contains an apparent horizon in equilibrium. Let us mention that this is far from being the first proposal for such a description. A method based on the spatial symmetries of the spacetimes can be found in [START_REF] Headrick | A New approach to static numerical relativity, and its application to Kaluza-Klein black holes[END_REF][START_REF] Adam | A numerical approach to finding general stationary vacuum black holes[END_REF][START_REF] Wiseman | Numerical construction of static and stationary black holes[END_REF]. It has been applied successfully to various situations (asymptotically anti-de Sitter, various dimensions). A method relying on the use of the Dirac gauge can be found in [START_REF] Vasset | Excised black hole spacetimes: Quasi-local horizon formalism applied to the Kerr example[END_REF] where it is applied to the Kerr spacetime.

The paper is organized as follows. In Sec. II, the formalism is presented. The bulk equations coming from the 3+1 setting are exhibited. The description of the apparent horizon is also investigated in detail and the resulting boundary conditions for the metric fields are given. Three different applications are then shown. In Sec. III, the Kerr black hole is recovered in this coordinate system (which is not analytic). In Sec. IV, a model of black hole studied by Martínez and collaborators is recovered numerically. The model has no angular dependence and contains a real scalar field minimally coupled to gravity. There is also a negative cosmological constant which causes the spacetime to be asymptotically anti-de Sitter. Section V shows the construction of black holes with complex scalar hairs as already obtained in [START_REF] Herdeiro | Kerr black holes with scalar hair[END_REF][START_REF] Herdeiro | Construction and physical properties of Kerr black holes with scalar hair[END_REF]. Future prospects are considered in Sec. VI.

Throughout this paper Greek indices are four-dimensional ones, ranging from 0 to 3 whereas Latin indices are spatial ones, ranging from 1 to 3. Units such that G = c = 1 are used.

II. FORMALISM

A. 3+1 formalism and gauge conditions

The 3+1 decomposition of Einstein's equations is widely used in the field of numerical relativity and it is at the core of this work as well. Basic features of this formalism are recalled (see for instance [START_REF] Gourgoulhon | 3+1 Formalism in General Relativity[END_REF] for more details). The four-dimensional metric g is decomposed as

g µν dx µ dx ν = (-N 2 + B i B i )dt 2 + 2B i dx i dt + γ ij dx i dx j . (1)
The hypersurfaces Σ t of constant time t are mapped by the purely spatial coordinates

x i . The 3+1 quantities are then a scalar function N the lapse, a vector field B i the shift and the metric induced on Σ t , γ ij . All indices of spatial quantities are manipulated by the induced metric. In the following ∇ denotes the covariant derivative associated with g µν and D the one associated with γ ij .

The normal to each slice Σ t is n µ = (-N, 0, 0, 0). In this framework, the second fundamental form, the extrinsic curvature tensor, reads as follows:

K ij = 1 2N (D i B j + D j B i -∂ t γ ij ) . (2)
Each index of Einstein's equations can then be projected either on the hypersurface Σ t or along the normal. It leads to the the 3+1 equations of general relativity:

H : R + K 2 -K ij K ij -2Λ = 16πE (3)
M i : D j K ij -D i K = 8πP i (4)
E ij : -∂ t K ij + L B K ij -D i D j N + N R ij + KK ij -2K ik K k j -Λγ ij (5) = 4πN 2S ij -γ kl S kl -E γ ij ,
where L denotes the Lie derivative, R ij and R the Ricci tensor and scalar and K the trace of the extrinsic curvature tensor. H denotes the Hamiltonian constraint, M i the momentum constraint and E ij the evolution equation. The equations are written here with a cosmological constant Λ and in the presence of matter. Matter terms contain the 3+1 projections of the stress-energy tensor E, P i and S ij .

As such, the system of equations (3-5) (supplemented with Eq. ( 2)), cannot be solved to find the fields N , B i and γ ij . Indeed the choice of coordinates (t, x i ) has not been prescribed yet and the general covariance of the theory would lead to an ill-posed problem.

In the following, one uses the same choice of coordinates that was successfully employed in [START_REF] Martinon | Gravitational geons in asymptotically anti-de Sitter spacetimes[END_REF] and which mathematical properties were assessed in [START_REF] Andersson | Elliptic hyperbolic systems and the Einstein equations[END_REF]. The slicing of spacetime (i.e.

the choice of Σ t ) is defined by the maximal slicing condition K = 0. One also demands that the spatial coordinates fulfill the spatial harmonic gauge (the 3D version of the well-known 4D harmonic one). It amounts to enforcing that

V k ≡ γ ij Γ k ij -Γk ij = 0. (6)
Γ k ij are the Christoffel's symbols of γ ij whereas the Γk ij correspond to a fixed background metric γij . Different choices of the background metric would lead to different choices of spatial coordinates. In standard cases the background metric is usually the flat metric one (in spherical or Cartesian coordinates). When a negative cosmological constant is present (as in [START_REF] Martinon | Gravitational geons in asymptotically anti-de Sitter spacetimes[END_REF] or in Sec. IV) it is convenient to consider spatial metrics linked to the anti-de Sitter (ADS) spacetime.

In [START_REF] Andersson | Elliptic hyperbolic systems and the Einstein equations[END_REF] it has been shown that the Ricci tensor can be expressed as

R ij = - 1 2 γ kl Dk Dl γ ij + 1 2 (D i V j + D j V i ) + terms(∂γ ∂γ), (7)
where D denotes the covariant derivative associated with the background metric. Using

V i = 0 thus ensures that the second order derivatives of the metric appearing in the Ricci tensor are all accounted for by γ kl Dk Dl γ ij , which is a Laplacian-like operator and thus well-behaved (see Sec. (3.2) of [START_REF] Martinon | Gravitational geons in asymptotically anti-de Sitter spacetimes[END_REF] for more details).

In order to enforce the gauge conditions, one removes all the occurrence of K in the 3+1 equations. R ij is also replaced by R ij -1/2 (D i V j + D j V i ), making the associated second order part Laplacian-like. In the following, only stationary black holes are considered so that one can also remove all the terms ∂ t . This leads to the following system of equations:

H : R -D k V k -K ij K ij -2Λ = 16πE (8)
M i : D j K ij = 8πP i (9)
E ij : L B K ij -D i D j N + N R ij - 1 2 (D i V j + D j V i ) -2K ik K k j -Λγ ij (10) = 4πN 2S ij -γ kl S kl -E γ ij ,
where V i is given by Eq. ( 6) and K ij by Eq. ( 2), which in that case reduces to:

K ij = 1 2N (D i B j + D j B i ) . (11)
Eqs. [START_REF] Abuter | Detection of the Schwarzschild precession in the orbit of the star S2 near the Galactic centre massive black hole[END_REF][START_REF] Akiyama | Event Horizon Telescope), First M87 Event Horizon Telescope results. I. The shadow of the supermassive black hole[END_REF][START_REF] Lisamission | Lisa consortium website[END_REF] are now an invertible system of 10 components, corresponding to the 10 unknown fields N , B i and γ ij . Once this system is solved, a very important check consists in verifying, a posteriori, that the gauge fields K and V i are indeed zero. If this is not the case, the solved equations do not coincide with the original ones, making the overall procedure fail. This check is of uttermost importance and it is very difficult to pass if anything is wrong somewhere in the equations. It is carefully monitored for the three different examples presented in this paper. The above procedure is the three-dimensional equivalent of the so-called De Turck procedure used to enforce the four-dimensional harmonic gauge [START_REF] Headrick | A New approach to static numerical relativity, and its application to Kaluza-Klein black holes[END_REF][START_REF] Deturck | Deforming metrics in the direction of their Ricci tensors[END_REF].

Validity of the procedure is discussed, for instance, in [START_REF] Figueras | Ricci solitons, Ricci flow, and strongly coupled CFT in the Schwarzschild Unruh or Boulware vacua[END_REF][START_REF] Figueras | On the existence of stationary Ricci solitons[END_REF].

B. Apparent horizon boundary conditions

It is well-known that event horizons are global objects that require the knowledge of the full spacetime in order to be located. It makes them difficult to use in the context of numerical relativity. In order to enforce the presence of black holes, one usually relies on the local notion of apparent horizon, first introduced in [START_REF] Hawking | Cambridge Monographs on Mathematical Physics[END_REF]. By local one means that apparent horizon can be defined on each slice Σ t by the sole knowledge of the geometry of the slice.

Apparent horizons are commonly used in black hole simulations (see for instance [START_REF] Gourgoulhon | Binary black holes in circular orbits. 1. A Global space-time approach[END_REF][START_REF] Grandclement | Binary black holes in circular orbits. 2. Numerical methods and first results[END_REF][START_REF] Caudill | Circular orbits and spin in black-hole initial data[END_REF] for applications in the binary context). In the following only relevant properties of apparent horizons are discussed (see [START_REF] Gourgoulhon | A 3+1 perspective on null hypersurfaces and isolated horizons[END_REF] for a detailed review). Let us also point out that, in the stationary cases, apparent and event horizon coincide.

Without loss of generality, the apparent horizon is assumed to be a sphere of constant radius r H . Spherical coordinates (r, θ, ϕ) are used to map the slices Σ t and tensors are given with respect to the associated (orthonormal) spherical tensorial basis. Let us denote si the unit, outward, spacelike, normal to the horizon. If need be, it can be shown that si = √ γ rr , 0, 0

= 1 √ γ rr γ rr , γ rθ , γ rϕ . ( (12) si
) 13
The very definition of an apparent horizon means that it is the outermost trapped surface.

So the expansion Θ of the outward future null vector must vanish. Using the 3+1 expression for Θ (see [START_REF] Gourgoulhon | A 3+1 perspective on null hypersurfaces and isolated horizons[END_REF]) it gives rise to the following equation

Θ ≡ D i si + si sj K ij = 0. ( 14
)
Moreover, following [START_REF] Caudill | Circular orbits and spin in black-hole initial data[END_REF][START_REF] Gourgoulhon | A 3+1 perspective on null hypersurfaces and isolated horizons[END_REF], one can ask that the coordinate system is stationary with respect to the horizon (i.e. the horizon location is fixed). It implies that, on the horizon

N = B i si . (15)
One can also show (see for instance [START_REF] Vasset | Excised black hole spacetimes: Quasi-local horizon formalism applied to the Kerr example[END_REF]) that the shear σ ab of the null light rays must vanish. It implies that the shift vector B i must be a conformal Killing vector of the sphere.

A possible choice together with Eq. ( 15) is

B i = N si -Ω BH m i , (16)
where Ω BH is a constant that captures the rotation velocity of the black hole and

m i = (∂ ϕ ) i .
The minus sign ensures that the angular momentum is positive.

Equations [START_REF] Headrick | A New approach to static numerical relativity, and its application to Kaluza-Klein black holes[END_REF][START_REF] Adam | A numerical approach to finding general stationary vacuum black holes[END_REF][START_REF] Wiseman | Numerical construction of static and stationary black holes[END_REF] have been successfully used in many publications, especially for computing binary black holes initial data [START_REF] Caudill | Circular orbits and spin in black-hole initial data[END_REF][START_REF] Uryu | New code for equilibriums and quasiequilibrium initial data of compact objects. II. Convergence tests and comparisons of binary black hole initial data[END_REF][START_REF] Papenfort | New public code for initial data of unequal-mass, spinning compact-object binaries[END_REF]. However in those papers the spatial metric is set to be conformally flat, an assumption that is not used in this work.

C. Differential gauges

Maximal slicing and the spatial harmonic gauges are differential gauges. By that, it is to be understood that they do not impose conditions directly on the metric fields but rather lead to partial differential equations. When regularity of the full spacetime is required, it leads to a unique choice of coordinates. However, as will be illustrated below, when a horizon is present, the coordinate system is fixed up to some boundary conditions which can be used to freely specify some quantities.

Consider an infinitesimal coordinate change of time of the form t = t + α (x i ). At first order it does induce a change on the lapse and the shift proportional to ∂ i α. Demanding that the new coordinate system obeys the maximal slicing condition K = 0 leads to a partial differential equation which is second order in terms of the coordinate change α, K containing first order derivatives of the shift. So in order to transform an arbitrary coordinate system into one with maximal slicing, one needs to solve a second order partial differential equation for α. It follows that α is determined up to two boundary conditions. The one at infinity is implicitly accounted for by demanding that the metric takes a fixed form (flat one in Sec III and V or ADS in Sec. IV). In the presence of a horizon, the inner boundary condition for α can translate in the free choice of the lapse. So on the horizon the lapse is a freely specifiable angular function N 0 (θ, ϕ).

The situation concerning the spatial coordinates is similar. Considering a coordinate change of the form x i = x i + ξ i (x j ), one can show that the spatial harmonic gauge equation ( 6) is of second order in terms of ξ i . Indeed the Christoffel's symbols contain first order derivatives of the metric which, in turn, contains first order derivatives of ξ i . However not all of the components of ξ i can be chosen freely. Remember that the location of the horizon has been chosen beforehand as being a sphere of constant radius. In order to maintain this location, one needs to have ξ r = 0. The two angular components of the coordinate change, however, can be freely chosen on the horizon and this choice translates into the possible free choice of some of the 3+1 quantities. The shift being fixed by Eq. ( 16), it is more convenient to enforce the value of some components of the metric. The most natural choice is to fix the components γ rθ and γ rϕ of the spatial metric (with the spectral methods used in this paper those components have the same spectral bases as ξ θ and ξ ϕ , leading to a well posed numerical system ; this is probably more profound than just technicalities of the numerical method). So, on the horizon one sets γ rθ = f and γ rϕ = g, where f and g are arbitrary angular functions.

D. Degeneracy of the equations

In the context of this work, an equation is said to be degenerate if the prefactor of the highest order derivative vanishes. As an illustration, consider the following equation

a (x) f +b (x) f +c (x) f +d (x) = 0 on [-1, 1]. If a (x)
does not vanish, then it can be solved with the imposition of two boundary conditions, at x = -1 and x = 1. If a (x = -1) = 0, the equation is degenerate and one can no longer choose any boundary condition at x = -1.

Indeed, at this point, the equation reduces to b (x) f + c (x) f + d (x) = 0 which is the only compatible choice and must be used as a boundary condition. In a sense, the equation is its own boundary condition.

This type of behavior is present in the 3+1 equations considered here and it must be dealt with carefully. This can be seen in the evolution equations [START_REF] Lisamission | Lisa consortium website[END_REF] where the factor in front of the second order radial derivatives of the metric is (B r ) 2 -N 2 γ rr /2N which vanishes on the horizon, given Eq. ( 15). However, as the full set of equations is coupled, one can not simply assign one equation to one unknown. The system should be considered as a whole.

In order to do so, one needs to isolate, in Eqs. [START_REF] Abuter | Detection of the Schwarzschild precession in the orbit of the star S2 near the Galactic centre massive black hole[END_REF][START_REF] Akiyama | Event Horizon Telescope), First M87 Event Horizon Telescope results. I. The shadow of the supermassive black hole[END_REF][START_REF] Lisamission | Lisa consortium website[END_REF] the terms involving the second radial derivatives of the fields.

                         0 0 0 0 -1 2 γ rr 0 1 2N (1 + γ rr γ rr ) 1 2N γ rr γ rθ 1 2N γ rr γ rϕ 1 2N γ rr B r 0 1 2N γ rθ γ rr 1 2N γ rr γ θθ 1 2N γ rr γ θϕ 0 0 1 2N γ rϕ γ rr 1 2N γ rr γ θϕ 1 2N γ rr γ ϕϕ 0 -1 1 N B r γ rr 1 N B r γ rθ 1 N B r γ rϕ 1 2N (B r ) 2 -N 2 γ rr 0 1 2N B r γ rθ 1 2N B r γ θθ 1 2N B r γ θϕ 0 0 1 2N B r γ rϕ 1 2N B r γ θϕ 1 2N B r γ ϕϕ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (17) -1 2 γ rθ -1 2 γ rϕ -1 2 γ θθ -1 2 γ θϕ -1 2 γ ϕϕ 1 2N γ rθ B r 1 2N γ rϕ B r 0 0 0 1 2N γ rr B r 0 1 2N γ rθ B r 1 2N γ rϕ B r 0 0 1 2N γ rr B r 0 1 2N γ rθ B r 1 2N γ rϕ B r 0 0 0 0 0 1 2N (B r ) 2 -N 2 γ rr 0 0 0 0 0 1 2N (B r ) 2 -N 2 γ rr 0 0 0 0 0 1 2N (B r ) 2 -N 2 γ rr 0 0 0 0 0 1 2N (B r ) 2 -N 2 γ rr 0 0 0 0 0 1 2N (B r ) 2 -N 2 γ rr                         
This can be done and the result is given by [START_REF] Vasset | Excised black hole spacetimes: Quasi-local horizon formalism applied to the Kerr example[END_REF]. The lines correspond to the equations of the system (i.e. line 1 to Eq. ( 8), line 2 to 4 to Eq. ( 9) and lines 5 to 10 to Eq. ( 10)).

The columns correspond to the metric fields (i.e. column 1 to N , columns 2 to 4 to B i and columns 5 to 10 to γ ij ). The quantities in [START_REF] Vasset | Excised black hole spacetimes: Quasi-local horizon formalism applied to the Kerr example[END_REF] are then the factors of the terms ∂ 2 rr of a given field, in a given equation.

As already noted, the factor appearing on the last lines vanishes on the horizon. Though it is difficult to find explicitly the eigenvalues of [START_REF] Vasset | Excised black hole spacetimes: Quasi-local horizon formalism applied to the Kerr example[END_REF], one can investigate them numerically after assigning some random but realistic values to the fields on the horizon. It appears that the multiplicity of the null eigenvalue is always three. It follows that the only degenerate equations correspond to the last three lines of [START_REF] Vasset | Excised black hole spacetimes: Quasi-local horizon formalism applied to the Kerr example[END_REF], which are the purely angular components of Eq. (10) (i.e. the components (θ, θ), (θ, ϕ) and (ϕ, ϕ)). Those components must be solved without any boundary conditions or equivalently as being their own boundary conditions.

E. Behavior of the expansion

At this point, it seems that there are enough boundary conditions to solve the problem. not impossible to achieve, the numerical system seems to exhibit some instabilities. It is observed that the code is much more stable when the boundary condition Θ = 0 is relaxed and replaced by the imposition of the value of γ rr on the horizon. One can then monitor the value of the expansion and check that it is indeed zero. This is shown in Fig. 1 where the value of Θ on the horizon is shown, as a function of the resolution, for three different values of γ rr on the horizon. When precision increases Θ goes to zero, for all three different values of γ rr . It shows that all the configurations correspond to valid non-rotating black holes, the different choices of γ rr on the horizon corresponding to different masses.

Indeed, one could use

N = N 0 , B i = N si -Ω BH m i , Θ = 0, γ rθ = f , γ rϕ = g and E θθ = E θϕ = E ϕϕ = 0.
The fact that there is no need to enforce directly Θ = 0 in the non-rotating case can actually be understood analytically. When Ω BH = 0, the only non-vanishing component of the shift is B r and the spatial line element reduces to ds 2 = Adr 2 + BdΩ 2 . All the quantities depend only on the radial coordinate. Given those expressions, one can compute the expansion and the trace of K ij and find that

Θ = rBB r ∂ r A + 2rAB∂ r B r + 2 (rN ∂B + 2BN ) √ A 2rABN (18) K = rBB r ∂ r A + 2rAB∂ r B r + 2rAB r ∂ r B + 4ABB r 2rABN . ( 19
)
Inserting K = 0 in Eq. ( 18) then leads to

Θ = N √ A -AB r (2r∂ r B + 4B) 2rABN , (20)
which is indeed zero on the horizon given the boundary condition Eq. [START_REF] Adam | A numerical approach to finding general stationary vacuum black holes[END_REF]. So it follows that if the system of equations ensures that K = 0, then Θ = 0 is also guaranteed. If the procedure presented in Sec. II A used to enforce maximal slicing works properly, then there is no need to explicitly impose that Θ = 0. Given the curves shown in Fig. 1, it appears to be the case.

When rotation is present (i.e. when Ω BH = 0) the situation is slightly more complicated.

If one tries the same procedure as in the non-rotating case, one can show that fixing γ rr on the horizon does not lead to Θ = 0. Moreover one can observe that the error does not contain any spherically symmetric component (for instance, if projected onto spherical harmonics, there is no component on Y 0 0 ). So in order to maintain some kind of continuity between the non-rotating and rotating cases, one is led to consider the following choice of boundary conditions:

• the spherical part of γ rr is chosen arbitrarily.

• the non-spherical part of Θ = 0 must be solved.

In this paper, spectral methods are used so that the splitting between spherical and nonspherical parts is essentially straightforward. Should other numerical methods be used, this may not be as simple but this is beyond the scope of this work.

F. Complete set of boundary conditions

The full set of boundary conditions on the horizon (here a sphere of fixed radius) is given by:

N = N 0 (θ, ϕ) (21)
B i = N si -Ω BH m i Spherical part : γ rr = γ 0 Non -spherical part : Θ = 0 γ rθ = f (θ, ϕ) γ rϕ = g (θ, ϕ) E θθ = 0 E θϕ = 0 E ϕϕ = 0.
In those equations, there are three freely specifiable angular functions N 0 , f and g and one free number γ 0 . A standard choice for those values is: N = 1/2, f = g = 0 and γ 0 = 8. Unless otherwise stated, this is what is used throughout this paper. The equations are solved numerically, using the Kadath library [START_REF] Grandclement | Kadath: A Spectral solver for theoretical physics[END_REF][START_REF] Grandclement | Kadath libray[END_REF]. This tool relies on spectral methods to solve systems of partial differential equations and it has been successfully applied to the study of various problems in general relativity and theoretical physics.

Let us recall that, once the equations are solved, it needs to be checked a posteriori that the gauge conditions K = 0 and V i = 0 are indeed verified. As already stated this is a very important test. There is also a need to verify that the spherical part of Θ vanishes, as it is not explicitly enforced by the set of boundary conditions Eqs. [START_REF] Martinon | Gravitational geons in asymptotically anti-de Sitter spacetimes[END_REF].

III. KERR BLACK HOLES

The simplest and most straightforward application of the equations presented above is the computation of a single rotating black hole in general relativity. It must lead to the famous Kerr spacetime [START_REF] Kerr | Gravitational field of a spinning mass as an example of algebraically special metrics[END_REF]. However, with the gauge choices used in this paper, the solution is found in coordinates that are not analytical.

The equations solved are (8-10) with the inner boundary conditions [START_REF] Martinon | Gravitational geons in asymptotically anti-de Sitter spacetimes[END_REF]. Λ is set to zero and there is no matter. The system is closed by demanding that flat spacetime is recovered at spatial infinity. It simply implies that N = 1, B i = 0 and γ ij = f ij , where f ij denotes the flat metric.

For the single black hole problem, the numerical spacetime is decomposed into several (typically 4) spherical shells. The last domain extends up to infinity by means of the variable 1/r so that boundary conditions are enforced at exact spatial infinity. The solutions are found iteratively by means of a Newton-Raphson iteration. The first computed configuration is the Schwarzschild one, for which Ω BH = 0. The angular velocity is then incremented in order to compute a sequence of rotating black holes with different Kerr parameters.

Once a given configuration has been computed, various global, coordinate-independent quantities can be computed. The ADM (Arnowitt-Deser-Misner) mass is given by a surface integral at infinity:

M ADM = 1 16π r=∞ f ik f jl Dj γ kl -Dk γ jl dS, ( 22
)
where D denotes the covariant derivative associated with the flat metric f ij and dS is the surface element at infinity. The spacetime being stationary, one can also define the Komar mass of the system by

M Komar = 1 4π r=∞ si D i N -K ij si sj dS. ( 23
)
Given Eq. ( 16), rotation is around the z-axis only and the angular momentum is given by

J = 1 8π r=∞ K ij m i sj dS. ( 24
)
Various error indicators are monitored in Fig. 2. The quantities are plotted as a function of the resolution. This corresponds to the number of coefficients of the spectral expansion, with respect to the coordinates r and θ (in that case the same number is used for both dimensions). Convergence of the various indicators is shown for a configuration with r H = 1

and Ω BH = 0.1, which corresponds to a mass M ADM ≈ 2.18 and a Kerr parameter a/M ≈ 0.73. The first panel shows that the quantities K and V i decrease exponentially as resolution increases. As already stated, this is an important test that ensures that the system solved coincides with Einstein's equations and that the gauge choices are indeed fulfilled. The second panel of Fig. 2 shows the maximal value of Θ (Eq. 14) on the horizon. As seen in Sec. II E the spherical part of Θ = 0 is not solved numerically and this curve shows that it is indeed verified as it goes to zero exponentially. The last two panels involve the computation of global quantities and so the errors are slightly bigger. Indeed the computation of those quantities involve surface integrals at infinity that introduce additional numerical errors, when compared to quantities like K or V i . This can explain the fact that convergence is less regular, as seen with the values for a resolution of 25. The saturation level is also somewhat higher, with a value of about 10 -7 -8 . This is not surprising as the Newton-Raphson algorithm was stopped at a threshold of 10 -8 .

That being said, the last two panels of Fig. 2 still show a good convergence with resolution.

The third one shows the relative difference between the ADM and Komar masses, as it is known that the two must be equal in that case (see Sec. 8.6.2 of [START_REF] Gourgoulhon | 3+1 Formalism in General Relativity[END_REF] for more details on this equality). The last panel shows the relative difference from the numerical angular velocity Ω BH and the one computed from M ADM and a by the analytic expression

Ω = a 2M M + √ M 2 -a 2 . ( 25
)
Using the highest resolution at hand [START_REF] Grandclement | Kadath: A Spectral solver for theoretical physics[END_REF] technical difficulty stems from the fact that Ω BH is not a monotonic function of the Kerr parameter, when the radius r H is fixed. This is clearly seen in the first panel of Fig. 3.

In order to be able to pass the maximum, the value of Ω BH can be made an unknown of the numerical problem and the condition that a/M has a given value can be added to the system. The Kadath library enables the use of such global unknowns (i.e. unknowns that are not fields but numbers). Doing so one can reach high values of the Kerr parameters, that correspond to small values of Ω BH . In the second panel the angular velocity is also shown, as a function of a/M , but this time scaled with the ADM mass. The fact that the mass is not constant along the sequence (it is r H that is fixed), explains the different behavior between the two panels. As is expected (see Eq. ( 25)), M Ω BH goes to 1/2 when one gets closer to the extremal case. Let us finally mention that very high values of a/M ≈ 0.99 can be reached without much trouble and while maintaining an accuracy of about 10 -7 . If need be, even higher values could be computed.

As an illustration, in Fig. 4, various contours of some fields are shown, in the xz-plane, for the configuration with the highest value of a/M ≈ 0.99.

IV. MTZ BLACK HOLES

In this section the formalism presented in Sec. II is applied to a class of non-rotating black holes with a negative cosmological constant and a minimally coupled scalar field. This solution was obtained analytically by Martínez, Troncoso and Zanelli (MTZ) in [START_REF] Martínez | Exact black hole solution with a minimally coupled scalar field[END_REF].

A. Analytic solution and adjustments to the system of equations

The field equations are the Einstein ones with a negative cosmological constant Λ = -3/ 2 ADS , as well as the Klein-Gordon equation for a real scalar field φ, minimally coupled to gravity.

The stress-energy tensor is given by

T µν = ∇ µ φ∇ ν φ - 1 2 g µν g αβ ∇ α φ∇ β φ -g µν V (φ), (26)
and the interaction potential is

V (φ) = - 3 4π 2 ADS sinh 2 4π 3 φ. ( 27
)
The 3+1 matter terms in that particular case are given in Appendix A 1.

In [START_REF] Martínez | Exact black hole solution with a minimally coupled scalar field[END_REF], a family of black hole solutions parameterized by an integration constant µ (linked to the black hole mass) is found :

ds 2 = r(r + 2µ) (r + µ) 2 - r 2 2 ADS -1 + µ r 2 dt 2 + r 2 2 ADS -1 + µ r 2 -1 dr 2 + r 2 dσ 2 (28) φ = 3 4π Arctanh µ r + µ . ( 29
)
The term dσ 2 represents the line element of a two-dimensional manifold with constant negative curvature. In this work, it is chosen as:

dσ 2 = 1 cos 2 θ (dθ 2 + sin 2 θ dϕ 2 ). (30)
This is to be contrasted with [START_REF] Martínez | Exact black hole solution with a minimally coupled scalar field[END_REF] where this two-dimensional manifold (and hence the horizon) is assumed to be compact. The form [START_REF] Gourgoulhon | A 3+1 perspective on null hypersurfaces and isolated horizons[END_REF] corresponds to coordinates on the hyperbolic plane.

Furthermore, the vector field (∂ ϕ ) i appearing in the boundary conditions ( 21) is not a conformal Killing field of the surface described by [START_REF] Gourgoulhon | A 3+1 perspective on null hypersurfaces and isolated horizons[END_REF]. However, as the solution represents a non-rotating black hole, for which Ω BH = 0, this is not an issue.

For r -→ ∞ or equivalently µ = 0, anti-de Sitter spacetime with constant negative (4D) curvature 4 R = 2Λ is recovered. The solution is said to be asymptotically anti-de Sitter (AADS).

The constant µ is bounded from below by µ > ADS 4 and r must be larger than the maximum of 0 and -2µ. Under those conditions, there is a horizon located at r = r + with

r + = ADS 2 1 + 1 + 4µ ADS . (31)
The corresponding apparent horizon is a surface of constant negative curvature.

In order to describe surfaces of negative curvature, one considers the following reference spatial metric:

γ ref ij dx i dx j = dr 2 + r 2 dσ 2 = dr 2 + r 2 cos 2 θ (dθ 2 + sin 2 θdϕ 2 ). ( 32
)
This reference metric is used when defining the Christoffel symbols of the real metric. More precisely, what is computed numerically is the difference between the Christoffel symbols of γ ij and those of γ ref ij . This differs from more usual situations where the reference metric is the flat one. The numerical description of the reference metric is done via the basis:

e r = ∂ r ; e θ = cos θ r ∂ θ ; e ϕ = cos θ r sin θ ∂ ϕ , (33)
which makes its orthonormal.

As the solutions are not asymptotically flat but AADS, the ADS boundary must be dealt with carefully, as some metric quantities diverge at this location. In order to do so, a particular form of the ADS spacetime is chosen

ds 2 = - cos 2 log R ADS r sin 2 log R ADS r dt 2 + 2 ADS r 2 sin 2 log R ADS r γ ref ij dx i dx j . ( 34
)
R ADS is a freely specifiable value that gives the position of the boundary of the ADS spacetime, which is located at r = R ADS . The form ( 34) is chosen because the spatial metric relates conformally to the reference one

γ ADS ij = 2 ADS r 2 sin 2 log R ADS r γ ref ij . (35)
The ADS spatial metric is the one which is used in defining the spatial harmonic gauge, meaning that γij = γ ADS ij (see Eq. ( 6)).

From Eq. ( 34) one can see that the lapse and spatial metric diverge at this boundary. In order to allow for a numerical treatment there is a need to regularize the divergences. This is done by defining the conformal factor

Ω ≡ sin log R ADS r . (36)
Ω can then be used to define the following regularized quantities (denoted by a tilde):

Ñ ≡ ΩN (37) Bi ≡ B i (38) γij ≡ Ω 2 γ ij (39) φ ≡ φ. ( 40
)
The outer boundary conditions, where r = R ADS , are then Ñ = 1, Bi = 0, γij =

R 2 ADS γ ref ij and φ = 0.
All the quantities appearing in the equations must be regularized near the ADS boundary, along with the equations themselves (see Appendix A 2 for explicit expressions).

This regularization procedure applies in the outermost numerical domain, which is a spherical shell extending up to the ADS boundary at r = R ADS . Proper continuity of the fields and their radial derivatives across the boundary with the inner domain is enforced.

This translates into the following non trivial conditions (where superscripts (I) and (O) stand for inner and outer respectively):

∂ r Ñ (O) = N (I) ∂ r Ω + Ω∂ r N (I) , (41)
∂ r γ(O) ij = Ω 2 ∂ r γ (I) ij + 2Ω(∂ r Ω)γ (I) ij . (42)
Last but not least, it must be noted that the the Klein-Gordon equation is degenerate on the horizon. From its expression, given by Eq. (A5), one can notice that the only contribution to the principal symbol for double r partial derivatives is γ rr -

(B r ) 2 N 2 ∂ 2 rr φ.
It is the same factor that appears in the Einstein evolution equations (see Sec. II D) and it vanishes on the horizon with the boundary condition [START_REF] Adam | A numerical approach to finding general stationary vacuum black holes[END_REF]. Therefore, the Klein-Gordon equation is degenerate on the horizon and there is no need to impose any additional boundary condition on the scalar field itself.

B. Results

The cosmological constant is chosen so that ADS = 20. Space is split into two spherical shells. The outer one goes from r = 30 to r = R ADS = 40 and so extends up to the ADS boundary. The radius of the horizon is varied in order to compute different configurations.

At the inner boundary, which is the apparent horizon, one sets N 0 = 1/2, γ 0 = 1 and f = g = 0 (see Eqs. ( 21)). Recall that the MTZ solution is a static one so that Ω BH = 0.

As in the Kerr black hole case, spectral convergence of the gauge quantities and of the expansion on the horizon is monitored. Figure 5 shows these quantities as a function of the radial resolution and convergence is clear. The saturation seen at 10 -8 is due to the Newton-Raphson solver that is stopped at this level. Profiles of various fields are shown in Fig. 6 for a configuration with r H = 8. Divergence of the lapse and of the spatial metric at the ADS boundary is noticeable (remember that this boundary is at R ADS = 40). The smoothness of the curves of N and γ ij illustrates the good behavior of the regularization procedure.. In order to compare the numerical results with the solution (28), coordinate-independent quantities must be compared. The four-dimensional Ricci scalar on the horizon 4 R H (φ H ) is one possibility. It can be computed in the 3+1 formalism by Eq. (3.75) [START_REF] Gourgoulhon | 3+1 Formalism in General Relativity[END_REF], written here in the case K = 0

4 R = R + K ij K ij - 2 N D i D i N. (43)
A sequence of MTZ black holes with various radii for the horizon, ranging between r H = 6.5 and r H = 10 is computed. Although each solution cannot directly be linked to a specific value of µ, the sequence consists of black holes with different masses. The four-dimensional Ricci scalar on the horizon, as a function of the value of the scalar field on the horizon is shown in Fig. 7. The circles denote the numerical results and the solid curve the analytical one. A very good agreement between the two is achieved. The smaller the horizon radius, the more intense the scalar field on the horizon and the larger the absolute value of the curvature on the horizon. This means that small radii correspond to more relativistic configurations.

V. BLACK HOLES WITH SCALAR HAIRS

In this section, the construction of Kerr black holes with scalar hairs is explained. This is the same system as computed by [START_REF] Herdeiro | Kerr black holes with scalar hair[END_REF][START_REF] Herdeiro | Construction and physical properties of Kerr black holes with scalar hair[END_REF] but using a different choice of coordinates.

A. Equations for the scalar field

The class of hairy black holes constructed here relies on the existence of a complex scalar field Φ. This field is minimally coupled to gravity and solutions containing a horizon can be found. The scalar field obeys the Klein-Gordon equation ∇ µ ∇ µ Φ = µ 2 Φ where ∇ denotes the covariant derivative of the four-dimensional metric. The right-hand-side corresponds to the case of a free massive field which has a potential V |Φ| 2 = µ 2 |Φ| 2 , where µ is the mass of the field.

In previous works the field is assumed to have the form Φ = φ (r, θ)

exp [i (ωt -kϕ)],
where ω is the angular velocity and k an integer dubbed the rotational quantum number.

Given the expression of the action (see [START_REF] Herdeiro | Kerr black holes with scalar hair[END_REF][START_REF] Herdeiro | Construction and physical properties of Kerr black holes with scalar hair[END_REF]), which has a U (1) symmetry, the resulting spacetimes are axisymmetric (the two Killing vectors being (∂ t ) µ and (∂ ϕ ) µ ) and the quantities ω and k appear as parameters of the solutions. When the ansatz is inserted into the Klein-Gordon equation, it leads to an expression of the form

(R KG + iI KG ) exp [i (ωt -kϕ)] = 0. ( 44
)
In [START_REF] Herdeiro | Kerr black holes with scalar hair[END_REF][START_REF] Herdeiro | Construction and physical properties of Kerr black holes with scalar hair[END_REF] the part R KG is the only one considered and it is the equation fixing the value of the amplitude φ. Indeed, given the coordinates used, and in particular the fact that the only non-vanishing component of the shift is B ϕ , it is easy to show that I KG is identically zero.

With the coordinate system introduced in this work, this is no longer the case (basically because the shift has a non-vanishing radial component). It means that the original ansatz does not pass through the equation and that a more general form must be used.

One considers the following form for the scalar field, introducing an additional imaginary component

Φ = (R Φ (r, θ) + iI Φ (r, θ)) exp [i (ωt -kϕ)] = 0. (45)
Inserting (45) in the Klein-Gordon equation and using the 3+1 decomposition of spacetime, one can find the expressions for R KG and I KG in terms of R Φ and I Φ . One gets

R KG = (ω + kB i δ ϕ i ) 2 N 2 R Φ + 1 N D i N γ ij - B i B j N 2 D j R Φ + γ ij - B i B j N 2 D i D j R Φ (46) -γ ij k 2 δ ϕ i δ ϕ j R Φ -µ 2 R Φ -ω B i N 2 D i I Φ - ω N D i B i N I Φ + 1 N D i N γ ij - B i B j N 2 kδ ϕ j I Φ + γ ij - B i B j N 2 kδ ϕ i D j I Φ ,
where δ ϕ i is the gradient of ϕ. Expressed in the spherical orthonormal basis used in this paper, it is given by δ ϕ i = 0, 0, 1 r sin θ . Divisions by sin θ -which vanishes on the z-axisare performed using the coefficients of the spectral expansion of the fields, in order to avoid divergences.

Eq. ( 46) is of the form A (R Φ ) + B (I Φ ) = 0 where A and B are linear differential operators of the scalar field (A is second order and B first order). Not surprisingly, given the ansatz used, the expression of I KG is very similar to the one for R KG . One finds that

I KG = A (I Φ ) -B (R Φ ).
As a first step, the metric fields are fixed and only the Klein-Gordon equation is solved.

It means that the backreaction of the field on the metric is neglected. This is known as a cloud solution. For the metric fields one uses the Kerr black hole configurations computed in Sec. III. Let us mention that, following [START_REF] Herdeiro | Kerr black holes with scalar hair[END_REF][START_REF] Herdeiro | Construction and physical properties of Kerr black holes with scalar hair[END_REF], the angular velocity of the black hole and of the field are linked by ω = kΩ BH , a condition that prevents the field from having a flux across the horizon. This condition is enforced throughout this paper.

The resolution of this linear problem needs to be dealt with carefully. First one can notice that the part of R KG that contains second order radial derivatives is

γ rr -B r B r N 2 ∂ 2 r R Φ .
As already seen in Sec. II and IV, with the boundary conditions used, this term vanishes on the horizon. The same is true for I KG , with respect to I Φ . It means that those equations do not require any inner boundary condition (see Sec. II D). So there is no need to enforce anything on the fields at the horizon.

Two other properties of the linear Klein-Gordon equations must be taken into account.

First notice that R Φ = I Φ = 0 is a solution so that one needs to prevent the code from converging to this trivial configuration. Also, it is expected that solutions can be found only for a discrete set of the physical parameters (k, Ω BH , µ in particular). The numerical procedure should be able to find those values. Those two features are typical of linear systems and have been dealt with successfully in the context of spacetimes with cosmological constants in [START_REF] Fodor | Scalar field breathers on anti-de Sitter background[END_REF][START_REF] Fodor | Self-gravitating scalar breathers with negative cosmological constant[END_REF] and the reader should refer to those publications for details about the procedure. It is only briefly sketched below.

The equations being second order, they should be solved by demanding the matching of the fields and their normal derivative across the boundaries of the various numerical domains. For one such boundary, the continuity of one angular spectral coefficient is relaxed and replaced by the condition that this coefficient has a given value (1 for instance). By construction, this prevents the code from going to the trivial solution. The value of the coefficient itself is unimportant, the problem being linear. However, in general, that would lead to a solution with a discontinuity on the derivative. The next step is then to scan the possible values of the parameter space and it appears that for some values, the error on the discontinuity vanishes. An example of that is shown in Fig. 8 where the error on the derivative of R Φ is plotted, as a function of the parameter µ, around the value for a true solution. The correct value of µ can be determined by a dichotomy algorithm. However, a good precision is required, as the error varies very fast with µ (notice the abscissa range in Fig. 8). As observed in [START_REF] Fodor | Scalar field breathers on anti-de Sitter background[END_REF][START_REF] Fodor | Self-gravitating scalar breathers with negative cosmological constant[END_REF], there are several possible values of the parameters that are admissible. They correspond to different number of nodes of the amplitude |Φ|. Only nodeless configurations are considered here.

An additional difficulty arises from the fact that the quantities R Φ and I Φ are determined up to a constant phase ; the fields defined as R Φ = R Φ cos α -I Φ sin α and I Φ = R Φ sin α + I Φ cos α are also a valid solution of the system. The value of this phase must be enforced when solving numerically the equations, otherwise the Newton-Raphson iteration would fail.

To do so, one follows that same procedure for I Φ as for R Φ : the continuity of the derivative is relaxed and replaced by a condition on the value of one coefficient of I Φ , at an arbitrary location. It appears that different values of I Φ all lead to valid solutions that correspond to different choices of α. As this technique is also implemented in the full system case (i.e. the one with gravitation), more details are given about this in Sec. V C.

The cloud solutions described in this section are only used as an initial guess to get the configurations of the coupled Einstein-Klein-Gordon system so that they are not investigated much here. However, their validity has been carefully checked by comparing with results coming from the direct resolution of the Teukolsky equation as done in [START_REF] García | Obstructions towards a generalization of no-hair theorems: Scalar clouds around Kerr black holes[END_REF] for instance.

B. Coupling with gravity

When the scalar field is coupled to gravity, one must take into account the stress-energy tensor of the field, which reads

T µν = ∇ (µ Φ∇ ν) Φ -1 2 g µν ∇ µ Φ∇ µ Φ + V |Φ| 2 .
Using the ansatz of the field (45) one can split the time and space components of the stressenergy tensor and from that deduce the 3+1 matter terms entering Einstein's equations.

Expressions of E, P i and S ij , as a function of the scalar field components R Φ and I φ are given in Appendix B.

Equations [START_REF] Abuter | Detection of the Schwarzschild precession in the orbit of the star S2 near the Galactic centre massive black hole[END_REF][START_REF] Akiyama | Event Horizon Telescope), First M87 Event Horizon Telescope results. I. The shadow of the supermassive black hole[END_REF][START_REF] Lisamission | Lisa consortium website[END_REF] are solved using the inner boundary conditions [START_REF] Martinon | Gravitational geons in asymptotically anti-de Sitter spacetimes[END_REF] and demanding that, at spatial infinity, flat spacetime is recovered.

Concerning the Klein-Gordon equation, the situation is different from the linear case exposed in Sec. V A. Indeed, with the coupling with gravity, the problem is no longer linear and admits solutions not for discrete values of the parameters but on a whole continuum.

It follows that the whole procedure about R Φ described in Sec. V A is irrelevant. However, the discussion about the arbitrary phase in the definition of R Φ and I Φ still holds and the same procedure as in V A is used.

The choice of the initial guess is of uttermost importance and that is where the knowledge of the cloud solutions of Sec. V A is needed. Consider a linear solution of the Klein-Gordon equation, corresponding to a set of parameters. When moving away from this solution, by changing the value of some of the parameters, one can construct a sequence of nonlinear solutions. The further away the parameters are from the linear values, the higher the amplitude of the scalar field. In this work, a sequence is constructed by varying the parameter Ω BH , maintaining the equality ω = kΩ BH . This is only a choice and other parameters could be varied (µ for instance).

In order to use the cloud solution as an initial configuration, one proceeds as follows. A small amplitude for the maximum of the field is chosen and the cloud solution is scaled so that it has this maximum value. The parameter Ω BH is considered as an unknown of the numerical problem. By this it is to be understood that the solver is allowed to change its value. As the system contains an additional unknown it needs to be supplemented with an additional condition. This condition is the fixing of the amplitude of the field. By increasing the amplitude of the field, different configurations can be computed, corresponding to different values of Ω BH . Alternatively, once the first full solutions are known, one can directly vary the parameter Ω BH along the sequence, in order to compute solutions with different amplitude of the field. The steps must be small enough so that the code does not converge to the trivial solution where the field vanishes everywhere.

C. Numerical results

A small modification in the numerical setting stems from the fact that the scalar field and black hole sizes are somewhat different. In the computations, the black hole is located at a coordinate radius of r H = 1. However, the toroidal shape of the scalar field has a maximum at about 30 times this value. In order to deal with those two scales, more numerical domains than in Sec. III are required. Typically, for each shell, the outer radius is twice the inner radius and a dozen of domains are used.

As already mentioned, a sequence of diverse values of Ω BH is exhibited. The starting point of the sequence is a cloud solution (see V A) constructed from a Kerr black hole with M ≈ 5 and a/M ≈ 0.95, which corresponds to Ω BH ≈ 0.0724. The rotational quantum number is chosen to be k = 1.

As in previous cases, one must check that the gauge quantities K and V i converge to zero when resolution increases. The value of the expansion [START_REF] Headrick | A New approach to static numerical relativity, and its application to Kaluza-Klein black holes[END_REF] on the horizon should also go to zero. In order to do so, an arbitrary configuration is chosen and computed with different As explained in Sec. V A, the quantities R Φ and I Φ are defined up to a constant phase. The continuity of the radial derivative of I Φ is relaxed, for one spectral coefficient. The matching of the radial derivative is replaced by a condition fixing the value of that same coefficient, for the field itself. It amounts to fixing the constant phase. In the following two possible choices are considered: one can demand that the coefficient of the field I Φ vanishes (referred to as the zero phase condition) or that the coefficient of I Φ is half that of R Φ (referred to as the half phase condition).

For the solution with Ω BH = 0.067, the fields resulting from those two cases are shown in Fig. 10 where profiles of the fields, along the x-axis are shown. As expected the decomposition of the complex scalar field in terms of R Φ and I Φ changes. However, it is easy to check that the two relate by a constant phase rotation, the parameter being α ≈ -0.464 in that particular case. It is also easy to check that the two configurations lead to the same global quantities (mass, charge etc...) as they should.

The continuity of the radial derivative of I Φ can also be used to monitor the validity of the procedure. As it has not been enforced numerically, it needs to be checked that it is actually satisfied by the solution. Figure 11 shows the error on the continuity of the radial derivative of I Φ , as a function of resolution. As expected with spectral methods, the error converges fastly to zero confirming that the solution is smooth as it should. Even if they are not explicitly shown here, the solutions passed several other tests (like equality of the

ADM and Komar masses for instance).

The ADM mass and angular momentum of the solutions can be computed using the standard formulae [START_REF] Andersson | Elliptic hyperbolic systems and the Einstein equations[END_REF] and [START_REF] Figueras | Ricci solitons, Ricci flow, and strongly coupled CFT in the Schwarzschild Unruh or Boulware vacua[END_REF]. As one moves along the sequence by going away from the Kerr black hole solution, those quantities increase, as the contribution of the scalar field to their values is more important. This can be seen in Fig. 12, where the mass and momentum are shown, as a function of Ω BH . The last point to the right (i.e. with Ω BH ≈ 0.0724) corresponds to the Kerr black hole from which the sequence is constructed.

In order to monitor the relative importance of the black hole compared to the scalar field, following [START_REF] Herdeiro | Kerr black holes with scalar hair[END_REF][START_REF] Herdeiro | Construction and physical properties of Kerr black holes with scalar hair[END_REF], one can compute the Noether charge. With the setting used in this paper it reads (see Sec. B for a precision about the use of the notations RΦ and ĪΦ )

Q = - Σt 1 N ω + kδ ϕ i B i R Φ RΦ + I Φ ĪΦ (47) + B i 2 I Φ D i RΦ + ĪΦ D i R Φ -R Φ D i ĪΦ -RΦ D i I Φ √ γd 3 x.
From the Noether charge one defines the quantity q = kQ/J. For a Kerr black hole q = 0 as the Noether charge vanishes. For a boson star, it is known that J = kQ, so that q = 1 [39]. It follows that the quantity q measures the relative importance of the hole and the field along the sequence, as it starts from 0 for the Kerr black hole and approaches 1 as the amplitude of the scalar field increases. The quantity q as a function of Ω BH is shown in the first panel of Fig. 13.

The Komar mass can also be used as a way of measuring the influence of the field. Indeed, if it can be computed as a surface integral at infinity by Eq. ( 23), it can also be obtained by the sum of two terms, M BH a surface integral on the horizon and M Φ a volume integral over matter terms. The general expressions for M BH and M Φ can be found in [START_REF] Gourgoulhon | 3+1 Formalism in General Relativity[END_REF] (Eqs. (8.70) and (8.71)) respectively. In the situation at hand they are obtained as

M BH = 1 4π r=r h √ h si D i N -K ij si B j dS (48) M Φ = Σt √ γ N E + γ ij S ij -2P i B i d 3 x, (49)
where h denotes the determinant of the metric induced on the horizon h ij = γ ij -si sj . In the second plot of Fig. 13 those quantities are shown, as a function of Ω BH . They are scaled by the total mass so that the sum is one (given that M Komar = M ADM ). As for the quantity q, one goes from a regime dominated by the black hole to a situation where the scalar field contributes more to the total mass. Last, as an illustration, some profiles of |Φ| = R Φ RΦ + I Φ ĪΦ 1/2 are shown in Fig. 14, along the x-axis. One can note that the maximum value of the field is located at about 40 times the radius of the black hole r H = 1. As expected, the further Ω BH is from the value of the Kerr black hole used to initiate the sequence, the higher the amplitude of the scalar field.

VI. LAST WORDS

A formalism to compute spacetimes containing stationary black holes is presented. The geometry is described by the 3+1 decomposition of spacetime and the various metric fields are found using the maximal slicing gauge for the choice of time and the spatial harmonic gauge for the spatial coordinates. The presence of the hole itself is enforced by demanding that a given sphere is an apparent horizon. Moreover, as stationarity is assumed, one demands that this horizon is in equilibrium. Using a combination of analytical and numerical studies, a set of boundary conditions for the various metric quantities is found. Let us point out that this leads to a choice of coordinates that is regular everywhere, even on the horizon itself. This is to be contrasted with the analytic Boyer-Lindquist coordinates of the Kerr black hole, which are singular on the horizon. This property may prove useful in numerical applications.

The whole procedure is applied to three different situations. The system of equations is solved numerically using spectral methods and the Kadath library [START_REF] Grandclement | Kadath: A Spectral solver for theoretical physics[END_REF][START_REF] Grandclement | Kadath libray[END_REF]. The first and most simple application is to recover the classical Kerr black hole (in coordinates that are not analytical). Configurations up to a Kerr parameter a/M of 0.99 are easily computed.

The second application concerns a static black hole where gravity is minimally coupled to a real scalar field, with a negative cosmological constant. It is called the MTZ black hole, from the names of the authors who obtained it analytically [START_REF] Martínez | Exact black hole solution with a minimally coupled scalar field[END_REF]. The formalism appears to work well also in the case of this spacetime which is asymptotically anti-de Sitter. The last application is devoted to the computation of a family of black holes with complex scalar hairs [START_REF] Herdeiro | Kerr black holes with scalar hair[END_REF][START_REF] Herdeiro | Construction and physical properties of Kerr black holes with scalar hair[END_REF]. This last example combines a coupling with matter (i.e. the scalar field) and the inclusion of rotation. For the three cases, errors are carefully checked, especially by showing fast convergence to zero of several error indicators when the numerical resolution increases.

In the future, it is hoped that the formalism presented here will be a valuable tool to study various models of black holes. A lot of the results in the field of black holes in alternative theories of gravity relies, at some degree, on the analytic choice of coordinate systems. This can be a difficulty to get the more general solutions possible. This is especially true when rotation is included. For instance, there are no rotating equivalent of the MTZ black hole known yet.

Another extension of this work concerns black holes that are not in exact equilibrium. In particular, this is the case for objects in binary systems. The boundary conditions proposed here could be applied to generate improved initial data for binary coalescence simulations.

Usually initial data are generated assuming that the spatial metric is conformally flat (see [START_REF] Gourgoulhon | Binary black holes in circular orbits. 1. A Global space-time approach[END_REF][START_REF] Grandclement | Binary black holes in circular orbits. 2. Numerical methods and first results[END_REF][START_REF] Caudill | Circular orbits and spin in black-hole initial data[END_REF][START_REF] Papenfort | New public code for initial data of unequal-mass, spinning compact-object binaries[END_REF]). This assumption could be relaxed using the techniques developed in this paper.

However, due to the emission of gravitational waves, the horizons are not in equilibrium and thus it is not expected that the boundary conditions will be exact. However their accuracy should improve as the separation increases. Equilibrium of the horizons are not the only place where deviation from stationarity must be accounted for. The terms in ∂ t present in the equations (3-5) must also be considered. Outer boundary conditions must also be carefully investigated, as simple spatial asymptotic flatness may not be sufficiently precise when gravitational waves are present. It is beyond the scope of this paper to implement the boundary conditions presented here in the binary black hole context. However this is something that is planned for the future and one can hope that the formalism exposed here will lead to significant improvement in the precision of the computation of initial data for binary black holes configurations.

  FIG. 1. For a non-rotating black hole, maximum value of the expansion |Θ| on the horizon, as a function of the resolution (i.e. number of radial collocation points). Three different values of γ rr on the horizon are displayed. Θ clearly goes to zero when precision increases.

FIG. 2 .

 2 FIG. 2. Various error indicators for a configuration with r H = 1 and Ω BH = 0.1. All the quantities are shown as a function of the number of spectral coefficients in both the r and θ directions. The first panel shows the gauge quantities K and V i , the second one the maximum value of Θ on the horizon, the third one the relative difference between M ADM and M Komar and the last one the relative difference between the numerical and analytical values of Ω BH .

FIG. 3 .FIG. 4 .

 34 FIG. 3. The first panel shows Ω BH as a function of a/M for a sequence of constant radius r H = 1 and the second M Ω BH as a function of a/M also. Configurations are computed with the highest resolution at hand.

FIG. 5 .

 5 FIG. 5. Maximum value of the gauge quantities K and V r and of the expansion on the horizon, as a function of the radial number of coefficients. Spectral convergence is clearly seen.

FIG. 6 .

 6 FIG. 6. Profiles of various fields (first panel N , second B r , third γ rr and γ θθ and fourth φ) for r H = 8 and 21 radial coefficients in each domains. As expected, the lapse and spatial metric diverge near the ADS border.

FIG. 7 .

 7 FIG. 7. Four-dimensional Ricci scalar on the horizon as a function of the value of the scalar field on the horizon. The solid line is computed from the analytic solution, while the circles correspond to the numerical results. Horizon radii r H range from 6.5 to 10 (from right to left). Large radii correspond to reduced scalar field and weaker curvature.

FIG. 8 .

 8 FIG. 8.

FIG. 9 .

 9 FIG. 9. Maximal values of K and V i (first panel) and of Θ on the horizon (second panel) as a function of resolution. The configuration emerges from a Kerr black hole with M ≈ 5 and a/M ≈ 0.95 and corresponds to Ω BH ≈ 0.069831. The various indicators exhibit a spectral convergence to zero.

  numbers of points. The chosen configuration corresponds to Ω BH ≈ 0.069831. The maximum values of K and V i in the whole space and of Θ on the horizon are shown in Fig. 9. The three different resolutions correspond to N r = N θ = 13, 17 and 21, where N r and N θ denote the number of collocation points in the r and θ directions. The two panels of Fig. 9 clearly show a spectral convergence of the various quantities proving the validity of the solution.

FIG. 10 .

 10 FIG. 10. Profiles along the x-axis of R Φ (blue curves) and I Φ (red curves) for the zero phase condition (solid lines) and the half phase one (dashed lines) (see text for details). The second panel shows the region close to the origin and the vertical line indicates the surface where the value of the phase is enforced.

FIG. 11 .

 11 FIG. 11. Discontinuity of the radial derivative of I Φ as a function of resolution (i.e. the number of radial and angular coefficients). The configuration is the same as the one of Fig. 9.

FIG. 12 .

 12 FIG. 12. ADM mass (first panel) and angular momentum (second panel), as a function of Ω BH , for a sequence emerging from a Kerr black hole with M ≈ 5 and a/M ≈ 0.95. The Kerr black hole corresponds to the point at the right end of the curves.

FIG. 13 .

 13 FIG. 13. The first panel shows the quantity q, as a function of Ω BH . The second panel shows the ratios M BH /M ADM and M Φ /M ADM , also as a function of Ω BH . Both panels illustrate that one goes from a pure Kerr black hole situation (right side) to a situation where the field is dominant (left side).

FIG. 14 .

 14 FIG. 14. Values of |Φ| along the x-axis, for four different values of Ω BH . The black hole horizon is located at r h = 1.

  points in r and θ), one can compute a sequence of Kerr black holes, with different values of the Kerr parameter a/M . It is well known that this parameter goes from 0 (Schwarzschild black hole) to 1 (extremal Kerr black hole). In order to compute the sequence, one can start from the non-rotating solution and increase Ω BH step-by-step. If this does work properly at first (i.e. for moderate values of a/M ), a

ACKNOWLEDGMENTS

The authors would like to thank Eric Gourgoulhon for many fruitful scientific discussions on the topics covered by this paper. This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche.

Appendix A: System of equations for the MTZ black holes 1. Matter terms for the MTZ black hole For a stationary real scalar field, the various 3+1 matter terms can be obtained as:

Regularization at the ADS border

For the MTZ black hole, the regularized quantities used close to the ADS boundary are given. The various quantities are multiplied by the appropriate power of Ω (Eq. ( 36)) in order to avoid divergences. It leads to

The regularized equations are thus:

Appendix B: Matter terms for the Einstein-Klein-Gordon system

The 3+1 matter terms derived from the ansatz (45) are computed. Given the dependence in terms of t and ϕ, one can compute the various following terms

Even if R Φ and I Φ are real quantities, the notations RΦ and ĪΦ are kept to differentiate terms that are in factor of exp [i (ωt -kϕ)] and those of exp [-i (ωt -kϕ)].

The terms (B1-B3) enter into the expressions of the components of the stress-energy tensor as

From those components, the 3+1 matter terms are expressed as

)