{'fraction_non_alphanumeric': 0.03020265003897116, 'fraction_numerical': 0.006966095089633671, 'mean_word_length': 4.612083105522143, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stuffed Black Holes

Initial data corresponding to spacetimes containing black holes are considered in the time symmetric case. The solutions are obtained by matching across the apparent horizon different, conformally flat, spatial metrics. The exterior metric is the vacuum solution obtained by the well known conformal imaging method. The interior metric for every black hole is regular everywhere and corresponds to a positive energy density. The resulting matched solutions cover then the whole initial (Cauchy) hypersurface, without any singularity, and can be useful for numerical applications. The simpler cases of one black hole (Schwarzschild data) or two identical black holes (Misner data) are explicitly solved. A procedure for extending this construction to the multiple black hole case is also given, and it is shown to work for all time symmetric vacuum solutions obtained by the conformal imaging method. The numerical evolution of one such 'stuffed' black hole is compared with that of a pure vacuum or 'plain' black hole in the spherically symmetric case.


I. INTRODUCTION AND OVERVIEW
Black holes are the most elementary objects in General Relativity. Allowing for the well known 'no hair' theorems, they can be characterized by their mass, charge and spin, like elementary particles in classical or quantum mechanics. This is because event horizons act as a one way membrane which separates the black hole interior from the exterior: there is no causal effect that can propagate from the inner to the outer region. Therefore, one does not need to know all the details about physical processes taking place at the interior in order to describe the overall evolution of the black hole. This also means that black hole metrics that differ only inside the holes can lead to the same exterior spacetime. The idea of this paper is to replace the singular vacuum solution for a black hole interior by a regular one corresponding to a non-vanishing energy density, but keeping the same exterior metric. The term 'stuffed' black hole is a good description of the resulting solution.
The 'no hair' property has been very useful in Numerical Relativity to avoid the interior black hole singularity appearing in the initial data. The standard practice is to excise the inner region from the computational domain, so that one can safely evolve the exterior (usually vacuum) region. This has been done for a single black hole (Schwarzschild case) by setting an internal boundary on the computational domain at the initial position r 0 of the apparent horizon. The same thing has been done for a system consisting of two non rotating black holes, starting the evolution from the Misner data [1,2], or even for an arbitrary number of them, by using the 'conformal imaging formalism' [3] to construct conformally flat initial data. We shall refer to all these pure vacuum solutions as 'plain' black hole metrics to distinguish them from the 'stuffed' ones we are presenting here.
In this paper we propose to take advantage of the 'no hair' property in a different way. We will keep the same exterior metric of the 'plain' black hole case, obtained by the conformal imaging formalism. But we will match this exterior metric to a regular conformally flat interior corresponding to a positive energy density. The exterior solution is not affected because the matching surfaces coincide in every case with the initial position r 0 of the black hole apparent horizon. This is done for a single non rotating black hole (Schwarzschild metric) in Section II. The interior metric in this case is the spatial part of a positive curvature Friedmann-Robertson-Walker (FRW) metric, so that the resulting construction can be interpreted as initial data for the Oppenheimer-Schneider dust collapse (Section III). In Section IV, we give in closed form the corresponding solution for two nonrotating stuffed black holes (the analogous of the Misner data). We provide in Section V that a step by step procedure to stuff time symmetric initial data containing any number of black holes with arbitrary sizes and locations.
Finally, in the Appendix, we compare the numerical evolution of a 'plain' black hole with that of a 'stuffed' black hole in the spherically symmetric case. This is a biased test, because spherical coordinates are singular at the origin and this coordinate singularity affects only the stuffed hole, because in that case there is no excised region and the computational domain includes the origin. In spite of that, we have evolved stuffed holes with the same accuracy and stability as plain holes. This opens the door to three-dimensional applications, where the 'stuffing' approach avoids putting internal curvilinear boundaries that would need a special treatment. The advantage of the 'stuffed' black holes versus the 'plain' ones will be manifest if the coordinate system contains a shift vector which allows the holes to move across the numerical grid: in the plain case one would need to treat a number of moving curvilinear boundaries, but in the stuffed case nothing special is to be done because there are no internal boundaries in that case.

II. SINGLE HOLE INITIAL DATA
The evolution formalism [4][5][6] is specially well suited for Numerical Relativity. In normal coordinates, the spacetime metric can be written as To start the evolution one needs both the spatial metric γ ij on the initial hypersurface and its second fundamental form (extrinsic curvature) These initial data, however, are constrained by the following equations (energy constraint), and (momentum constraint), where both the Ricci tensor R ij and the covariant derivatives are the ones corresponding to the three-dimensional geometry of the initial slice. These constraint equations are first integrals of the Einstein equations: they are verified at every time slice by the spatial part of any spacetime metric. One can learn a lot about the spacetime by just looking at the geometry of the spatial hypersurfaces. Let us consider for instance a closed two-surface σ in one such a threedimensional manifold. The expansion θ of a congruence of outgoing light rays starting at σ is given by [7] where the three-dimensional vector n i is the unit normal to σ. In the 'time symmetric' case (K ij = 0) the minimal surfaces of the three-geometry are also apparent horizons (θ = 0) and vice versa.
In the case of non-rotating black holes, one usually starts with time symmetric initial data (and therefore zero initial momentum density S i ) and conformally flat initial metrics, namely This simplifies the constraint equations [8], so that (3) can be written in terms of a flat space Laplace operator ∆ In the spherically symmetric case, the regular vacuum solution of (8) can be written as which corresponds to the space part of the Schwarzschild metric in isotropic coordinates. Allowing for (7), the minimal surface (apparent horizon) condition (6) reads which holds for r 0 = m/2. Note that the combination r ψ 2 is precisely the 'area radius' R of the spherical geometry, so that the apparent horizon (minimal surface) is placed at Condition (10) can also be obtained by realizing that the metric given by (7,9) is invariant under the discrete symmetry which can be interpreted as a coordinate inversion at r = r 0 . This inversion symmetry provides suitable boundary conditions when one excises the spherical region r < r 0 from the computational domain. The resulting inner boundary condition is very easy to implement numerically in spherical coordinates, but it is much more difficult to manage in the threedimensional case (where one would get a curvilinear boundary in cartesian coordinates) even for a single black hole.
As an alternative to this inner boundary aproach, we will construct a 'stuffed' black hole by replacing the initial data (9) with for r > m/2 so that in the exterior region one recovers the 'plain' black hole solution (9), whereas in the interior region one gets, allowing for (7), a homogeneous three-dimensional metric which is the space part of a closed Friedman-Robertson-Walker (FRW) model (see Fig. 1).

III. PHYSICAL INTERPRETATION
Up to now, we have just constructed time symmetric data on the initial slice. This means that we have used only the energy and momentum constraints, but not the remaining (evolution) field equations. From the energy constraint (8), one can compute the energy density whereas the momentum constraint (4) plus the time-symmetry condition (K ij = 0) imply that the momentum density should vanish. However, in order to get a physical interpretation of the matter content of the spacetimes generated by these initial data, we need to say something about the remaining (space) components of the stress-energy tensor T µν .
To do this, we can study the time evolution of the matching conditions between the constant density interior and the vacuum exterior, namely where Φ = 0 is the equation of the matching hypersurface. Allowing for the fact that the momentum density vanishes, one gets easily from (15) that Φ can not depend on time. This means that the matching surface can be taken to be the sphere for any value of time. It follows then from (15) that the radial direction gives an eigenvector of the stress-energy tensor (computed at the matching surface) with zero eigenvalue. If we assume for simplicity a barotropic perfect fluid matter content, it follows that the pressure should vanish (incoherent matter). The initial data (13) can then be understood as corresponding to a particular case of the well known Oppenheimer-Schneider dust collapse: a constant density spherical star which is initially at rest. In our case, the initial star radius actually coincides with the position of the apparent horizon.
The former is just one of the many possibilities allowed by the no-hair theorems. A very interesting case is the 'string perfect dust' matter content [9] T µν = S µρ S ν ρ , where is a simple surface-forming bivector.
There are basically two different possibilities to get a zero radial eigenvalue in the string case [9]: either S µr = 0 (the bivector has no radial component) or the vector λ µ = S µr is isotropic (and, being antisymmetric, has no radial component). In both ways, it follows that the string stress tensor breaks the spherical symmetry of the initial data, and this will lead to a non-spherical time evolution. This 'dynamical symmetry breaking' is in contrast with the more predictable behaviour that one gets in the dust case, where the full stress-energy tensor is spherically symmetric, so that the spherical symmetry of the initial data will be preserved during evolution.

IV. TWO BLACK HOLES INITIAL DATA
The Misner initial data [1] are axially symmetric and describe two identical non-rotating black holes which are initially at rest. They can be obtained by linear superposition of spherically symmetric solutions of the Laplace equation, with centers distributed along the symmetry axis: where In order to see how this solution is obtained, let us notice that in the time symmetric vacuum case the constraint (8) on the conformal factor ψ is the flat space Laplace equation. We know from electrostatics that we can take advantage of the invariance of the Laplace equation under discrete symmetries, such as (12), by using the 'imaging method' [10]. This technique was adapted to Relativity by Misner [12] in order to obtain initial data which are invariant by inversion across a number of spherical surfaces, which will become at the end minimal surfaces (apparent horizons) of the resulting solution. This solution can be interpreted as describing time symmetric initial data for a number of black holes.
The apparent horizons (minimal surfaces) in (19) are the two spheres σ ± given by The term Φ + n+1 has been constructed as the image of Φ − n under inversion across σ + . Also, the term Φ − n+1 has been constructed as the image of Φ + n under inversion across σ − . This means that the dipole combinations are invariant under inversion across the sphere σ ± , respectively. We have seen in the previous section how to stuff a single black hole. Now we have two identical holes and we will look for regular interior solutions for each one. Let us begin by considering the interior region to the first sphere σ − . The infinite sum in the Misner solution (19) can here be written in a more convenient form so that every term is now invariant under inversion across the sphere σ − . We will make use of this symmetry property to match every term in (24) to (the conformal factor of) a constant curvature metric.
The conformal factor for a constant (positive) curvature metric (closed FRW) can be written as where F is an arbitrary scale factor and λ and b are parameters related to the conformal transformations of Euclidean three-space. If we impose inversion symmetry (Eq. 10) across the sphere σ − , we get This single condition ensures, allowing for (10), that the normal derivatives of the vacuum and FRW conformal factors will coincide on σ − if and only if both conformal factors actually coincide there. It follows that, in order to complete the matching, we need only to look at the conformal factor values at σ − and tune the free parameters in (25). The FRW conformal factor can be easily evaluated on σ − : where we have noted On the other hand, a straightforward calculation shows that Now it is easy to obtain the values of the arbitrary parameters in (25) that ensure the matching between the vacuum dipole term Λ − n and the conformal factor across σ − for every value of n. Their actual values are: Allowing for the symmetry of the solution across the equatorial plane, it is easy to see that the conformal factor ψ + n = F n λ n 1 + 2 b n (z + a coth µ 0 ) + sinh 2 µ 0 (r + 1 /a) 2 will also match the vacuum dipole term Λ + n across σ + . It follows that the complete solution for the stuffed two hole problem reads The resulting solution is regular and smooth everywhere (see Fig. 2). Following Bowen et al [11], we can define an effective energy density starting from a 'Newtonian potential' linearly related to Ψ. In our case, allowing for (8), this amounts to This effective energy density can be easily expressed as the sum of the effective energy density of every FRW component, namely The positivity of the factors ψ ± n ensures that the energy density (37) is positive inside the holes and it is bounded by the maximum energy density of the FRW components.
As in the single hole case, the energy density has a jump at the matching surfaces due to the discontinuity of the second metric derivatives there. It follows from the previous section that a reasonable physical description of the solution could be that of two balls of incoherent matter. Notice, however, that the energy density of every ball is not constant, as it was in the single hole case. The gravitational interaction between the two balls accounts for their distortion.

V. MULTIPLE BLACK HOLE CASE
Let us consider now the time symmetric multiple black hole case. As is well known, the vacuum exterior solution can be obtained by the conformal imaging method [12]. As an input for this method, one must provide the size and location of an arbitrary number N of spheres, which will become at the end the apparent horizons of N black holes. The resulting solution, by construction, will be then inversion symmetric across the apparent horizon of every black hole. It could be written as a linear superposition of poles, Our goal is to provide a suitable interior solution for every black hole. As far as the holes do not overlap, we can consider them separately. Let us begin with the first one: its horizon σ (1) is a sphere centered at x 1 , Note that, as (39) is invariant under inversion across the sphere σ (1) , half of the poles in (39) are outside and half inside. Moreover, every outside pole is the image under inversion of an inside pole and vice versa. This means that we can combine every pole with its image to form the invariant dipoles where the exterior point x j ′ is the image of the interior one x j under inversion across the sphere σ (1) . The sum in the solution (39) can then be rearranged as follows Our strategy will be then to match separately every invariant dipole (41) to a closed FRW factor, given by As in the previous section, we must impose first that the interior solution (43) have a 'throat' (10) at the sphere σ (1) , obtaining the condition which generalizes (26). As in the previous section, this condition ensures that the radial derivatives of the vacuum and the FRW conformal factors will coincide on the sphere σ (1) if and only if both conformal factors actually coincide there. Alowing for (45), the interior factor (43) at the horizon is On the other hand, allowing for the fact that the two images in the invariant dipole (41) coincide by construction on the sphere σ (1) , we have It follows that the choice of parameters ensures the matching between (41) and (43) at σ (1) . The interior solution for the first hole can then be obtained by summing the corresponding FRW factors (43) The interior solution for every other hole can be obtained in exactly the same way.

APPENDIX: 1D NUMERICAL EVOLUTION
We will compare here the numerical evolution of a single spherically symmetric (1D) stuffed black hole against a 'plain' one. We use in both cases the same finite difference code with a 400 point evenly space numerical grid. The metric quantities and their evolution equations are also the same. The matter density and speed are computed in the stuffed case by using the standard 'upwind' method to model the continuity equation and the Euler equation for dust, respectively (they are equivalent to the stress-energy tensor conservation).
In the 'plain' case, we have used the throat inversion symmetry to provide the inner boundary condition. In the general (3D) stuffed case there is no inner boundary. The only concern will arise from the discontinuity of the energy density (second derivatives of the metric) at the throat. We have not seen any problem with the matter variables (density and speed) in our 1D numerical evolution.
The singularity of spherical coordinates at r = 0, however, demands a special treatment of the origin in the 1D stuffed case. We have done it by adding two virtual points as the mirror image of the first two grid points across the origin. In the 'plain' hole case the origin is contained in the excised region and such special treatment is not needed. This is why this 1D test is biased: it is the best case for the plain hole (inversion symmetry in spherical coordinates) and the worst one for the stuffed case (need for a special treatment of the origin).
We have compared the code performance, in terms of accuracy and stability, in both cases for three different slicing conditions: harmonic, '1+log' and maximal. In the last two cases, we can barely notice any significant difference between the plain and the stuffed hole evolution in the exterior region (see Fig. 3). In the harmonic slicing case, however, the code for the stuffed black hole crashes as the dust ball in the interior region collapses, evolving towards a singularity. This is a consequence of the weaker singularity avoidance properties of harmonic slicing.
These results are promising with regards to 3D applications in rectilinear grids, where inversion symmetry across the spherical horizons is not so easy to implement (plain case), whereas there is no internal boundary of any kind (the coordinate system is regular everywhere) in the stuffed case. Plot of the conformal factor describing initial data for a single stuffed black hole with mass m = 2 (apparent horizon at r = 1 in isotropic coordinates). The smoothness and regularity of the solution is evident.
FIG. 2. Surface plot of the conformal factor describing initial data for two stuffed black holes with µ 0 = 2, a = 1. The plot is in cylindrical coordinates with the axis along the line joining the two holes. The azimuthal angle φ is suppressed, as the solution is axially symmetric. The smoothness and regularity of the solution is evident.
FIG. 3. Evolution of plain (solid line) and stuffed (dashed line) black holes in the spherically symmetric case ('1+log' slicing, 400 grid points). Left: the radial metric coefficient is plotted after t = 300m. Right: the maximum error for the black hole mass near the horizon is plotted as a function of time. Both the stuffed and the plain case show the same accuracy.