What is a singular black hole beyond General Relativity?

TL;DR

This paper investigates singular black hole spacetimes in quadratic $f(R)$ gravity coupled with an anisotropic fluid, revealing diverse configurations and challenging the notion of singularities in quantum gravity.

Contribution

It introduces minimal extensions of GR using metric-affine $f(R)$ gravity with anisotropic fluids, exploring their impact on black hole singularities and horizon structures.

Findings

Discovery of configurations with and without wormholes

Identification of solutions with de Sitter interiors

Presence of up to four horizons in certain models

Abstract

Exploring the characterization of singular black hole spacetimes, we study the relation between energy density, curvature invariants, and geodesic completeness using a quadratic $f(R)$ gravity theory coupled to an anisotropic fluid. Working in a metric-affine approach, our models and solutions represent minimal extensions of General Relativity (GR) in the sense that they rapidly recover the usual Reissner-Nordstr\"{o}m solution from near the inner horizon outwards. The anisotropic fluid helps modify only the innermost geometry. Depending on the values and signs of two parameters on the gravitational and matter sectors, a breakdown of the correlations between the finiteness/divergence of the energy density, the behavior of curvature invariants, and the (in)completeness of geodesics is obtained. We find a variety of configurations with and without wormholes, a case with a de Sitter interior, solutions that mimic non-linear models of electrodynamics coupled to GR, and configurations with up to four horizons. Our results raise questions regarding what infinities, if any, a quantum version of these theories should regularize.

Tables

Table 1. Table 1: Summary of the features of the four families of configurations studied in Sec. V (the case { σ > 0 , β < 0 } formulae-sequence 𝜎 0 𝛽 0 \{\sigma>0,\beta<0\} with γ = 1 𝛾 1 \gamma=1 hides some peculiarities, see Sec. V.1 , so it is not contained in this table). The metric components, energy density and curvature scalars refers to the behaviour at the wormhole throat (when it exists) or otherwise at the innermost region of the solutions. Incomplete geodesics refer to the existence of (at least) a single incomplete null or timelike geodesic curve. The breakdown of the correlations among these three concepts is clear.

	Wormhole	Metric components	Energy density	Curvature scalars	Geodesics
Case I	YES if $\gamma >1$	Divergent	Finite	Divergent	Complete if $\gamma >1$
()	NO if				Incomplete if
Case II	NO	Finite	Divergent	Finite	Incomplete
($\sigma >0,\beta >0$)
					Incomplete? if ${\delta }_1>{\delta }_c^{\left(\gamma \right)}$
Case III	NO	Divergent if ${\delta }_1\ne {\delta }_c^{\left(\gamma \right)}$	Finite	Divergent if ${\delta }_1\ne {\delta }_c^{\left(\gamma \right)}$	Complete if ${\delta }_1={\delta }_c^{\left(\gamma \right)}$
()		Finite if ${\delta }_1={\delta }_c^{\left(\gamma \right)}$		Finite if ${\delta }_1={\delta }_c^{\left(\gamma \right)}$	(de Sitter core)
					Incomplete if
Case IV	YES	Divergent	Finite	Divergent	Complete
()

Equations

$$S=\frac{1}{2{\kappa }^2}\int d^{4}x\sqrt{-g}f(R)+S_m(g_{\mu \nu },{\psi }_m),$$

$$f_R{R}_{\mu \nu }-\frac{f}{2}{g}_{\mu \nu }$$

$${\nabla }_{\lambda }^{\Gamma }(\sqrt{-g}{f}_R{g}^{\mu \nu })$$

$$R{f}_R-2f={\kappa }^{2}T,$$

$$G_{\mu \nu }={\kappa }^2{\tau }_{\mu \nu },$$

$${\tau }_{\mu \nu }$$

$${\nabla }_{\lambda }^{\Gamma }(\sqrt{-h}{h}^{\mu \nu })=0,$$

$${\Gamma }_{\mu \nu }^{\lambda }=\frac{h^{\lambda \rho }}{2}\left[{\partial }_{\mu }{h}_{\rho \nu }+{\partial }_{\nu }{h}_{\rho \mu }-{\partial }_{\rho }{h}_{\mu \nu }\right].$$

$$h_{\mu \nu }=f_R{g}_{\mu \nu };h^{\mu \nu }=f_R^{-1}{g}^{\mu \nu },$$

$${R^{\mu }}_{\nu }(h)=\frac{1}{f_R^2}\left(\frac{f}{2}{{\delta }^{\mu }}_{\nu }+{\kappa }^2{T^{\mu }}_{\nu }\right),$$

$${T_{\mu }}^{\nu }=\text{diag}(-\rho ,P_r,P_{\theta },P_{\phi }).$$

$${T_{\mu }}^{\nu }=\text{diag}[-\rho ,-\rho ,K(\rho ),K(\rho )].$$

$${T_{\mu }}^{\nu }=\frac{1}{8\pi }\text{diag}[\phi -2(X{\phi }_X+Y{\phi }_Y),\phi -2(X{\phi }_X+Y{\phi }_Y),\phi ,\phi ],$$

$$d{s}^2=-C(x)d{t}^2+B^{-1}(x)d{x}^2+r^2(x)(d{\theta }^2+{\sin }^2\theta d{\phi }^2),$$

$$r^2(x)=r_0^2\exp \left[-{\int }^{\rho }\frac{d\tilde{\rho }}{\tilde{\rho }+K(\tilde{\rho })}\right],$$

$$K(\rho )=\alpha \rho +\beta {\rho }^2,$$

$$\rho (r)=\frac{(1+\alpha ){\rho }_0}{{\left(\frac{r}{r_0}\right)}^{2(1+\alpha )}-\beta {\rho }_0},$$

$${\rho }_m=\frac{(1+\alpha )}{|\beta |},$$

$$\rho =\frac{{\rho }_m}{z^4-s_{\beta }}.$$

$$f(R)=R-\sigma R^2,$$

$$f_R=1+s_{\beta }{s}_{\sigma }\frac{\gamma }{(z^4-s_{\beta }{)}^2},$$

$$d{s}_h^2=-e^{2\Phi (x)}A(x)d{t}^2+\frac{1}{A(x)}d{x}^2+x^2(d{\theta }^2+{\sin }^2\theta d{\phi }^2),$$

$${R^{\theta }}_{\theta }(h)=\frac{1}{x^2}(1-A-x{A}_x),$$

$$A(x)=1-\frac{2M(x)}{x},$$

$$2\frac{M_x}{x^2}=\frac{1}{f_R^2}\left(\frac{f}{2}+{\kappa }^2{T^{\theta }}_{\theta }\right),$$

$$M(x)=M_0(1+{\delta }_{1}G(x)),$$

$${\delta }_1\equiv \frac{{\kappa }^2{\rho }_m(r_0|\tilde{\beta }{|}^{\frac{1}{4}}{)}^3}{r_S}.$$

$$d{s}^2=-\frac{A(x)}{f_R}d{t}^2+\frac{d{x}^2}{A(x)f_R}+r^2(x)(d{\theta }^2+{\sin }^2\theta d{\phi }^2).$$

$$x^2=f_R(r)r^2,$$

$${\tilde{x}}^2=f_R(z)z^2,$$

$$\frac{dz}{d\tilde{x}}=\frac{1}{f_R^{1/2}\left[1+\frac{1}{2}\frac{z{f}_R,z}{f_R}\right]},$$

Full text

I Introduction

One of the most serious drawbacks associated to Einstein’s theory of General Relativity (GR) is the unavoidable existence, under reasonable physical assumptions, of spacetime singularities deep inside black holes, as well as in the early universe Theorems . This is due to the fact that at such singularities the predictability of physical laws comes to an end because measurements are no longer possible. The underlying reason is that the existence of incomplete geodesics implies the destruction/creation of observers and/or information (light signals) as some limiting boundaries are approached. As a way out of this problem, Penrose introduced CCC the cosmic censorship conjecture, by which singularities emerging out of gravitational collapse are assumed to be hidden behind an event horizon, so they cannot causally affect physical processes taking place in the portion of universe accessible to far away observers. Since sweeping the problem under the carpet does not solve it, finding a consistent description of the interaction between gravity and matter, where the resolution of spacetime singularities may be naturally achieved, has become a major goal from different perspectives (classical and quantum, fundamental and phenomenological).

It is typically argued that spacetime singularities should be resolved by a quantum theory of gravity. This is supported by the idea that the quantum degrees of freedom of the gravitational field are expected to be non-negligible in regions of very high curvature. This view, inherited from the effective field theory approach to quantum theory, is very appealing but should be taken with care in gravitational scenarios, where the notion of singularity is not necessarily tied to the divergence of some quantities in some regions Geroch ; Wald ; Hawking-Ellis . For geometric theories of gravity (classical theories), the very existence of observers is more fundamental than the possibility of obtaining absurd results in a measurement, as the latter is not possible without the former. It is for this reason that the existence of incomplete geodesics, for which the affine parameter is not defined over the whole real line, appears as the key element in the singularity theorems.

In the context of GR, the incompleteness of geodesics usually occurs simultaneously with the divergence of scalar quantities, such as the energy density of the matter sources or certain curvature invariants. These divergences appear as a reason for the incompleteness of the geodesics, leading to a rule of thumb for the identification of singular spacetimes phy (see HS for a critical viewpoint on this issue). Indeed this has shaped many approaches to the singularity problem based on the idea that such quantities should remain bounded (see e.g. Ansoldi for a review).

One of such approaches is given by classical non-linear models of the electromagnetic field. This is supported on the success of Born-Infeld theory of electrodynamics, where a square-root modification of the Maxwell action gets rid of the divergence of the self-energy of Coulomb’s field by imposing a maximum bound on the electric field at the center BI . It is natural to wonder whether a similar mechanism for the removal of singularities could occur in the context of gravitation. In this sense Born-Infeld electrodynamics, though successful in making the energy density of the electromagnetic field finite, fails to keep at bay divergences on the curvature scalars when coupled to gravity, which comes alongside with the incompleteness of (some) geodesics BI-grav . In this regard, similar attempts using other well defined non-linear electrodynamics models have failed as well NED-grav . Nonetheless, it is worth mentioning that some examples of non-linear electrodynamics do regularize curvature divergences AB , but such models are constructed in an ad hoc way and yield unphysical features, as shown by Bronnikov Bronnikov (see also Novello ). This strategy has been extended to the case of gravitational actions going beyond the Einstein-Hilbert Lagrangian of GR, such as Gauss-Bonnet and, more generally, Lovelock theories GB , where similar disappointing results have been obtained (see e.g. GBNED for some attempts in this context). Consequently, it is fair to say that such models have been unable to find a fully consistent way out of the singularity problem in GR.

In this work we shall examine in detail the relation between energy density, curvature invariants, and geodesic completeness in some theories of gravity beyond GR. This will allow us to see if the correlations observed in GR among those quantities still persist in other gravitational theories (see RegularBGR for related ideas explored in this context). In other words, can matter/curvature infinities be seen as the reason for the incompleteness of geodesics? This study is relevant in order to understand what problems, if any, a quantum version111Note that we are assuming that any classical theory of gravity should admit a quantum version. of those theories of gravity should solve.

In our approach, we interpret gravitation as a geometric phenomenon, but geometry as something more than just curvature. In the metric-affine (or Palatini) formulation of classical gravitation, geometric properties such as non-metricity and torsion, besides curvature, are allowed by construction. The lack of these freedoms in the usual Riemannian approach could be an excessive constraint with a potentially non-negligible impact on the problems that gravity theories typically exhibit at high-energy. It should be noted that non-metricity and torsion are necessary to deal with different kinds of geometric defects in continuum systems with a microstructure, such as Bravais crystals or graphene SSP . For this reason, metric-affine geometry is commonly used in the study of condensed matter physics PL . Nonetheless, for operational convenience, in this work we shall neglect torsion (see, however, ORtorsion for a discussion on the role of torsion in metric-affine theories) and focus on non-metricity only Pal . Indeed, the question of whether gravity as a manifestation of the curvature222As a matter of fact, gravity could be interpreted as a manifestation of torsion in a flat background, such as in the teleparallel formulation of general relativity (see e.g. TEGR ), but also it could belong to a more general picture where curvature and torsion are both required to properly describe the gravitational interaction as in the case of Einstein-Cartan theories EC . of spacetime is purely a matter of metrics or if the affine structure of spacetime is on equal footing as the metric one has been at debate since soon after the establishment of GR (see e.g. Zanelli for a pedagogical discussion). Certainly, when GR is formulated à la Palatini, the variation of the action with respect to the independent connection yields a set of equations that simply express the metric-connection compatibility condition. The fact that this approach yields the same dynamics as that of considering the metric as the only independent degree of freedom (metric approach) has frequently lead to regard the Palatini variation as merely an alternative way to deriving the field equations of GR. For other theories of gravity, however, the compatibility between metric and connection is broken and the peculiarities of the metric-affine approach become manifest.

The scenario considered here corresponds to a simple quadratic $f(R)$ gravity extension of GR (for which many applications have been investigated in the literature, see e.g. fRlit ), formulated in a metric-affine framework. It should be pointed out that with the advent of the gravitational wave astronomy following the discovery of GW150914 by LIGO LIGO , both gravitational extensions of GR and exotic compact objects in such models can be put to experimental test ECOs . As the matter sector, in our setup we consider an anisotropic fluid (constrained to satisfy standard energy conditions), which has been recently investigated in some detail in a number of astrophysical/cosmological scenarios af . Such fluids include a number of particularly interesting cases, such as that of non-linear electrodynamics. The resulting spacetimes are split into four different cases, depending on the combinations of the signs of the coupling constant of the quadratic gravity contribution and of a constant associated to the matter sector. A noteworthy feature of many of the solutions obtained is the emergence of a finite-size wormhole structure [see Visser for detailed account on wormhole physics] replacing the point-like singularity typically found at the center of GR black holes. It is worth pointing out that wormholes have been suggested as solutions to spacetime singularities in approaches to quantum gravity such as loop quantum gravity LQG and shape dynamics SD (see also BroFab and references therein, where wormholes are linked to regularization mechanisms.)

The main aim of the present work is to determine when the typically assumed correlation between divergence of curvature scalars and geodesic incompleteness is broken. In this sense, we note that the concept underlying the formulation of the singularity theorems Wald is that of geodesic completeness, namely, whether a geodesic curve can be extended to arbitrarily large values of its affine parameter or not. This is a logically independent and more primitive concept than that of curvature divergences [see phy for a nice discussion on this issue], with the latter playing no role on such theorems. As already mentioned, the widespread identification between them in the literature is explained as due to the fact that in many cases of interest (particularly in GR) those spacetimes having (some) incomplete geodesics, also yield (some) divergent curvature scalars Ansoldi . In some of the spacetimes found here we explicitly show that the presence of wormholes yield geodesically complete spacetimes, though curvature scalars may blow up at the wormhole throat. In other cases without wormholes, we meet the incompleteness of geodesics despite the finiteness of curvature scalars. The relation of these magnitudes with the (boundedness of the) energy density of the matter fields is also discussed.

The paper is organized as follows: in Sec. II we introduce the action and main equations of $f(R)$ gravity formulated à la Palatini. In Sec. III we specify the matter sector of our theory under the form of an anisotropic fluid and introduce a number of constraints on it. Next, in Sec. IV, we focus our discussion upon a quadratic $f(R)$ model and solve the field equations for the metric. Sec. V contains the main results of this work, where we study the four different classes of spacetimes, and discuss in detail the relation between energy density, curvature scalars, and geodesic completeness. We conclude in Sec. VI with a summary and some perspectives.

II Action and main equations

The action of $f(R)$ gravity can be written as

[TABLE]

with the following definitions and conventions: $\kappa^{2}$ is Newton’s constant in suitable units (in GR, $\kappa^{2}=8\pi G/c^{4}$), $g$ is the determinant of the spacetime metric $g_{\mu\nu}$, $f(R)$ is a given function of the curvature scalar, $R\equiv g_{\mu\nu}R^{\mu\nu}$, where the Ricci tensor, $R_{\mu\nu}\equiv R_{\mu\nu}(\Gamma)$, which follows from the Riemann tensor as $R_{\mu\nu}\equiv{R^{\alpha}}_{\mu\alpha\nu}$, is entirely built out of the affine connection, $\Gamma\equiv\Gamma^{\lambda}_{\mu\nu}$, which is a priori independent of the metric (metric-affine or Palatini approach). Finally, $S_{m}$ is the matter action, which is assumed to depend only on the matter fields, collectively denoted as $\psi_{m}$, and on the metric $g_{\mu\nu}$.

Performing independent variations of the action (1) with respect to metric and connection one gets two systems of equations

[TABLE]

where $f_{R}\equiv df/dR$ and $T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta S_{m}}{\delta g^{\mu\nu}}$ is the stress-energy tensor of the matter. It is worth mentioning that Eq. (3) simply states that the independent connection fails to be metric or, in other words, that a non-metricity tensor ${Q_{\lambda\mu\nu}}\equiv\nabla_{\lambda}^{\Gamma}g_{\mu\nu}\neq 0$ is present. In the GR case, $f_{R}=1$ and Eq. (3) becomes $\nabla_{\lambda}^{\Gamma}(\sqrt{-g}g^{\mu\nu})=0$, which is fully equivalent to $\nabla_{\lambda}^{\Gamma}g_{\mu\nu}=0$ and thus $\Gamma_{\mu\nu}^{\lambda}$ becomes the Levi-Civita connection of the metric $g_{\mu\nu}$, while the field equations (2) boil down to those of GR with possibly a cosmological constant term. This is the underlying reason for the equivalence between the Palatini and metric formulations of GR. For more general $f(R)$ Lagrangians, however, non-metricity becomes an inherent feature of the field equations.

It is also important to understand the intimate relation existing between matter and gravity in Palatini theories of gravity. Tracing with $g^{\mu\nu}$ in Eq. (2) yields the result

[TABLE]

where $T$ is the trace of the stress-energy tensor. This is not a differential equation, but instead it just establishes an algebraic, non-linear relation between curvature and matter. Given an $f(R)$ theory, solving Eq. (4) yields a solution $R=R(T)$, which generalizes the GR relation, $R=-\kappa^{2}T$. This algebraic relation explains the absence of extra dynamical degrees of freedom in our theory as compared to the usual metric approach, where the scalar curvature satisfies a second-order differential equation, thus implying the presence of propagating scalar degrees of freedom. In the Palatini case, the additional curvature terms are just nonlinear functions of $T$ and can be collected as extra pieces in an effective stress-energy tensor. This way, the Palatini field equations for the metric (2) can be simply written as

[TABLE]

where the effective stress-energy tensor is written as

[TABLE]

However, from a practical point of view, in many cases of interest it is easier to solve the field equations by noting that the result $R=R(T)$ allows us to introduce in Eq. (3) a rank-two tensor $h_{\mu\nu}$ satisfying

[TABLE]

such that the independent connection $\Gamma_{\mu\nu}^{\lambda}$ can be expressed as the Christoffel symbols of the metric $h_{\mu\nu}$, i.e.,

[TABLE]

Comparing this with Eq. (3), it is immediately seen that the physical metric $g_{\mu\nu}$ can be obtained out of $h_{\mu\nu}$ according to the conformal transformations

[TABLE]

where, recall, $f_{R}$ is a function of the matter, $f_{R}\equiv f_{R}(T)$.

An alternative representation of the field equations is now possible in terms of $h_{\mu\nu}$ by contracting Eq. (2) with $h^{\alpha\mu}$ and using the relations (9) to obtain

[TABLE]

where ${R_{\mu\nu}}(h)\equiv{R_{\mu\nu}}(\Gamma)$ is the Ricci tensor constructed with the Christoffel symbols of the metric $h_{\mu\nu}$, see Eq. (8). Note that due to the fact that $f\equiv f(R(T))$ all the objects on the right-hand-side of Eq. (10) are just functions of the matter. Thus Eq. (10) represents a set of second-order field equations for $h_{\mu\nu}$ and, since the conformal transformations (9) depend only on the matter sources, the field equations for $g_{\mu\nu}$ will be second-order as well. In vacuum, ${T^{\mu}}_{\nu}=0$, one has $g_{\mu\nu}=h_{\mu\nu}$ (up to a trivial re-scaling of units) and the field equations (10) reduce to those of GR with a cosmological constant term, which confirms the absence of ghost-like propagating degrees of freedom in these theories.

III Anisotropic fluids

In this work we are interested on obtaining black hole solutions in Palatini $f(R)$ theories, and to compare their structure with that of electrically charged black holes of GR. However, due to the fact that the non-linear corrections appearing on the right-hand-side of the new gravitational field equations (either in Eq. (5) or Eq. (10)) depend just on the trace of the matter, $f(R)\equiv f(R(T))$, the new dynamics encoded in Palatini $f(R)$ theories can only be excited when non-traceless stress-energy tensors are considered. This implies that considering a classical Maxwell electromagnetic field, whose trace is zero, would yield electrovacuum solutions identical to those of GR with a cosmological constant (Reissner-Nordström-Anti-de Sitter black holes). Thus, in order to explore new physics in these scenarios, we must consider stress-energy tensors with a non-vanishing trace. One can then assume that a trace anomaly or other types of corrections are generated by quantum effects and propose a stress-energy tensor of the following form:

[TABLE]

This corresponds to an anisotropic fluid, where $\rho$ is the energy density and $\{P_{r},P_{\theta},P_{\varphi}\}$ are the (different, in principle) pressures. This class of fluids has been recently considered in Refs. fluids1 ; fluids2 ; Tamang where, working in slightly different scenarios, it was found that wormhole solutions can be constructed333Here the word “constructed” means that the wormhole geometry is given first, and then the gravitational field equations are driven back in order to find the matter sources threading the geometry. This is a widely spread strategy in the context of wormhole physics Visser . in Eddington-inspired Born-Infeld theories of gravity without violation of the energy conditions. In contrast to that approach, as we shall show below, in the Palatini $f(R)$ scenario considered here, wormholes can be obtained directly as solutions of the field equations without a priori designer approach.

Fluid model

To simplify the analysis and obtain analytically accessible scenarios, let us constrain the functions defining our model. First we restrict the fluid to satisfy $P_{r}=-\rho$ and $P_{\theta}=P_{\varphi}=K(\rho)$, where $K(\rho)$ is a free input function whose form will be specified later. Thus, the stress-energy tensor for this fluid reads

[TABLE]

A motivation for considering these constraints is the fact that the form of the stress-energy tensor (12) exactly matches that of some non-linear theories of electrodynamics. Indeed, in such a case, defining the matter model as a given function $\varphi(X,Y)$ of the two field invariants $X=-(1/2)F_{\mu\nu}F^{\mu\nu}$ and $Y=-(1/2)F_{\mu\nu}F^{*\mu\nu}$, that can be built out of the field strength tensor $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ and its dual $F^{\mu\nu*}=\frac{1}{2}\epsilon^{\mu\nu\alpha\beta}F_{\alpha\beta}$, the corresponding stress-energy tensor is written as

[TABLE]

where $\varphi_{X}\equiv d\varphi/dX$ and $\varphi_{Y}\equiv d\varphi/dY$. Identifying $-8\pi\rho=\varphi-2(X\varphi_{X}+Y\varphi_{Y})$ and $8\pi K(\rho)=\varphi$, it is clear that specifying a function $K(\rho)$ allows to solve these equations to determine the function $\varphi(X,Y)$, at least in implicit form, associated to the anisotropic fluid under consideration.

To obtain additional information on the fluid described by the stress-energy tensor (12), using the fact that the independent connection $\Gamma_{\mu\nu}^{\lambda}$ does not couple to the matter in the action (1), one finds that the standard conservation equation, $\nabla_{\mu}{T^{\mu}}_{\nu}=0$, holds in these theories. Now, considering static spherically symmetric spacetimes, we can write a line element for the spacetime metric $g_{\mu\nu}$ as

[TABLE]

where the functions $C(x)$, $B(x)$ and $r(x)$ are to be determined by integration of the gravitational field equations.

With this line element, the conservation equation above just reads $\rho_{x}+2[\rho+K(\rho)]r_{x}/r=0$, where $\rho_{x}\equiv d\rho/dx$ and $r_{x}\equiv dr/dx$, which can be integrated to give a relation between $r(x)$ and $\rho(x)$ as

[TABLE]

where $r_{0}$ is an integration constant with dimensions of length and $\tilde{\rho}$ is the energy density without dimensions. To proceed further and integrate explicitly this equation, we need to specify a function $K(\rho)$. Let us take the choice

[TABLE]

where, for dimensional consistency, $\alpha$ is a dimensionless constant and $\beta$ has dimensions of inverse density. This choice covers a number of interesting cases and allows us to obtain analytical solutions. Indeed, in this case, from the expression (16), the relation between $\rho(x)$ and $r(x)$ in Eq. (15) is explicitly written as

[TABLE]

where $\rho_{0}$ is a reference energy density that arises as an integration constant and can be fixed from the asymptotic behavior of the fluid. In particular, for $\alpha=1$, the fluid density and the metric far from the center tend to those generated by a Maxwell field, namely, $\rho r^{4}=\frac{q^{2}}{8\pi}$, which allows to relate $\rho_{0}r_{0}^{4}$ with the electric charge, $q$. Moreover, if $\beta=0$, the stress-energy tensor of the fluid exactly becomes that of a Maxwell field with a vanishing trace and, as already mentioned, this yields the same dynamics as that of GR. However, non-trivial combinations of $\alpha$ and $\beta$ provide modified field equations and generate new solutions.

The analysis now requires to be split into the cases $\beta<0$ and $\beta>0$, since their properties are very different. For $\beta>0$ there is a critical radius $r_{\star}=(|\beta|\rho_{0})^{1/[2(1+\alpha)]}r_{0}$ at which the energy density blows up. Thus the location of the standard divergence in the density of the fluid (Maxwell case) shifts from $r=0$ to the finite radius $r_{\star}$. On the other hand, for the case $\beta<0$ the energy density is finite everywhere, having a maximum value

[TABLE]

at the center. This is quite a similar result as that found in certain models of non-linear electrodynamics, such as the one of Born and Infeld BI , where the electric field attains a maximum value at the center and regularizes the energy density. In Sec.V we will study the implications and impact of the finiteness (or not) of the energy density, via the bound (18), on the regularity of the corresponding spacetimes. Note in this sense that the particular case with $\beta=0$ and $0<\alpha<1$ was studied in detail in Ref.Universe .

To simplify the analysis and the notation let us fix $\alpha=1$ from now on and define $\tilde{\beta}=s_{\beta}|\beta|\rho_{0}$, with $s_{\beta}=\pm 1$ denoting the sign of $\beta$, and introduce the dimensionless variable $z=r/r_{\star}$, with $r_{\star}$ the critical radius defined above. Then, we get $z^{4}=s_{\beta}r^{4}/\tilde{\beta}r_{0}^{4}$ so that the energy density of the fluid simply reads

[TABLE]

To conclude this section, we emphasize that we are only considering matter sources satisfying the energy conditions. For instance, the weak energy condition (WEC) states that the following conditions have to be fulfilled Visser : $\rho>0$ and $\rho+p_{i}>0$ ($i=r,\theta,\varphi$) in Eq. (11). For the particular ansatz (12) with the choice (16) and the expressions for the energy density (17) and (18), it follows that the WEC will be satisfied whenever $\alpha>0$, which is consistent with the choice $\alpha=1$ above.

IV Gravity model and formal solutions

To work with the simplest possible scenario, let us consider the quadratic $f(R)$ model

[TABLE]

where $\sigma$ is a constant with dimensions of length squared. This model is particularly amenable for calculations because the trace equation (4) yields $R=-\kappa^{2}T$, which is the same linear relation as in GR, this result being just an accident related to the functional form of the quadratic model in four dimensions. With this choice, we find that the quantity $f_{R}$, which will play a key role in the characterization of the solutions, takes the simple form444If in the gravity Lagrangian we allow $\sigma$ to take positive and negative values, then $\gamma$ should be parameterized as $s_{\gamma}|\gamma|$. This leads to four types of models depending on the different combinations of $s_{\beta}$ and $s_{\gamma}$.

[TABLE]

where $\gamma\equiv\rho_{m}/\rho_{\sigma}$ (and we have introduced $\rho_{\sigma}\equiv 1/(8\kappa^{2}|\sigma|)$ to denote the energy scale associated to the gravitational coupling constant $\sigma=s_{\sigma}|\sigma|$) represents the relative strength between the matter and gravitational sectors, such that the GR limit is recovered when $\gamma\to 0$. Note that the parametrization of $\sigma$ with $s_{\sigma}$ and of $\beta$ with $s_{\beta}$ leads to four different configurations, which will be studied separately in Sec.V .

IV.1 The metric

To solve the field equations (10) we introduce a static, spherically symmetric line element for the auxiliary metric $h_{\mu\nu}$ as

[TABLE]

where $\Phi(x)$ and $A(x)$ are two functions to be determined using the field equations (10). From the symmetry ${T^{t}}_{t}={T^{x}}_{x}$ one finds that ${R^{t}}_{t}-{R^{x}}_{x}=0$, which implies that $\Phi(x)=$ constant, which can be put to zero by a redefinition of the time coordinate without loss of generality. The remaining field equation follows from the component

[TABLE]

which can be simplified by introducing the mass ansatz

[TABLE]

leading to the first-order equation

[TABLE]

where $M_{x}\equiv dM/dx$. To handle the integration of the mass function $M(x)$ it is useful to take a parametrization

[TABLE]

with $2M_{0}\equiv r_{S}$ representing the Schwarzschild radius and $\delta_{1}$ a dimensionless constant defined as

[TABLE]

This puts forward that $M(x)$ is made out of a constant contribution, $M_{0}$, plus a term generated by the fluid and represented by the function $G(x)$ (see Eq. (32) below). The resulting solution allows to construct the physical metric $g_{\mu\nu}$ by means of the conformal relations (9). This way, the physical line element can be written as

[TABLE]

Taking now into account that such conformal transformations also imply that

[TABLE]

whose dimensionless version using $z$ and $\tilde{x}\equiv x/r_{\star}$ is

[TABLE]

and then we obtain the relation

[TABLE]

which allows us to express (25), by means of Eq. (26), as a differential equation involving only the variable $z$:

[TABLE]

with $G_{z}\equiv dG/dz$. Therefore, by formally integrating $G_{z}$, the metric component $g_{tt}$ in Eq. (28) is obtained in terms of the radial function $z$ as

[TABLE]

IV.2 Geodesic completeness

The non-trivial modified dynamics induced by the gravitational $R^{2}$ corrections necessarily modifies the geodesic structure of the corresponding geometry as compared to GR solution. This is a question of utmost interest, given the fact that geodesic completeness, namely, whether any (null and timelike) geodesic can be extended to arbitrarily large values of the affine parameter, is the most fundamental and generally accepted criterion to determine whether a spacetime is singular o not Wald . Since timelike geodesics are associated to physical observers and null geodesics to the propagation of information, this criterion captures the intuitive idea that in a physically well behaved spacetime nothing can suddenly cease to exist and that nothing can emerge out of nowhere. Nonetheless, as discussed in the introduction, there is frequently a misunderstanding in the literature, taking curvature divergences as an equivalent concept to that of geodesic completeness in order to detect the presence of spacetime singularities. As we shall show in Sec.V, such an identification explicitly breaks in many of the geometries considered in this work. Thus we are mainly interested in studying the geodesic structure in those cases where the GR geodesics are incomplete and consequently yield a singularity, regardless of the presence or not of curvature divergences. To this end, in this section we shall specify the geodesic equation for Palatini $f(R)$ theories and solutions of the form studied here.

In a coordinate system, a geodesic curve $\gamma^{\mu}=x^{\mu}(\lambda)$ associated to a given connection $\Gamma^{\mu}_{\alpha\beta}$ is defined by the equation Wald

[TABLE]

where $\lambda$ is the affine parameter. Since in the action (1) defining our model, the matter part couples to the metric but not to the connection, we will focus on the geodesics associated to the physical metric $g_{\mu\nu}$, which are the ones that the matter fields follow according to the Einstein equivalence principle (see OlmoBook for an extended discussion on geodesics in metric-affine spaces).

The analysis can be largely simplified by writing the geodesic equation using the tangent vector $u^{\mu}=dx^{\mu}/d\lambda$, which satisfies $u_{\mu}u^{\mu}=k$, with $k=1,0,-1$ corresponding to spacelike, null, and timelike geodesics, respectively. Taking advantage of spherical symmetry, without loss of generality we can rotate the angular plane in such a way that it coincides with $\theta=\pi/2$, which further simplifies the problem. From the line element (28) we can, in addition, identify two conserved quantities of motion, $E=(A(x)/f_{R})dt/d\lambda$ and $L=r^{2}(x)d\varphi/d\lambda$. For timelike geodesics, these quantities carry the meaning of the total energy per unit mass and angular momentum per unit mass, respectively. For null geodesics $E$ and $L$ lack a proper meaning by themselves, but the quantity $L/E$ can be identified as an apparent impact parameter as seen from the asymptotically flat infinity Chandra .

Under these conditions, the geodesic equation (34) for the above geometries simply reads

[TABLE]

By using the relation of coordinates (30) (and also the associated Eq. (31)), we rewrite the geodesic equation (35) as

[TABLE]

where $\lambda$ is measured in units of $r_{\star}$, and the sign $\pm$ corresponds to outgoing/ingoing geodesics, with

[TABLE]

as one can deduce by following the steps of Sec. IV.1. Equivalently, the geodesic equation can be written in the more convenient form

[TABLE]

In the next section we shall study in detail the properties of the four different cases of configurations, corresponding to the combinations of the signs of $\sigma$ and $\beta$, and their respective features regarding the behaviour of the energy density, the curvature scalars, and geodesic completeness.

V Analysis of the solutions

V.1 Case I: $\sigma>0$, $\beta<0$

Let us now particularize the above equations to the case in which $s_{\sigma}=1$ and $s_{\beta}=-1$, for which we obtain

[TABLE]

The function $G(z)$ determined by Eq. (41) can be easily solved using power series expansions, and the resulting solutions can be classified in terms of the values of the parameter $\gamma$ defined in Eq. (21). Depending on whether $\gamma$ is greater or smaller than unity, one finds different families of solutions. In this sense, the behavior of the function $z=z(\tilde{x})$, which arises from the resolution of Eq. (30), contains valuable information. Note that according to Eq. (66) the function $f_{R}$ vanishes at

[TABLE]

which sets a critical value for $\gamma=1$. When $\gamma\geq 1$, the radial function $z(x)$ has a minimum at $z_{c}$ where, according to Eq. (30), $\tilde{x}=0$. (From now on, we drop the tilde from $\tilde{x}$ to lighten the notation). Though a compact expression for $z=z(x)$ is not easy to find in general, a series expansion around $z=z_{c}$ yields the result

[TABLE]

From this expression one finds that $z\approx z_{c}+x^{2}(1+z_{c}^{4})/(8z_{c}^{5})$, which shows that for $x>0$ the area of the $2$-spheres decreases with decreasing $x$, but for $x<0$ increases with decreasing $x$, with a minimum at $x=0$ ($z=z_{c}$). This behavior is clearly seen in Fig. 1 where Eq. (30) has been inverted numerically for several values of $\gamma>1$. The interpretation of this minimal area in the two-spheres is well known in the literature: it represents a wormhole Visser , a topologically non-trivial bridge connecting two asymptotically flat spacetime regions, where $z_{c}$ ($x=0$) sets the location of the throat. As it has been found in other cases of Palatini $f(R)$ theories coupled to various matter sources Universe ; bcor16 , the emergence of this structure is directly related to the existence of zeros in the function $f_{R}$.

To study in more detail the geometry around $z_{c}$, it is useful to consider the following expansions

[TABLE]

where $C=\frac{1}{\sqrt{32}}\frac{z_{c}^{7}}{[z_{c}(1+z_{c}^{4})]^{3/2}}$ is a constant. Upon integration, one finds that

[TABLE]

which diverges at $z=z_{c}$. The $g_{tt}$ component of the physical metric appearing in Eq. (33) can thus be approximated as

[TABLE]

Due to the divergence in this metric component, curvature scalars generically diverge near $z=z_{c}$. Obviously, this behavior is shared by all those models in which the function $f_{R}$ has a single pole at $z=z_{c}$. It should be noted, however, that curvature divergences are not synonyms with spacetime singularities, as mentioned in the introduction (see Sec.IV.2 below).

Now, let us analyze the case $0<\gamma<1$ for which $z_{c}$ in (42) has no real solutions. In this case, $z$ belongs to the range $(0,+\infty)$ and no wormhole geometries are found. Near the center, $z=0$, the relation (30) can be expanded as

[TABLE]

while $f_{R}$ becomes there

[TABLE]

and the function $G_{z}$ is finite

[TABLE]

In this case, near the origin $G(z)$ can be approximated as

[TABLE]

where $\delta_{c}^{(\gamma)}$ is a constant (different for each value of $\gamma$) whose value guarantees that the Reissner-Nordström solution of GR is recovered in the far limit, $z\gg 1$. The expansion of the metric component $g_{tt}$ around the center is

[TABLE]

Note that for the choice $\delta_{1}=\delta_{c}^{(\gamma)}$, the metric is finite everywhere. On the other hand, for $\delta_{1}>\delta_{c}^{\gamma}$ the metric at the center is divergent and timelike, while for $\delta_{1}<\delta_{c}^{\gamma}$ it becomes spacelike. Nonetheless, no matter the behaviour of the metric at the center, in all cases curvature invariants such as $K\equiv{R_{\alpha}}^{\beta\mu\nu}{R^{\alpha}}_{\beta\mu\nu}$ do always have divergences at $z=0$. The behaviour of the metric at the center also determines the number (and type) of the horizons, mimicking the basic description of some models of nonlinear electrodynamics BI-grav ; NED-grav ; dr : two, one (degenerate) or no horizons for $\delta_{1}>\delta_{c}^{\gamma}$; a single non-degenerate horizon if $\delta_{1}<\delta_{c}^{\gamma}$, and no horizons if $\delta_{1}=\delta_{c}^{(\gamma)}$.

Finally, the critical case $\gamma=1$ must be treated separately, leading to

[TABLE]

Thus the metric diverges at $z=0$, which induces the presence of curvature divergences there. Whether a wormhole exists in this case or not is a matter of taste, as its throat would have vanishing area:

[TABLE]

To summarize the results obtained so far, we can say that when the matter density scale $\rho_{m}$ is larger than the gravity scale $\rho_{\sigma}$, i.e., $\gamma=\rho_{m}/\rho_{\sigma}>1$, the theory yields wormhole solutions. Whether this wormhole is hidden behind an event horizon or not depends on the combination of parameters $\gamma,\delta_{1},r_{S}$ characterizing the solutions. However, a detailed analysis of the horizon structure of these solutions is beyond the purpose of the present work. We just mention that the geometry is almost identical to the Reissner-Nordström solution of GR everywhere except in the region within the inner horizon555Note, in this sense, that in the asymptotic limit $z\gg 1$, and for arbitrary $\gamma$, we have $f_{R}\simeq 1$ , so that $z^{2}\approx\tilde{x}^{2}$ (we explicitly reintroduce the tilde here) and the role of $z$ as the radial coordinate in GR is restored, while the function $G(z)$ in Eq. (32) quickly converges to the GR solution, $G_{z}\approx 1/z^{2}$, thus recovering the Reissner-Nordström geometry of GR., where some departures arise and modify the structure of horizons. When the matter density scale $\rho_{m}$ is lower than the gravity scale $\rho_{\sigma}$, the wormhole throat closes $\gamma=1$ and no wormhole solution exist anymore $\gamma<1$.

Geodesic structure

Let us begin by noting that regardless of the value of $\gamma$, far from the center ($z\to+\infty$), $f_{R}\approx 1$ and $G_{z}\approx 1/z^{2}$. This means that in that region the GR solution is recovered and the geodesics are essentially coincident with those of GR there. One can verify numerically that this approximation is valid (almost exact!) for all configurations with $r_{S}/r_{\star}\geq 10$ and arbitrary $\delta_{1}$.

Let us focus first on the wormhole configurations, $\gamma>1$, for which our main concern is to study the deviations in the behavior of geodesics near the throat, located at $z=z_{c}$. Consider first radial null geodesics ($k=L=0$). Near the wormhole throat $z\to z_{c}$ the geodesic equation (38) becomes

[TABLE]

where $\hat{C}=\sqrt{1+z_{c}^{4}}/(4E\sqrt{2z_{c}})$. By direct integration, we find

[TABLE]

From this expression it follows that as $z\to z_{c}$ one has $\lambda\to\pm\infty$. Stated in words, this means that ingoing light rays, emitted from $z\to+\infty$ when $\lambda\to-\infty$, approach the wormhole at $z\to z_{c}$ as $\lambda\to+\infty$, while outgoing light rays, which propagate to $z\to+\infty$ as $\lambda\to+\infty$, set off from the wormhole at $z\to z_{c}$ and $\lambda\to-\infty$. A complete representation of the radial null geodesics is shown in Fig. 2. From this plot one verifies that the far limit recovers the GR behavior while near the throat the affine parameter diverges, guaranteeing in this way the completeness of these geodesics, in agreement with the analysis of the asymptotic behaviors provided above.

Let us now consider nonradial geodesics ($L\neq 0$) and/or timelike geodesics ($k=-1$). Since the left-hand side of (35) is positive by construction, physical trajectories must preserve the positivity of the right-hand side. From the expansions (44) and (46) it follows that

[TABLE]

which diverges to $+\infty$ as $z\to z_{c}$. As a result, the right-hand side must vanish at some $z>z_{c}$, forcing in this way the bounce of these curves and preventing them from reaching the wormhole throat. This is analogous to the behavior observed in the Reissner-Nordström solution of GR, where all such geodesics meet an infinite potential barrier generated by the central object Chandra and never reach the central singularity. We thus conclude that all null and timelike geodesics are complete in this wormhole spacetime666One can check that spacelike geodesics are also complete in these wormhole spacetimes. Some of them never reach the wormhole but others can go through it. The latter correspond to $E=0=L$..

Let us recall that at the wormhole throat, $z=z_{c}$, curvature divergences arise. However, due to the fact that radial null geodesics take an infinite affine time to reach the throat, this implies that they lie at the boundary of the spacetime and do not belong to the physically accessible region. This way such divergences have no influence upon physical observers and there is no need to invoke any cosmic censorship conjecture or similar arguments to hide such configurations behind an event horizon. Since the wormhole throat cannot be causally reached in finite affine time, these results put forward the existence of explicit examples where the presence of curvature divergences do not unavoidably entail singular solutions.

Let us now consider those cases with $0<\gamma<1$ for which no wormhole structure was found. From Eqs. (49) and (51), we obtain for $z\to 0$ that

[TABLE]

The full discussion of the geodesic structure would proceed now in much the same way as in the case of certain models of non-linear electrodynamics coupled to GR, where the nature of the central region (spacelike or timelike), which depends on the ratio $\delta_{1}/\delta_{c}^{(\gamma)}$, will determine the type of geodesic able to approach the innermost region. Nonetheless, it is enough to consider radial null geodesics, for which the geodesic equation (36) reads

[TABLE]

This equation can be readily integrated, $\lambda(z)=\lambda_{0}\pm\frac{z}{E(1-\gamma)^{1/2}}$, implying that the origin can be reached in a finite affine time, without possibility of further extension. This result is identical to that found in the GR case, which is regarded as singular, but is in sharp contrast with the previous results for the wormhole case.

Finally, for the transition case $\gamma=1$ (with $z_{c}=0$), the expansion of the metric as $z\to 0$ yields

[TABLE]

Radial null geodesics satisfy the following equation

[TABLE]

and one can easily find that $\lambda(z)\approx\mp\frac{r_{\star}}{E}\frac{1}{\sqrt{2}z}+\lambda_{0}$, which puts forward that these geodesics are complete, as they take an infinite affine time to reach the center. However, this model hides an unusual complexity (see Fig. 3). Indeed, if one considers nonradial and/or timelike geodesics, configurations with $\delta_{1}\geq\delta_{c}^{(1)}$ lead to a bounce at some $z>0$, while for those with $\delta_{1}<\delta_{c}^{(1)}$ geodesics take a finite affine time to reach the origin. Thus, despite the completeness of radial null geodesics and of those with $\delta_{1}\geq\delta_{c}^{(1)}$, the case $\gamma=1$ may lead to geodesically incomplete configurations, depending on the values of the parameters.

V.2 Case II: $\sigma>0$, $\beta>0$

Let us now shift our attention to the case in which both $s_{\sigma}$ and $s_{\beta}$ are positive. Then we find

[TABLE]

An important difference as compared to the previous case is that, regardless of the value of $\gamma$, $f_{R}$ and $G_{z}$ diverge as $z\to 1$, where the energy density $\rho$ becomes infinite. Whether these divergences imply that the spacetime is singular or not is something nontrivial which must be determined after a careful scrutiny of the geometry and its geodesic structure. But before getting into that, one should note that the relation between the coordinates $x$ and $z$, determined by $x^{2}=z^{2}f_{R}$, now is not monotonic, having a minimum as shown in Fig. 4.

Unlike in the $\beta<0$ case, now the minimum is in the function $x=x(z)$ rather than in $z(x)$. This puts forward that now it is the auxiliary geometry which is of the wormhole type. This minimum represents the throat of the wormhole and its location is given by $x_{min}=2z_{min}^{3}/\sqrt{1+3z_{min}^{4}}$, where $z_{min}$ is related to $\gamma$ via $\gamma=(z_{min}^{4}-1)^{3}/(1+3z_{min}^{4})$. It should be noted that this wormhole is not symmetric, having an asymptotically Minkowskian region as $x\approx z\to\infty$ and a non flat region as $x\to\infty$ when $z\to 1$. In fact, as $z\to 1$ in the internal region, $G_{z}\approx\gamma^{\frac{1}{2}}/64(z-1)^{4}$ and $f_{R}\approx\gamma/16(z-1)^{2}$ imply that $G(z)\approx-\gamma^{\frac{1}{2}}/[192(z-1)^{3}]$, which leads to

[TABLE]

Using this relation and noting that as $z\to 1$ we have $x^{2}\approx\gamma/16(z-1)^{2}\to\infty$, one gets

[TABLE]

which is timelike and divergent as $x\to\infty$ ($z\to 1$). The physical metric, on the other hand, has a completely different behavior. Given that $g_{tt}=h_{tt}/f_{R}$, and that $g_{rr}dr^{2}=g_{xx}dx^{2}$, expanding about $z\to 1$ yields

[TABLE]

which are always finite at $z=1$ where, recall, the energy density diverges. In Fig. 5 the behavior of the $g_{tt}$ component is shown for a configuration which exhibits up to four horizons in the $z>1$ interval.

By proceeding in the same way as in the previous case, we now perform the analysis of the geodesic structure. In this sense, Fig. 6 puts forward that a region of infinite energy density is reached by null and timelike radial geodesics in a finite affine time. If the divergence in the matter sector is interpreted as defining a limiting boundary of the physical spacetime, where the equations no longer make sense, then the fact that geodesics can reach it in a finite affine time would imply that this geometry is singular.

From the numerical results shown in Figs. 5 and 6 and, from the above analytical approximations, it is evident that nothing special happens to the physical metric at the points $z=1$ or $z=z_{min}$. This can be further emphasized by looking at the whole line element in the $z\to 1$ region, whose form is

[TABLE]

Using this line element, one readily verifies that all curvature invariants are finite at $z=1$ despite the energy density being divergent at that point. Though this divergence in the matter sector must be seen as a breakdown in the description of the fluid model considered, it serves to illustrate that divergences in the matter sector do not necessarily imply divergent curvature invariants. At the same time, the finiteness of curvature invariants is unrelated to the completeness of geodesics.

Analytical extension to $z<1$

The fact that the energy density diverges at $z=1$ and that it changes sign in the $z<1$ region somehow suggests that the physical region should be restricted to the open interval $z>1$. Given that both null and timelike radial geodesics reach the $z=1$ surface in a finite affine time, if one wants to have a nonsingular spacetime, an artificial wormhole extension attached at $z=1$ should be considered to complete the geodesics. This construction, though mathematically admissible, seems a bit unnatural as compared to the wormhole solutions found in the $\{\sigma>0,\beta<0\}$ case, where the energy density is always finite.

Nonetheless, let us note that the numerical integration of the equations that define the geometry is well defined everywhere except at the point $z\to 0$. Moreover, the fact that $f_{R}$ and $G_{z}$ develop divergences at $z=1$ is not an obstacle to extend the integration across the $z=1$ surface, as a simple change of variables avoids the numerical difficulties associated to the parametrization in terms of $G(z)$. In fact, since $g_{tt}$ is well defined even at $z=1$, one can write a smooth differential equation for that function. The point is that the divergence in $G_{z}$ is somehow compensated by the divergence in $f_{R}$. Considering the ansatz (33) to isolate the function $G(z)$, one can compute its derivative and express it in terms of $g_{tt}$ and its first derivative. The resulting equation can be multiplied by $(z^{4}-1)^{4}$ on both sides to get rid of all the divergent terms. This new equation can be numerically integrated from the region where it coincides with GR down to $z=0$. The result is shown in Fig. 7 and is in complete agreement with Fig. 5 in the overlapping region $z\geq 1$. The corresponding extension of the geodesics appears in Fig. 8. Obviously, this region with $z<1$ does not make any physical sense, as it implies a negative energy density that changes from $+\infty$ on $z\gtrsim 1$ to $-\infty$ on $z\lesssim 1$ and remains negative until $z\to 0$. Nonetheless, the numerical problem is well defined all over the $z>0$ domain.

V.3 Case III: $\sigma<0$, $\beta<0$

When $s_{\sigma}=-1$ and $s_{\beta}=-1$, the model is characterized by

[TABLE]

As follows from these expressions, the energy density and the functions $f_{R}$ and $G_{z}$ are everywhere smooth and finite. One can verify that for $\gamma>4$ the function $x(z)$ has a minimum at $z>1$ given by $\gamma=\frac{(1+z_{min}^{4})^{3}}{3z_{min}^{4}-1}$. For this critical value of $\gamma$, the function $G_{z}$ has a zero at $z_{min}$, while for $\gamma>4$ it has two. The former occurs at $z_{min}$, while the latter appears at $0<z_{max}\leq 1\leq z_{min}$ and represents a local maximum (see Fig. 9).

Unlike in the case of Sec.V.2, the minimum in $x(z)$ cannot be regarded as representing a wormhole in the auxiliary metric because it is not an absolute minimum in the sense that $x$ does not bounce back to infinity after crossing $z_{min}$. In the current case, below $z_{min}$ (see Fig. 9) $x(z)$ grows just until it reaches a new local extremum, now a maximum, after which it goes monotonically towards $x\to 0$ as $z\to 0$. The interpretation of this structure in the auxiliary geometry is unclear and will not be explored further in this paper. Since the existence of an internal wormhole or other structures in the auxiliary geometry does not impose any restriction on the radial function $z$ that defines the area of the two-spheres in the physical spacetime, we will assume that $z$ is naturally defined over the whole region $z\geq 0$.

An expansion of the metric in the $z\to 0$ region leads to

[TABLE]

which shows that $g_{tt}$ diverges at the center if $\delta_{1}\neq\delta_{c}^{(\gamma)}$. When $\delta_{1}=\delta_{c}^{(\gamma)}$, the above expansion is not valid and must be re-evaluated, leading to

[TABLE]

which is completely regular at $z=0$. In this particular case, one verifies that the geometry at $z=0$ becomes of de Sitter type, with $R_{\mu\nu}=\frac{\delta_{c}^{(\gamma)}}{\delta_{2}}g_{\mu\nu}$ and finite curvature scalars. It is also worth noting that by a simple rescaling of the time coordinate, $t\to\sqrt{1+\gamma}\tilde{t}$, the metric $g_{\mu\nu}$ becomes Minkowskian at the center, showing that these coordinates locally represent free falling observers.

By inspecting the geodesic equation (38) shows that if $\delta_{1}>\delta_{c}^{(\gamma)}$ then non-radial and timelike geodesics bounce before reaching to the center (like in the usual Reissner-Nordström black hole of GR). However, if $\delta_{1}<\delta_{c}^{(\gamma)}$ (Schwarzschild like configuration) nothing prevents those geodesics from getting there in a finite affine time. Radial null geodesics also reach the center regardless of the value of $\delta_{1}$. It should be noted that curvature divergences exist at $z=0$ as long as $\delta_{1}\neq\delta_{c}^{(\gamma)}$. But note that these divergences arise despite the fact that the energy density is finite everywhere, as in Case I of Sec. V.1.

On the other hand, when $\delta_{1}=\delta_{c}^{(\gamma)}$, we have seen above that the geometry near the center is of de Sitter type and, therefore, the geodesics reaching there should not experience any pathological effect, being able to cross the center and continue their path777Note in this sense that replacing the point-like singularity by a de Sitter core is a standard strategy in looking for regular solutions in the context of GR, see e.g. AdScore ; Ansoldi .. Configurations with $\delta_{1}=\delta_{c}^{\gamma}$ should thus be regarded as nonsingular. In this sense, the fact that radial null geodesics are insensitive to the value of $\delta_{1}$ suggests that they should always be able to go through the apparently pathological region $z=0$ even when $\delta_{1}\neq\delta_{c}^{(\gamma)}$. This view is further reinforced by the lack of correlation between the behavior of curvature scalars, divergent for $\delta_{1}\neq\delta_{c}^{(\gamma)}$, and the energy density, which is always finite. In view of all this, it is unclear whether one should regard any of these solutions as singular, as is typically assumed in the case of GR. Therefore, further analysis on the impact of curvature divergences on the transit of physical observers and on the scattering of waves in these spacetimes to get deeper into their singular/nonsingular character seems necessary and will be carried out elsewhere.

V.4 Case IV: $\sigma<0$, $\beta>0$

We now consider the last case, $s_{\sigma}=-1$ and $s_{\beta}=1$, which is characterized by

[TABLE]

In this model the energy density diverges at $z=1$ but, fortunately, this surface lies beyond the physically accessible region. Therefore, it can be considered finite beyond $z>1$. Indeed, a glance at the function $f_{R}$ above puts forward that $z(x)$ has a minimum of magnitude $z_{c}=(1+\gamma^{1/2})^{1/4}$ at $x=0$, confirming in this way that this model describes wormholes (see Fig. 10) and that $z>1$ always.

A series expansion of the metric about $z=z_{c}$ leads to

[TABLE]

which is always negative and divergent at the throat, implying the existence of curvature divergences on that surface. With this behavior near the throat, a glance to the geodesic equation (38) indicates that all timelike and non-radial geodesics ($L\neq 0$) bounce before reaching the wormhole, just like in the Reissner-Nordström black holes of GR. For radial null geodesics, however, one finds

[TABLE]

which implies that

[TABLE]

diverges as $z\to z_{c}$, confirming that such geodesics take an infinite affine time to reach (or come out from) the wormhole throat, similarly as in the case $\{\sigma>0,\beta<0\}$ of Sec.(V.1). This model, therefore, always yields geodesically complete, nonsingular spacetimes. Note also that, though the geometry has curvature divergences at $z=z_{c}$, the energy density is finite there, taking the value $\rho(z_{c})=\sqrt{\rho_{m}\rho_{\sigma}}$.

VI Conclusions

VI.1 Overview of results

In this work we have considered the problem of classical, non-rotating black holes in a quadratic Palatini $f(R)$ extension of GR. As matter fields we have taken an anisotropic fluid whose stress-energy tensor, besides satisfying the energy conditions, covers a number of interesting cases, in particular, that of non-linear models of electrodynamics. The latter have been frequently employed in the context of GR in order to take care of the singularity problem, but the strategy applied has been unsatisfactory so far. In contrast, the non-trivial interaction between gravity and matter in our case, which is encoded in the (gravity and matter) parameters $\sigma$ and $\beta$, yields a number of novelties on the corresponding spacetimes. These geometries quickly recover the Reissner-Nordström solution of GR far from the center, but drastically modify the innermost structure. The corresponding analysis is split into four different cases, according to the signs of $\sigma$ and $\beta$.

• •

In the case $\{\sigma>0,\beta<0\}$ two classes of configurations are found, according to a certain critical scale, parameterized by $\gamma=1$. If $0<\gamma<1$, the nature of the central point-like singularity (timelike, spacelike or null), is qualitatively similar to that found in certain families of models of non-linear electrodynamics, and the same applies to the structure of horizons dr . In such a case one finds both divergence of curvature scalars and incompleteness of geodesics, hence these spacetimes are regarded as singular. But when $\gamma>1$ the point-like singularity is shifted to a spherical surface of radius $z_{c}(\gamma)>0$, which corresponds to the minimum of the radial function $z(x)$ and represents the throat of a wormhole. Despite curvature scalars being divergent on this surface, we have explicitly shown that this does not prevent the completeness of geodesics for null and timelike observers, in such a way that the wormhole always lies on the future (or past) boundary or the manifold and cannot be reached in finite affine time by geodesic observers. The limiting case $\gamma=1$ has a number of peculiar features, such as complete radial null geodesics but incomplete non-radial and timelike geodesics for certain configurations (with $\delta_{1}<\delta_{c}^{(1)}$). This allows us to conclude that the cases with $\gamma\leq 1$ represent singular spacetimes, while those with $\gamma>1$ are nonsingular and possess a wormhole structure. Nonetheless, note that in all these cases (both nonsingular and singular) the energy density of the fluid is everywhere finite.

• •

For the case $\{\sigma>0,\beta>0\}$ one finds that it is the auxiliary metric which has a wormhole at $z=1$, while the physical metric becomes finite there, exhibiting up to four horizons. The $z=1$ surface is a region of divergent energy density of the fluid, which is reached in finite affine time by null and timelike geodesics, but nonetheless the curvature scalars are finite there. The metric itself admits an analytical extension to the region $z<1$, but the spacetime would still be geodesically incomplete and, besides, the breakdown of the fluid model at $z=1$ seems to indicate the unphysical character of such an extension.

• •

In the case $\{\sigma<0,\beta<0\}$ no wormhole solution is found and the metric is naturally extended down to $z=0$. An expansion of the physical metric there reveals that curvature divergences are present as long as $\delta_{1}\neq\delta_{c}^{(\gamma)}$, but do disappear for the choice $\delta_{1}=\delta_{c}^{(\gamma)}$, for which the metric at the center has a de Sitter behavior. Nonetheless, in both cases the energy density is everywhere finite. In the finite curvature cases, $\delta_{1}=\delta_{c}^{(\gamma)}$, radial null geodesics reaching the $z=0$ region can clearly continue their path. Given that these geodesics are insensitive to the value of $\delta_{1}$, one is tempted to interpret all such radial geodesics as complete, including the cases with $\delta_{1}\neq\delta_{c}^{(\gamma)}$. This interpretation is appealing but risky, as it would logically lead to the conclusion that the Reissner-Nordström black hole of GR is nonsingular (as far as geodesic motion is concerned). To deepen into this question, an analysis on the impact of curvature divergences on congruences of observers and on the propagation of waves would be necessary along the lines of Refs.sw and congruences .

• •

For the case $\{\sigma<0,\beta>0\}$, it always yields wormhole structures that naturally provide geodesically complete spacetimes, much in the same way as in the $\{\sigma>0,\beta<0\}$ case for $\gamma>1$. Note that this is so despite the fact that curvature divergences are always present at $z=z_{c}$, but the energy density is finite there.

In table 1 we display the most relevant features of the four classes of configurations classified according to the combinations of signs of $\sigma$ and $\beta$. Such features are i) the existence or not of wormhole configurations and the behaviour at the wormhole throat (or at the center when no such wormhole exists) of the following objects: ii) the metric components , iii) the energy density, iv) the curvature scalars, and v) the completeness (or not) of all null and timelike geodesics on such spacetimes. From this table it is clear the breakdown of the correlations between the behaviour of energy density, curvature scalars and completeness of geodesics in different ways.

VI.2 Final comments

The research presented in this work is added to the growing set of results within metric-affine geometries supported by $f(R)$ gravity Universe ; bcor16 , quadratic $f(R,R_{\mu\nu}R^{\mu\nu})$ and Eddington-inspired Born-Infeld gravities sw ; BIreg , where the point-like singularity of GR is replaced by a wormhole structure, which allows geodesics to be complete. The results obtained so far seem to point out two different mechanisms for the resolution of spacetime singularities in this context. (i): For $f(R)$ theories the wormhole lies on the future (or past) boundary of the spacetime, but the fact that geodesics can be indefinitely extended means that these spacetimes are nonsingular according to the criteria employed in the singularity theorems Theorems . ii): For the cases of quadratic and Born-Infeld gravity (which admit extensions to higher hd and lower 3D dimensions with similar results), the wormhole may be reached in a finite affine time by some geodesics, but they are naturally extended through it. It should be stressed that in the latter case one may wonder about the physical impact of curvature divergences on physical (extended) bodies crossing the wormhole throat, since the presence of large tidal forces could rip them apart. In this sense, in four-dimensional Born-Infeld gravity a separate analysis using a congruence of geodesics to model such extended bodies has revealed that curvature divergences occurring at the wormhole throat do not pose any destructive threat, being the transit smooth and weakly affected by the large tidal forces at the throat congruences . The problem of scalar wave scattering off this wormhole turns also to be well posed sw .

To conclude, the results obtained in this paper provide a fauna of examples on which the (typically assumed) correlations among curvature divergences, energy density, and geodesic (in)completeness are explicitly broken. The emergence of structures such as wormholes appears as a key element for the extendibility of geodesics (note that we also found a de Sitter interior which guarantees the extendibility of geodesics in a particular example). These non-perturbative features should be taken into account in the interpretation of the existence of divergences in curvature scalars or in the matter fields. Those infinities are significant in perturbative frameworks, where they signal the end of validity of a certain approximation, but in a non-perturbative scenario their interpretation is unclear. It is worth noting, in this sense, that in systems such as graphene, wormholes have been constructed in the discrete description through a careful design of the lattice, adding a series of heptagonal rings that help join two flat sheets with a carbon nanotube grapheneWH . In the continuized description, the wormhole exhibits a curvature divergence at the throat which does not prevent the study of fermion propagation in the resulting effective geometry. This example is, in our view, very representative of what may correspond to a quantum version of the models presented here: geodesic completeness is guaranteed by non-perturbative aspects (the wormhole), though infinities (in curvature and energy density) may arise because of the lack of control of the underlying microstructure in a certain limit crystals . The behavior of classical and quantum fields in these backgrounds is currently underway.

Acknowledgments

C. B. is funded by the National Scientific and Technical Research Council (CONICET). G. J. O. is supported by a Ramon y Cajal contract, the i-COOPB20105 grant of the Spanish Research Council (CSIC), the Consolider Program CPANPHY-1205388, and the Severo Ochoa grant SEV-2014-0398 (Spain). D. R.-G. is funded by the Fundação para a Ciência e a Tecnologia (FCT, Portugal) postdoctoral fellowship No. SFRH/BPD/102958/2014 and the FCT research grant UID/FIS/04434/2013. This work is also supported by the Spanish grant FIS2014-57387-C3-1-P from MINECO and the CNPq (Brazilian agency) project No. 301137/2014-5. C. B. and D. R.-G. thank the Department of Physics of the University of Valencia for their hospitality during the elaboration of this work. This article is based upon work from COST Action CA15117, supported by COST (European Cooperation in Science and Technology).