Multicenter solutions in Eddington-inspired Born-Infeld gravity

We find multicenter (Majumdar-Papapetrou type) solutions of Eddington-inspired Born-Infeld gravity coupled to electromagnetic fields governed by a Born-Infeld-like Lagrangian. We construct the general solution for an arbitrary number of centers in equilibrium and then discuss the properties of their one-particle configurations, including the existence of bounces and the regularity (geodesic completeness) of these spacetimes. Our method can be used to construct multicenter solutions in other theories of gravity.


I. INTRODUCTION
Among the many exact solutions of interest known within General Relativity (GR), we find the class of gravitating configurations in self-equilibrium. Such is the case, for instance, of geons [1], gravitating solitons [2,3] and skyrmions [4,5], and other long-lived configurations [6]. There is a sub-class of solutions of such a family corresponding to those where the attractive gravitational field of a system of particles is exactly balanced by a repulsive electrostatic force. Revolving around some previous considerations by Weyl restricted to axial symmetry [7], the first explicit solution of this kind, corresponding to the sourceless Einstein-Maxwell field equations, was first found by Majumdar [8] and, independently, by Papapetrou [9]. It was latter shown by Hartle and Hawking [10] that the Majumdar-Papapetrou solution can actually be interpreted as a collection of extremal black holes in static equilibrium [11], raising further interest on the topic from different communities [12][13][14][15]. A particularly interesting property of this solution is its lack of symmetry in the distribution of black holes. Many other solutions of this family and generalizations have been subsequently found and characterized in the literature regarding their astrophysical features [16][17][18][19][20][21][22].
In the quest for counterparts of gravitating configurations in extensions of GR a number of solid scores have been hit, for instance, within rotating black holes [23], compact stars [24], further horizonless compact objects [25], and so on. The search for such configurations is suggested from the point of view of finding alternatives to canonical GR objects (either black holes or compact stars) whose properties can be tested against the present and future stream of data from multimessenger astronomy [26], and act as observational discriminators with respect to GR predictions. In implementing this program one faces the fundamental difficulty of the inherently more involved equations of motion of most extensions of GR, which largely prevents the construction of solutions of theoretical and observational interest (among other troubles). It is therefore timely and of great relevance the development of novel methods and algorithms to shortcut the structure of such field equations to find explicit solutions.
One such method has been developed recently [27]. It works for theories of gravity including scalar objects built out of contractions of the (symmetric part of the) Ricci tensor with the metric, and formulated in metricaffine (Palatini) spaces, where metric and affine connection are regarded as independent entities. The resulting family of models are the so-called Ricci-Based Gravities (RBGs), which yield second-order, ghost-free equations of motion which are compatible with all solar system and gravitational wave observations so far. For these theories it is possible to introduce new variables such that the corresponding field equations are cast in the Einstein frame, in such a way that the nonlinearities are transferred to the matter sector. From this frame one can take (when known) the corresponding GR solution, and find the counterpart in the RBG frame via purely algebraic transformations. This is known as the mapping method, whose reliability was originally proven for anisotropic fluids [27], and has allowed to obtain new solutions for electromagnetic [28] and scalar fields [29,30].
The main aim of this paper is to progress further in the analysis of the capabilities of this mapping by working out the counterpart of the Majumdar-Papapetrou (MP) solution using one particular RBG, namely, the so-called Eddington-inspired Born-Infeld (EiBI) gravity theory, popularized by Ferreira and Banados [31], though first studied in the metric-affine formalism by Vollick [32] inspired on a work by Deser and Gibbons [33]. This choice is motivated on the plethora of applications of this theory within astrophysics, black holes physics and cosmology, as found in the last few years (see [34] for a recent review). We shall explicitly cast the mapping method for this theory and the combination of fluid and electromagnetic fields supporting the MP solution, and then we will build its generalization within EiBI gravity. We shall show that, depending on the sign of the EiBI parameter, the corresponding solutions can be understood as i) a collection of point-like objects with bounded stress-energy density everywhere and that lose their (extremal) horizon in the lower part of the mass spectrum, and ii) as a set of extremal black holes in equilibrium, each of which contains a non-traversable wormhole (or black bounce [35][36][37]) in its interior. We discuss the regularity properties of both families of solutions regarding the completeness of geodesics and the behavior of curvature scalars.
The paper is organized as follows: in Sec. II we introduce the main elements of the MP solution in GR. The RBG family of theories and the mapping with electromagnetic fields is presented in Sec. III. We then combine these two ingredients to present the counterpart of the MP solution in Sec. IV, discussing its horizon structure and geodesic equation for one-particle configurations. We finally conclude in Sec. V with a discussion and some perspectives.

II. EXTREME COUNTERPOISED DUST AND THE MAJUMDAR-PAPAPETROU SOLUTION
The Majumdar-Papapetrou family is a particular class of solutions solving Einstein equations where κ 2 = 8πG is Newton's constant, while the contributions to the energy-momentum tensors read which corresponds to a pressureless dust component, where ρ is the energy density and the unit time-like vector u µ u µ = −1, and which corresponds to a standard (Maxwell) energymomentum tensor associated to the electromagnetic field, where F µν = ∂ µ A ν − ∂ ν A µ is the field strength tensor of the vector potential A µ . The Einstein field equations (1) must be complemented with the matter field equations, which read where the charge density J ν = ρ e u ν . This is nothing but Maxwell field equations sourced by a current generated by the pressureless fluid (2).
Assuming purely electric fields, A µ = δ t µ φ, Weyl showed [7] that the most general relation between the metric and the electric potential solving both the Einstein and electromagnetic equations must be of the form (A, B some constants) supported by axially symmetric spatial symmetry. The MP solution generalizes the proposal of Weyl to any spatial symmetry by imposing the following more stringent restriction on the relation between the metric and the electrostatic potential Enforcing this condition, it can be proved that the solution is static with the density of the dust that equals the electric charge distribution, namely ρ e = ρ, which is known as an extreme counterpoised dust. Moreover, the spatial sections are conformally flat when the pressure is set to zero [38]. Under these conditions, the MP background spacetime metric (in Cartesian coordinates) can be expressed as where U is a function characterizing the geometry that is related to the electrostatic potential through the MP condition (6). Since the Einstein-Maxwell field equations only depend on derivatives of φ, this potential can be redefined to cancel the constant C so that (6) simplifies as This is the choice of [10] up to a choice of units of the electric charge that can reabsorb the √ 2 factor. On the other hand, from the line element (7) the matter field equations (4) result in the following Poisson equation that reduces to Laplace equation for the electrovacuum solution, recovering the original MP solution.
In summary, MP solutions are a particular subset of the Einstein-Maxwell-dust system in which the mass is exactly tuned to the electric charge. For instance, any collection of extreme Reissner-Nordström black hole solutions (m 2 = q 2 ) located at will is a particular MP solution [11], without any need to impose additional symmetries. For this reason, these configurations are sometimes called multicenter solutions. In any case, of particular interest are those configurations enjoying spherical symmetry, in which case Eq.(7) can be expressed as and thus Eq.(9) remains formally the same but now the Lagrangian and its functional dependencies are expressed in terms of R. A further coordinate change may bring the line element into standard Schwarzschild form, that is ds 2 = −A(r)dt 2 + B(r)dr 2 + r 2 (dθ 2 + sin 2 θdφ 2 ) , (11) via the identifications r = RU (R), A(t) = U −2 (R) and B(r) −1/2 = 1 + R U dU (R)/dR. Either under form (10) or (11), there are typically two paths followed in the literature to solve the corresponding field equations: i) one assumes a functional form for U (R) and solves (9) to find the matter energy density ρ(R) threading the geometry [39,40] or ii) a function ρ(R) for the inner region is set [41], and resort to numerical methods to resolve the corresponding differential equation in order to get U (R). Perhaps the most well known example of the first path is the so-called Bonnor stars [42], where one sets two different regions as where the exterior solution, U E corresponds to an extreme Reissner-Norström black hole, which is matched to the interior solution at a certain r = r 0 .

A. Ricci-based gravities
To describe the mapping procedure for electromagnetic fields and the main elements required for our analysis, let us begin by defining the action of RBGs as where g is the determinant of the spacetime metric g µν . The functional dependence of the gravitational Lagrangian L G must be through traces of powers of the object M µ ν ≡ g µα R (αν) , where R (αν) (Γ) is the symmetric part of the Ricci tensor 1 , which is solely built out of the affine connection Γ ≡ Γ λ µν (not necessarily symmetric [45]), the latter being a priori independent of the metric (Palatini approach). Regarding the matter Lagrangian L m (g µν , ψ m ), it depends on a set of matter fields ψ m minimally coupled to the spacetime metric.
The field equations for RBGs can be conveniently written as [27,44] where G µ ν (q) is the Einstein's tensor of a new rank-two tensor related to the spacetime metric via the fundamental relation Here Ω α ν is dubbed as the deformation matrix, and can always be written on-shell as a function of the stressenergy tensor T µν (g) ≡ 2 √ −g δLm δg µν . The relation between T µ ν (g) andT µ ν (q) follows from the original RBG field equations and takes the form where vertical bars denote a determinant and T ≡ g µ T µν is the trace of the stress-energy tensor. This effective stress-energy tensorT µ ν (q) can be similarly derived from another Lagrangian densityL m , that is, This procedure establishes a correspondence or mapping between RBGs coupled to L m and GR coupled toL m , as established in [27][28][29][30]. Note that, in general, both L m andL m will contain fields of the same kind (that is, scalar fields map into scalar fields, electromagnetic into electromagnetic, and so on), though the functional dependence will be different, yielding in general non-canonical Lagrangians in one (or both) sides.

B. The mapping for Eddington-inspired Born-Infeld gravity
In this work we shall be interested on the case where we choose the following RBG Lagrangian density where ǫ is a constant with dimensions of length squared, and the theory features an effective cosmological constant Λ ef f = λ−1 ǫ . This is the well known and well studied Eddington-inspired Born-Infeld (EiBI) gravity [31,34]. We shall first perform a few manipulations in order to apply the mapping above (i.e. Eq. (16)) to this setting. The equations of motion for EiBI gravity are obtained by variation of the action (17) with respect to the connection and the metric as the two sets of nonlinear equations (here S ν αγ represents the torsion tensor) respectively. Following [45] we can introduce an auxiliary metric q µν defined by the equation and after performing a projective transformation 2 it is possible to reduce the field equations associated to the variation of the connection to the metric compatibility condition,∇ ρ q µν = 0 3 . Eq. (19) applied to (17) reveals that the auxiliary metric appearing in Eq. (15) is which is a well known result since the seminal paper [31]. Using the following general recipe of [48] to generate the metric for the EiBI gravity (17) coupled to any matter theory it is straightforward to find the following relation between the metrics in the case of any nonlinear electrodynamics where we have introduced the two invariants of the electromagnetic field the expression above boils down to This equation provides a direct shortcut to find any solution on the EiBI side (as given by g µν ) starting from a seed solution on the GR side (as given by q µν ). In the next section we shall use this powerful result in order to generate the counterpart of the MP solution within EiBI gravity coupled to BI electrodynamics. Let us now focus on electromagnetic fields. Following the procedure detailed in [48,49], the field equations (18) can be reduced to the standard Einstein equations written in terms of the auxiliary metric and the tilted connection if and only if the matter sector in the EiBI frame is related to the matter sector in the GR frame through the following parametrization [48] provided that we choose the matter content on the GR side as given by Maxwell electrodynamics in (25). Our next goal is to express (28) in terms of quantities in the EiBI frame by writing the invariants of the q µν frame, which are those appearing in the Lagrangian density (28), in terms of those of the RBG frame (the untilted variables). Using the inverse mapping between metrics (21), a little algebra allows to find the following relations between the field invariants in the GR and RBG frames [48]:K Replacing these expressions into the parametrization (28) one gets the form of the matter Lagrangian in the EiBI frame as Now, if we take from now on, for simplicity, asymptotically flat solutions, λ = 1, and make the identification β 2 = 4π/(ǫκ 2 ) choosing the minus sign solution, then the Lagrangian density above becomes which is the well known Born-Infeld (BI) theory of electrodynamics. Therefore, we have just seen that GR coupled to Maxwell electrodynamics maps into EiBI gravity coupled to BI electrodynamics. Note, however, that the sign of ǫ is not a priori restricted to be positive, which implies that β 2 could also be negative despite being written as the square of a parameter β.

A. The Majumdar-Papapetrou solution in EiBI
Now that we have all the necessary elements of the mapping for this scenario under control, the derivation is remarkably simple. Indeed, in order to extract the counterpart of the MP solution (7) within EiBI gravity we just need to compute the extra corrections appearing in Eq.(26) via the ansatz (10). First we find that where {i, j} = x, y, z, which allows to find so that after a bit of algebra Eq.(26) reads which is the counterpart of the MP solution, given by a certain U ( x), within EiBI gravity (17) coupled to Born-Infeld electrodynamics (32), as obtained via the mapping. Nice and easy.
To construct multicenter solutions out of the above solution one can consider a collection of N solutions of this type such that the metric function (in GR) is written as where m i labels the mass of each object, and ( where (x i , y i , z i ) are the coordinates labeling the center of each solution. If instead of Cartesian coordinates one labels the centers by means of spherical coordinates, each vector x i becomes x i = r i (sin θ cos ϕ, sin θ sin ϕ, cos θ) such that For the sake of the discussion of the features of the generalized MP solutions, it is much more convenient to rewrite the generalized MP solution (35) in terms of the electrostatic potential. This can be done after noting that the expression of the electrostatic potential A µ (x, y, z) = (φ(x, y, z), 0) with φ(x, y, z) = 1/U (x, y, z) allows to write ∇U/U 2 = − ∇φ, which simplifies many expressions. Replacing this into the line element (35), and suitably rearranging terms one finds where dΩ 2 = dθ 2 + r 2 sin 2 θdϕ 2 is the usual unit volume element of the two-spheres. In order to understand the properties of this multicenter solution, it is useful to have a careful look at individual centers, which will provide a reasonable approximation of the geometry for sufficiently separated objects. Corrections depending on the separation could be computed by perturbative methods and the general case should take into account the exact line element (38).
In the case of having only one center, the angular contributions to the anisotropy vanish and only the radial part remains, allowing us to write ∇φ · d r = φ r dr. The line element (38) can then be written as Given that for one-particle configurations one has the electrostatic potential φ = r/(r+m), where m is its mass, it is easy to go from the current isotropic coordinates to Schwarzschild-like ones by redefining the factor r 2 /φ 2 = R 2 = (r + m) 2 , which turns the above one-particle line element into where now φ(R) = 1 − m/R, the parameter s = ±1 denotes the sign of ǫ, and R 4 c ≡ |ǫ|κ 2 m 2 16π becomes a fundamental scale that characterizes the features of these solutions.
A glance at the factor that multiplies the two-spheres in the above line element shows that if s = +1 then the area of the two-spheres vanishes at R = R c , which suggests to introduce a new coordinate ρ 2 = R 2 − R 4 c /R 2 that allows to rewrite the line element in such a way that the center would be located at ρ = 0. By contrast, for s = −1 we find that the two-spheres have a non-vanishing minimum area at R = R c , which can be interpreted as representing the throat of a wormhole (for a detailed account of wormhole physics see [50]) or, alternatively, as a black bounce [35][36][37]. The two signs, therefore, describe quite different objects, namely, point-like particles if s = +1 and wormholes if s = −1. In what follows, we will study the properties of these once-center extremal solutions 4 .

Curvatures of individual centers
A look at the curvature scalars provides a useful comparison with the GR solutions. In GR (s = 0), the Ricci scalar vanishes and the Kretschmann diverges at R = 0 as ∼ m 4 /R 8 . When s = +1, we find that as R → R c the Ricci scalar goes like ∼ 1/(R − R c ) 2 and the Kretschmann as ∼ 1/(R − R c ) 4 , having a softer behavior if R c = m, where they become (2m) −1 /(R − R c ) and (2m) −2 /(R − R c ) 2 , respectively. On the other hand, if s = −1, we find that as R → R c the Ricci scalar goes like ∼ 1/(R − R c ) 3 and the Kretschmann as ∼ 1/(R − R c ) 6 , having also a softer behavior if R c = m, where they become − 3 2 (2m) −1 /(R − R c ) and 9 4 (2m) −2 /(R − R c ) 2 , respectively. Since this s = −1 case represents a wormhole, with its throat at R = R c , it is also relevant to look at the curvature scalars in the limit R → 0, where the area of the two-spheres goes to infinity again. In this region, we find that both the Ricci and the Kretschmann scalars are finite, taking the values 36m 2 /R 4 c and 408m 4 /R 8 c , respectively. We thus see that in all cases the relevant curvature divergences of these solutions are much weaker than in the GR case.

Horizons of individual centers
For static individual centers, spherical symmetry allows us to identify the location of horizons by finding x =constant hypersurfaces with vanishing norm.
Regarding such horizons, the case s = +1 has the same structure as the GR solution, with a degenerate horizon located at R = m, whenever m > R c (equivalently m > (|ǫ|κ 2 /16π) 1/2 ). If m < R c (equivalently m < (|ǫ|κ 2 /16π) 1/2 ) then there is no degenerate horizon because the area of the two-spheres vanishes at R = R c > m and the geometry cannot be extended further below. Note that for such small mass configurations if one assumes that |ǫ| ∼ l 2 P lanck , then m m P lanck . In the case s = −1, it is important to note that precisely at R = R c , where the two-spheres reach their minimum, the g tt component vanishes while g RR diverges, Now, in order to avoid the coordinate singularity at R = R c due to the vanishing of the off-diagonal term dvdR when s = −1, a further change of variables is necessary to absorb the factor 1 − with the plus sign corresponding to R > R c and the minus sign to R < R c to guarantee that y(R) is a monotonic function on all the domain R ∈]0, ∞[ (see Fig.1). The explicit relation between y and R can be integrated as with y(R c ) ≡ y c = 4R c /3. This construction guarantees the continuity and derivability of the change of coordinates, as can be seen in Fig.1. Note that in the neighborhood of R c , one has y ≈ y c ± 2 Rc (R − R c ) 2 . As a result of this change of coordinates, (41) becomes which covers the whole range y ∈] − ∞, +∞[ with a single chart and avoids coordinate metric singularities everywhere (except, obviously, in the angular sector). In these coordinates one finds that the normal vector to the hypersurfaces y =constant becomes null at R = m (if m ≥ R c ) and at R = R c (always). The time-like Killing vector χ = ∂ v also has vanishing norm there, confirming that these two locations represent Killing horizons (if m ≥ R c ). While at R = m one generically finds zero surface gravity, this quantity diverges at R c as The reason lies on the fact that the surface gravity is defined as the limit |κ| = lim c R 4 (y) φ 2 and y H represents the location of the horizon. Since near the R = R c horizon one finds that A(y) ≈ 4(R − R c )/R c = (y − y c )/2R c has a square root dependence on (y − y c ), then its derivative necessarily induces a divergence in the denominator of κ. Only when R c = m does κ vanish at R c . The divergence of the surface gravity in this one-center solution is a generic property that only disappears when a specific charge-tomass relation is satisfied, which includes the case R c = m but also other cases with non-vanishing temperature [51].
To deepen into the physical meaning of this infinite surface gravity, it might be useful to have a look at the properties of the matter field there. Given that for the GR solution one hasK = m 2 /R 4 , one readily finds that where ρ represents the field energy density. When s = −1, it is evident that both the electric field intensity squared and the electromagnetic Lagrangian (28) diverge at R = R c , but the energy density is finite and reaches its maximum value there. On the other hand, for s = +1 all those quantities have the reversed behavior, namely, the electric field intensity squared and electromagnetic Lagrangian (which is related to the traversal pressures of the field) are finite and well behaved, though the energy density diverges. We thus see that the divergence of curvature scalars at x = R c is totally uncorrelated with the behavior of the matter fields, which may or may not be divergent at that location. It is thus unclear if there is any physical reason or implication for the divergence of the surface gravity at the wormhole throat when R c = m.

Geodesic structure of individual centers
For spherically symmetric spacetimes the geodesic equation can be written in a simple form [52], which for the line element (40) becomes where u is here the affine parameter and E the energy per unit mass. As usual, this geodesic equation is akin to the motion of a single particle in a one-dimensional effective potential, which in the present case reads explicitly with k = 0, −1 for null and timelike geodesics, respectively, and L denotes the angular momentum per unit mass. Let us consider first the case s = +1. As we approach the center, R → R c , the effective potential is dominated by the term which represents a divergent barrier for any nonzero angular momentum [see Fig.2]. Thus, light rays and massive particles will bounce before reaching the center and will never get there. For exactly radial motions, L = 0, then V ef f ≈ −2k(R c − m) 2 /R 2 c , which vanishes for null rays and becomes a finite positive barrier for time-like observers (k = −1). Thus, particles with E 2 > 2(R c − m) 2 /R 2 c will be able to reach the center in finite affine time. If m > R c it is unclear if the curvature divergence at the center will do any harm to infalling observers because physical objects have a finite extension and, as we have just seen, any nonzero L will experience a bounce before reaching the center. One expects that the internal forces that keep the object cohesive will slow down the infinitesimal elements falling radially (L = 0) and make them bounce with the rest of the body. Nonetheless, a detailed analysis of the behavior of geodesic congruences and/or the interaction of waves with these objects would be necessary to further explore the regularity of these solutions [53]. In any case, the fact that the energy density is divergent at the center suggests a physical obstruction to the extension of null radial geodesics in these configurations.
Let us now consider the wormhole case, s = −1. The shape of the effective potential is depicted in Fig.2 for the case of time-like (k = −1) radial (L = 0) geodesics, though we point out that for geodesics with non-zero angular momentum the qualitative behaviour is similar. Only those geodesics with energy E larger than the maximum of the effective potential will be able to go to the interior region of these solutions and interact with the wormhole throat. In the limit m ≫ R c this maximum is located at R ≈ 3 1/4 R c and grows as V max ≈ 2m 2 In Fig.2], leading there to (u c being an integration constant) for all null (with angular momentum) and time-like geodesics (recall Eq. (42)). This solution is exactly the same as the one corresponding to null radial geodesics (k = 0, L = 0), for which V ef f = 0 everywhere, and indicates that null rays and massive particles reach the surface R = R c in finite affine time. Whether geodesics can be extended beyond this point is a subtle issue with no straightforward answer. On the one hand, Eq. (52) shows that the relation between the radial coordinate y and the affine parameter u is smooth across the throat. However, one could argue that the divergence of curvature scalars (and of the surface gravity) at R c should preclude the extensibility of geodesics, though there are examples which contradict this view [53][54][55][56] (see also [57][58][59]) based on the fact that the energy density at R c is finite. In addition, it has been shown in [51] that radial null geodesics in the one-center solutions are insensitive to the charge and mass parameters, existing a specific combination of them for which all curvature scalars are finite. Since there are no reasons to believe that geodesics should not be extended in that specific case and the geodesics satisfy exactly the same equation, it seems natural to conclude that they are always extensible across x = R c . In Fig.3 we illustrate the case of time-like (k = −1) radial (L = 0) geodesics assuming that they can be extended across the wormhole throat R = R c (note that there is no numerical problem in doing so).
Provided that we have continued the geodesics across R = R c , to see what happens with them in the asymptotic region R → 0 we need to study the effective poten-   (49) with the effective potential (50), taking E = 6, m/Rc = 10 and uc = 2. Note that in this figure the wormhole throat is located at R = Rc = 1. As is apparent from this figure, u can be indefinitely extended across the wormhole throat towards u → +∞ at R = 0.
tial (50) that for massive particles goes like while for null non-radial geodesics (k = 0, L = 0) one has In both cases, the effective potential is attractive and divergent. In the latter case, one finds that the geodesic equation (49) can be integrated as while for the former (k = −1) we have instead where the minus sign in (55) corresponds to outgoing geodesics and to ingoing geodesics in (56). Both of the above expressions show that as R → 0 the affine parameter diverges, u → ±∞ [see Fig.3], confirming in this way that all such geodesics are complete from that side since they would take an infinite affine time to get there. Let us finally point out that the case of radial null geodesics is trivial, since from Eq.(49) one immediately finds that and y ∈] − ∞, ∞[ (as follows from (43)). Thus, the region corresponding to R → 0 (or y → −∞) represents a regular boundary of this spacetime, since it cannot be reached in finite affine time by any geodesic.

V. CONCLUSIONS
In this work we have considered a particular class of the Majumdar-Papapetrou family of solutions of GR, given by a collection of extreme black holes in equilibrium (multicenter solutions), to build a new family of multicenter solutions within modified gravity. The gravitational model considered here is the Eddington-inspired Born-Infeld gravity, a particular member of the Ricci-based gravities class, while to generate such solutions we have made use of a new powerful tool dubbed as the mapping method. We have shown that via this method, solutions of the Einstein-Maxwell system are mapped into those of the EiBI gravity theory coupled to a Born-Infeld-type electrodynamics. We then used the MP multicenter solution on the GR side as the seed to generate a new family of solutions in our modified gravity model. The resulting solutions represent a collection of exotic new objects in equilibrium with important novelties as compared to those of GR, whose properties critically depend on the sign of the EiBI parameter ǫ.
On the one hand, when ǫ > 0, one finds a family of point-like objects which look like extremal black holes for masses above the Planck scale but which loose their horizon in the lower part of the mass spectrum. These objects have bounded electric field at their centers but divergent energy density.
On the other hand, the interpretation of the case ǫ < 0 requires to play a bit to cast the corresponding line element in suitable coordinates. There we have shown that for the case of individual centers they represent a kind of extremal black hole with an internal wormhole, whose throat represents a Killing horizon with divergent surface gravity and curvature. In the interior region, geodesics take an infinite affine time to get to R → 0, which represents a boundary of infinite area. We have interpreted these solutions as geodesically complete because of results previously obtained in [51], though the coexistence of a finite energy density at the throat with a divergent electric field intensity is certainly inconvenient. Let us point out that asymmetric wormhole solutions with incomplete geodesics have been found in the EiBI theory coupled to scalar fields [30,60], but in such a case the matter fields at the wormhole throat are well behaved while it is the innermost region which becomes problematic. These results defy the intuitive notion that wormholes always allow for the completeness of geodesics, and raises further questions on the suitability of different markers to characterize pathologies in the characterization of spacetimes (for a broad discussion on singularity regularization see [61,62]).
Using the methods presented in this work, it is possible to generate multicenter solutions in other gravity theories, such as f (R) and others, coupled to nonlinear electrodynamics theories whose form depends on the specific target gravity theory chosen. One thus may wonder if it would be possible to find multicenter solutions in new gravity theories coupled to Maxwell electrodynamics. This question is relevant because, in particular, in the EiBI case coupled to Maxwell electrodynamics, individual wormhole solutions are known and are better behaved than those found here [53][54][55][56] (some of them are traversable, not hidden behind an event horizon). Thus there is the hope of generating multiwormhole traversable configurations on which ideas related to quantum entanglement and the ER=EPR correlation could be tested [63]. Attempts to find such solutions have also been carried out by the authors but without success so far. The key reason for this negative result seems to lie on the lack of specific symmetries in the field equations of the Ein-stein+nonlinear electrodynamics system. On the other hand, applications of these new multicenter configurations within the context of supersymmetric theories is yet to be explored. Research in these directions is ongoing and we hope to be able to report on them elsewhere.