Eddington-inspired Born-Infeld gravity. Phenomenology of non-linear gravity-matter coupling

Eddington-inspired Born-Infeld gravity. Phenomenology of non-linear gravity-matter coupling

Paolo Pani CENTRA, Departamento de Física, Instituto Superior Técnico, Universidade Técnica de Lisboa - UTL, Av. Rovisco Pais 1, 1049 Lisboa, Portugal.    Térence Delsate CENTRA, Departamento de Física, Instituto Superior Técnico, Universidade Técnica de Lisboa - UTL, Av. Rovisco Pais 1, 1049 Lisboa, Portugal. Theoretical and Mathematical Physics Dpt.,Université de Mons, UMons, 20, Place du Parc 7000 Mons, Belgium.    Vitor Cardoso CENTRA, Departamento de Física, Instituto Superior Técnico, Universidade Técnica de Lisboa - UTL, Av. Rovisco Pais 1, 1049 Lisboa, Portugal. Department of Physics and Astronomy, The University of Mississippi, University, MS 38677, USA.
Abstract

Viable corrections to the matter sector of Poisson’s equation may result in qualitatively different astrophysical phenomenology, for example the gravitational collapse and the properties of compact objects can change drastically. We discuss a class of modified non-relativistic theories and focus on a relativistic completion, Eddington-inspired Born-Infeld gravity. This recently proposed theory is equivalent to General Relativity in vacuum, but its non-trivial coupling to matter prevents singularities in early cosmology and in the non-relativistic collapse of non-interacting particles. We extend our previous analysis, discussing further developments. We present a full numerical study of spherically symmetric non-relativistic gravitational collapse of dust. For any positive coupling, the final state of the collapse is a regular pressureless star rather than a singularity. We also argue that there is no Chandrasekhar limit for the mass of non-relativistic white dwarf in this theory. Finally, we extend our previous results in the fully relativistic theory by constructing static and slowly rotating compact stars governed by nuclear-physics inspired equations of state. In the relativistic theory, there exists an upper bound on the mass of compact objects, suggesting that black holes can still be formed in the relativistic collapse.

pacs:
04.50.-h, 98.80.-k

I Introduction

The beauty and the beast of Einstein’s General Relativity (GR) are encoded in the non-linearity of its field equations. Already in vacuum, GR describes the dynamics of non-linear objects, like black holes Penrose (1996). In order to study the formation of black holes in dynamical situations, e.g. during a stellar collapse, one needs to couple the vacuum theory to matter. How to include this coupling in a proper way was one of Einstein’s main concerns. The requirement of stress-energy tensor conservation, , which in turn implies geodesics motion, together with Bianchi’s identities, , naturally suggests a linear coupling between the Einstein tensor and the stress-energy tensor, . Indeed, under quite generic assumptions (see e.g. D’Inverno (1992)) Einstein’s equations are the most general field equations which involve the stress-energy tensor linearly. Nonetheless, it comes as a surprise that a highly non-linear theory as GR is just linearly coupled to matter. In this paper we investigate modified theories of gravity which account for a non-linear coupling to matter, while retaining many appealing features of Einstein’s theory.

In the weak-field regime, GR has brilliantly passed all experimental and observational scrutinies so far Will (2005); Everitt et al. (2011), but new effects are still viable in the strong-field, non-linear regime. Furthermore, while current experiments have confirmed the weak-field behavior of GR in vacuum or in orbital motion, performing null tests of Einstein’s theory inside matter may be extremely challenging because, due to the equivalence principle, purely gravitational effects are hard to disentangle from those due to non-standard matter coupling. In fact, one usually assumes a minimal coupling, then solves the full Einstein equations in some relevant situation – e.g. in cosmological settings, in the interior of Sun-like stars or inside more compact objects like neutron stars (NSs) – and finally compares theoretical models with observations.

On the other hand, GR suffers from severe theoretical problems and long-standing observational puzzles, which may be precisely related to our poor understanding of the gravity-matter coupling. For example, the dark matter problem (see e.g. Ref. Bertone et al. (2005) for a review) may be explained by invoking new fundamental interactions Bekenstein and Milgrom (1984); Milgrom (1983), rather than assuming exotic particles. Similarly, the cosmological acceleration of the universe may be explained in terms of more complicated interactions, rather than postulating the existence of a mysterious form of dark energy (see e.g. Ref. Clifton et al. (2011) for a recent review on cosmology in modified gravities). These postulates are somehow the modern versions of the aether, suggesting that perhaps something is missing in our understanding of gravitational interactions inside matter.

Likewise, it is well known that the dynamical evolution of matter fields in GR is generically plagued by the formation of singularities (e.g. the Big Bang or those forming in the gravitational collapse of matter fields), which signal a break-down of the theory. According to Penrose’s cosmic censorship Penrose (1969), these singularities must be covered by an event horizon, i.e. a black hole must form as an outcome of the gravitational collapse. However, the cosmic censorship remains a conjecture and its validity (or even its precise formulation) is an open issue.

In a spherically symmetric stellar collapse, locally naked singularities can be generically formed Lemos (1992); Shapiro and Teukolsky (1991a); Joshi and Dwivedi (1993). Initial density gradients can produce shear in the collapsing fluid which, in turn, delays the formation of an apparent horizon, leaving the singularity locally naked Joshi et al. (2002). In some cases, the singularity can even be globally naked, i.e. connected to distant observers by timelike or null-like geodesics, thus suggesting that the censorship can be evaded (see Refs. Joshi (2000); Gundlach (2003); Joshi (2008) for detailed discussions on this topic). However, these “shell-focusing” singularities may be limited to the spherically symmetric case Wald (1997). Relaxing the assumption of spherical symmetry, the situation is even more controversial. While Shapiro and Teukolsky reported numerical evidence for formation of naked singularities in the collapse of highly prolated gas spheroids Shapiro and Teukolsky (1991b), Wald and Iyer Wald and Iyer (1991) pointed out that similar properties – namely the absence of trapped surfaces lying on their maximal slices in a portion of singular spacetime – are also present in a Schwarzschild spacetime with a particular choice of the time slice. These studies show that, although there appears to be growing evidence in support of the cosmic censorship, its validity in GR remains an open issue which is far to be solved Wald (1997); Lehner (2001).

However, we stress here that the issue of singularity formation and the related cosmic censorship strongly depend on how gravity interacts with matter. Thus, the outcome of a gravitational collapse may be different if the gravity-matter coupling is modified.

In this paper, following previous works Vollick (2004, 2005, 2006); Kerner (1982); Banados and Ferreira (2010); Pani et al. (2011a), we take a different perspective and investigate phenomenologically viable modifications to GR that account for non-trivial coupling to matter while being equivalent to Einstein’s theory in vacuum. As a consequence, this class of corrections results in a qualitatively different phenomenology once matter is included (see also Ref. Bertolami et al. (2007)). For concreteness, we will focus on a recent proposal by Bañados and Ferreira Banados and Ferreira (2010) (see also some previous attempt by Vollick Vollick (2004, 2005, 2006)) which we shall call Eddington-inspired Born-Infeld (EiBI) gravity, not to be confused with Eddington-Born-Infeld Banados (2008) or Born-Infeld-Einstein theory Deser and Gibbons (1998). While being conceptually different from each other, these theories share the same determinantal form of Born-Infeld non-linear electrodynamics Born and Infeld (1934); Tseytlin (1999). In this sense, similarly to non-linear electrodynamics, EiBI gravity can be thought as an effective prototypical theory in which resummation of higher order curvature terms – which are qualitatively similar to those expected by a quantum gravity completion – is advocated in order to resolve curvature singularities Klebanov and Strassler (2000). It turns out that these corrections effectively account for some non-linear matter coupling to usual Einstein’s gravity.

Interestingly enough, such corrections can even affect the non-relativistic gravitational regime. We shall discuss the non-relativistic phenomenology in detail and most of our results apply to any theory whose non-relativistic limit is described by the following modified Poisson equation Banados and Ferreira (2010)

 ∇2Φ=4πGρ+κ4∇2ρ, (1)

where is the usual gravitational potential, is the gravitational constant and is the matter density. Within the parametrized post-Poissonian approach proposed in Ref. Casanellas et al. (2011), the equation above is the most general spacially covariant Poisson equation, which is first order in and and reduces identically to standard Laplace’s equation, , in vacuum. While possible second order terms are strongly constrained by tests of the equivalence between gravitational and inertial mass to one part in (see e.g. Ref. Su et al. (1994) and references therein), the term proportional to in Eq. (1) is presently mild constrained.

To our knowledge, null tests of Poisson equation inside matter, in order to constrain the extra term in Eq. (1), have not been conceived yet. Mild observational constraints on come from solar physics observations Casanellas et al. (2011), . This bound has been derived by comparing the observed solar neutrino fluxes and helioseismology observables to the predictions of modified solar models. Such constraints are quite mild, due to some inherent uncertainties on the Sun interior. Very recently, considerably stronger observational constraints were derived in Ref. Avelino (2012), provided direct measures of the central density or of the radius of a NS are available (see also Ref. Pani et al. (2011a)). However, some uncertainty exists on the interior of compact stars and usually GR is assumed to infer the magnitude of the central density, given some “observable” quantity such as the mass and the radius. Thus, table experiments in some controlled setting would be highly desirable.

The scope of this paper is twofold. First, we wish to discuss the phenomenology of Eq. (1) in relation with stellar collapse and static stellar configurations in the non-relativistic limit. Our purpose is to show that, while being compatible with current observations, a simple modified Poisson model leads to a new interesting phenomenology, even for arbitrarily small values of . Secondly, we shall discuss a well-motivated relativistic completion of Eq. (1), EiBI gravity Banados and Ferreira (2010). In Ref. Pani et al. (2011a), we reported some results on the gravitational collapse and compact star solutions in this theory. Here, we give further details and discuss new developments.

The plan of this work is the following. In Sec. II we briefly review some known and new features of EiBI gravity Banados and Ferreira (2010). Sections III and IV are devoted to the phenomenology of the modified Poisson equation (1), which includes the non-relativistic limit of EiBI gravity. Section III is devoted to solve the hydrodynamics equations describing the non-relativistic collapse of non-interacting particles. In pass, we shall correct a wrong result in Ref. Pani et al. (2011a) which, however, does not affect the final picture. In Sec. IV we show that the end-point of the gravitational dust collapse is a regular, pressureless star and we discuss other non-relativistic stellar equilibrium configurations. Finally, in Sec. (V) we discuss relativistic stellar models using nuclear-physics based equations of state. In Sec. VI we draw our conclusion and discuss open issues and future developments.

Ii Eddington-inspired Born-Infeld gravity

EiBI gravity is described by the following action Banados and Ferreira (2010)

 Sg=2κ∫d4x(√|gab+κRab|−λ√−g)+SM[g,ΨM], (2)

where is the matter action, generically denotes any matter field, is the symmetric part of the Ricci tensor built from the connection and is related to the cosmological constant, , so that asymptotically flat solutions are obtained when .

EiBI gravity is reminiscent of a Born-Infeld action for non-linear electrodynamics Born and Infeld (1934); Tseytlin (1999), where here the Ricci tensor plays the role of the field strength . Moreover, the metric and the connection are considered as independent fields Banados and Ferreira (2010), as in Palatini’s approach to GR. Similarly to gravities (see e.g. Sotiriou and Faraoni (2010)) the Palatini formulation is not equivalent to the metric one. In the metric approach, EiBI theory contains ghosts, which must be eliminated by adding extra terms to the action Deser and Gibbons (1998); Vollick (2004). Furthermore, it is assumed that matter minimally couples to the metric tensor only, i.e. the matter action in Eq. (2) only depends on the metric and on the matter fields, but not on the independent connection .

The field equations are conveniently written as Banados and Ferreira (2010)

 Γcab(q) = 12qcd(∂aqbd+∂bqad−∂dqab), (3) qab = gab+κRab(Γ), (4) √−qqab = λ√−g(gab−κTab) (5)

where is an auxiliary metric, is the standard stress-energy tensor and is the symmetric part of the Ricci curvature tensor. Crucially, since matter is minimally coupled to the metric , the field equations above imply the usual conservation of the stress-energy tensor, , where is the covariant derivative written in terms of the metric tensor . Note that is the matrix inverse of and differs inside matter from .

When , the auxiliary metric coincides with the physical metric . Thus, in vacuum is the metric connection. Indeed, in absence of matter fields EiBI theory is equivalent to GR Banados and Ferreira (2010). In presence of matter, the small limit of the field equations reads

 Rab(Γ)=Tab−12Tgab+κ[Sab−14Sgab]+O(κ2). (6)

where . Notice that Einstein’s theory is recovered as . In the field equation above, two qualitatively different corrections to GR are manifest at : those depending on , which are quadratic in the matter fields, and those hidden in , which implicitly depend on derivatives of matter fields. The latter survive even in the non-relativist limit of the theory, which is described by a modified Poisson equation (1).

Since the auxiliary metric is algebraically related to the physical metric , EiBI gravity does not introduce any extra dynamical degree of freedom with respect to GR. In this sense, the theory is similar in spirit to gravities in the Palatini approach Sotiriou and Faraoni (2010).

ii.1 Two formulations of EiBI gravity

Let us now investigate the field equations of EiBI theory in more detail. Written in terms of the new rank-two auxiliary tensor , the theory is reminiscent of a particular bi-metric theory (see e.g. Refs. Drummond (2001); Milgrom (2009) for some specific attempt) which, in vacuum, degenerates into a single metric theory, GR. On the other hand, inside matter the two metrics are different and the connection is metric with respect to , and not , cf. Eq. (3). Nonetheless, the conservation of the stress-energy tensor – and geodesic motion therein – is guarantee by a minimal gravity-matter coupling in the action (2). This is precisely the essence of the theory.

It should be noted however that there are at least two ways of interpreting the theory to which we shall refer to, with some abuse, as the “Jordan frame” and the “Einstein frame”, in accordance to scalar-tensor theories. In the Jordan frame matter is minimally coupled to the metric and Einstein equations are modified according to Eq. (6). In the Einstein frame the metric is treated as fundamental and matter is non-minimally coupled to it. This appears quite naturally in the small limit of the theory, Eq. (6), where can be written solely in terms of and its derivatives by using Eq. (4).

Let us however start by pointing a particular symmetry of the action. In the Jordan frame, the action reads as in Eq. (2). In what we define as the Einstein frame, the action can be written explicitly by defining a new metric ,

 Sγ = −2λκ∫d4x(√|γab−κRab(Γ)|−1λ√−γ) (7) +SM[γab−κRab(Γ),ΨM].

Written in this form, the metric is now non-minimally coupled to matter. Furthermore, the field equations impose , but the connection is not anymore the -metric connection, since enters explicitly in the matter action through .

Neglecting the matter action in and , the following duality holds

 gab↔γab,κ↔−κ,λ↔1/λ,Γ↔Γ, (8)

up to an overall factor . This duality can be seen in the vacuum field equations, where the rescaling follows from the algebraic relation (5). Note that only in the case , i.e. for a vanishing cosmological constant, . This also shows that the case is special not only because it corresponds to asymptotically flat space in both frames, but it is a fixed point of the duality (8). If , negative values of in the Einstein frame correspond to positive values of in the Jordan frame, provided one modifies the matter coupling.

However, the coupling to matter breaks the symmetry above, introducing a non-minimal coupling in the Einstein frame. Interestingly, the action (7) in the Einstein frame reduces, in the small limit, to Einstein’s theory with non-minimally coupled matter (see also Ref. Bertolami et al. (2007)),

 Sγ = λ∫d4x√−γ(R−2λ−1κλ) (9) +∫d4x√−γ(1−12κR(Γ))Lm(γab−κRab(Γ)) +O(κ2).

where we have defined .

Furthermore, if the matter Lagrangian is constructed in terms of the metric , the second term in Eq. (9) involves other corrections and some higher order coupling to the Ricci tensor and to the independent connection. For example, let us consider the small limit of a minimally coupled free scalar field in the Jordan frame

 Lm(gab,ϕ)=gab∂aϕ∂bϕ. (10)

In the Einstein frame, the matter Lagrangian reads

 Lm(γab,Γ,ϕ)=γab∂aϕ∂bϕ−κRab(Γ)∂aϕ∂bϕ, (11)

leading to the following matter action in the Einstein frame

 Sm = ∫d4x√−γ[γab∂aϕ∂bϕ −κ(Rab+12Rγab)∂aϕ∂bϕ+O(κ2).

Note that the action above involves second derivatives of the scalar field in the gravity equations hidden inside the Ricci tensor evaluated on shell.

A similar transformation as the one presented above persists in the non-relativistic limit. Indeed, if we define , Eq. (1) reduces to the standard Poisson equation,

 ∇2ϕ=4πGρ. (13)

However, in this case the gravitational force would read , being thus dependent on the matter field.

The modified Poisson equation (1) has an interesting interpretation in terms of standard Poisson equation with a modified source term. The second term on the right hand side of Eq. (1) leads to an effective force:

 Feffρ≡−∇Peffρ=−κ4∇ρ, (14)

from which we obtain that . Therefore, the modified Poisson equation is formally equivalent to the standard one supplemented by an effective polytropic fluid with equation of state (EOS) , where the polytropic index and . In fact, a similar property exists also in the relativistic case. In the weak field limit of EiBI theory, one can show that an effective pressure and effective internal energy appears. This will be explained in detail in a forthcoming work Steinhoff and Delsate (2012).

ii.2 Linear structure of EiBI theory

The linear dynamics of EiBI theory is equivalent to linearized Einstein’s theory. To prove this, we expand the field equations (4)-(5) around a vacuum Minkowski background,

 qab=ηab+δqab,gab=ηab+δgab, (15) qab=ηab−δqab,gab=ηab−δgab, (16) Tab=δTab,Tab=δTab, (17)

where indeces are raised and lowered by the Minkowski metric, . In order for the Minkowski metric to be solution, we set in the field equations. Linearization of Eq. (4) follows immediatly from the standard linearized Ricci tensor,

 δqab−δgab=κδRab≡κ2(δqca,bc+δqcb,ac−□δqab−δqab), (18)

where a coma denotes partial derivative with respect to the flat background and is the trace of . In order to linearize Eq. (5), we recall that, at linear order, , so we get

 δqab−ηab2(δq−δg)=δgab+κδTab (19)

By taking the trace of the equation above, we get , which can be substituted back into Eq. (19). Finally Eq. (18) can be written as

 δRab=δTab−ηab2δT, (20)

which does not depend on . This equation has exactly the same form as in GR, but the Ricci tensor is written in terms of . However, at linear order, we note that

 δTab=δTabδΨM∣∣∣vacuumδΨM+δTabδglm∣∣∣vacuumδglm, (21)

where the second term vanishes when evaluated in vacuum, since matter fields are coupled to the metric. Therefore, and do not depend explicitly on the metric and Eq. (20) is exactly equivalent to the linearized Einstein equations to all orders in , but for the auxiliary metric .

In particular, the linearized version of EiBI theory admits the same differential structure as linearized GR. For example the principal symbols, i.e. terms in the linear equations involving highest order derivatives, are the same as GR. Although a detailed analysis of the well-posedness of the theory is beyond our scope, these results suggest that there exists a well-posed formulation of EiBI theory, similarly to GR. For a unitarity analysis of general BI gravity theories see Ref. Gullu et al. (2010).

Finally, in the small limit, our initial assumption of Minkowski background can be relaxed. In this case, similar conclusion about the linear structure of EiBI around any background metric can be drawn.

Iii Non-relativistic stellar collapse

In this section, we describe the non-relativistic collapse of dust in EiBI theory. Hereafter, we focus on asymptotically flat solutions, setting . The collapse of incoherent dust in the Newtonian limit shares many properties with its relativistic analogue Oppenheimer and Snyder (1939); Florides (1977). Here, we shall solve the hydrodynamics equations in the case of spherical symmetry, showing that the end-state of the evolution is the pressureless star found in Ref. Pani et al. (2011a) and described in detail in Section IV below. This confirms and extends the analytical argument presented in Ref. Pani et al. (2011a). However, at the end of this section we show that our analytical computation in Ref. Pani et al. (2011a) is partially flawed, due to a typo in the series expansion close to the origin, and we shall discuss why the numerical results are still in agreement with the overall picture given in Ref. Pani et al. (2011a).

The hydrodynamical equations can be equivalently solved using either the Eulerian or the Lagrangian approach. In the next two sections, we shall briefly review these formulations, together with a standard procedure to avoid shock wave formation, and the modifications related to the problem at hand.

iii.1 Eulerian formulation

In the non-relativistic limit, the collapse of non-interacting () particles is governed by the mass conservation (continuity equation) and momentum conservation (Euler equation), the latter being modified in EiBI gravity due to Eq. (1).

Following the seminal work VonNeumann and Richtmyer (1950), we supplement the hydrodynamics equations by an artificial viscosity term, in order to avoid divergences due to shock wave formation during the evolution. In the Eulerian formulation, the relevant equations in the pressureless case read Florides (1977); Iwakami et al. (2007)

 ∂ρ∂t+u⋅∇ρ+ρ∇⋅u = 0, (22) ∂u∂t+u⋅∇u+∇Φ+∇⋅Qρ = 0, (23) m−∫d3xρ = 0, (24)

where is the fluid velocity, is the density and is the mass function. Note that the second equation above is vectorial, is the viscosity tensor and is a vector.

The viscosity tensor is defined as Iwakami et al. (2007)

 Q=ℓ2ρ∇⋅u[∇u−e3∇⋅u] (25)

if , otherwise . In the equation above and is the unit tensor. Note that this equation only contains scalar invariant quantities so that can be directly specialized to the spherically symmetric case, contrarily to Eq. (8) in Ref. VonNeumann and Richtmyer (1950), which is only valid in Cartesian coordinates.

In indicial form,

 Qij=ℓ2ρ∇⋅u[∂jui+∂iuj2−δij3∇⋅u], (26)

where is a constant with dimension of length.

We are interested in the spherical symmetry case, where all the dynamical variable only depends on and the only non-vanishing component of vectorial quantities above is on the radial direction, e.g. . In this case only the radial component of Eq. (23) is nontrivial and , where . Due to the symmetry, the vector , and

 q(t,r)=ℓ2ρ(t,r)∇⋅u[u′(t,r)−13∇⋅u]. (27)

Finally, in spherical symmetry, and the artificial viscosity term in Eq. (23) is , so that we are left with the set of PDEs

 ∂ρ∂t+uρ′+ρu′+2rρu = 0, (28) ∂u∂t+uu′+Gmr2+κ4ρ′+Qρ = 0, (29) m′−4πr2ρ = 0, (30)

where a dot and a prime denote derivatives with respect to and , respectively. The viscosity term reads

 Q = 2ℓ23r2[rρ(u+2ru′)u′′ (31) −(u−ru′)(3ρu′+(2u+ru′)ρ′)],

if , otherwise . In practice, , where is the grid spacing in the radial direction and is some dimensionless constant. This allows the shock to be “smeared” in a region of width , while leaving the rest of the dynamics unchanged. In our simulations, we have checked that the results do not depend on the precise value of . Actually, the results do not even depend on the details of the artificial viscosity term, as long as it introduces a smearing effect on the shocks while keeping the rest of the physics unaltered. We have explicitly checked this fact, by comparing test simulations with different artificial terms.

Finally, one can derive a Bernoulli-like relation by multiplying Eq. (23) by , leading to

 ddt(12u(t,r)2+Φ)=∂Φ∂t. (32)

where we neglected the artificial viscosity term and where is the total time derivative. This expression is valid locally and the term on the right hand side takes into account the local variation due to matter flow. Indeed in the pure Newtonian case, this term can be rewritten as a surface integral of the density current.

iii.2 Lagrangian formulation

The basic Lagrangian equations are given in VonNeumann and Richtmyer (1950). They are expressed in terms of the comoving volume and the velocity field , where is the initial density. The continuity equation for the volume is given by

 ρ0(x)˙V(t,x)=˙X(t,x), (33)

which, in terms of the density , reads

 ˙ρ=−ρu′(t,x)X′(t,x)=−ρ∂u(t,x)∂X(t,x). (34)

If is a radial variable in spherical coordinates, this equation becomes

 ˙ρ(t,x)+2R(t,x)ρ(t,x)u+ρ(t,x)u′(t,x)R′(t,x)=0. (35)

Note that the equation above can be simply obtained from the corresponding Eulerian formulation by replacing the Eulerian time derivative by the Lagrangian time derivative according to

 DDt=∂∂t+u⋅∇, (36)

where is the time derivative in the Lagrangian formulation.

Following the same principle, the equation for the velocity field is given by

 ˙u=−Gm(t,x)R(t,x)2−κρ′(t,x)R′(t,x)+Fq(t,x)ρ0(x), (37)

where is the artificial viscosity term given by VonNeumann and Richtmyer (1950) which explicitly reads

 Fq=∂∂x(ρ0(x)(cΔx)2u′(t,x)|u′(t,x)|R′(t,x)), (38)

where is the grid spacing and is a constant.

iii.3 Results

We have solved the initial-value problem defined by the above system of PDEs using the method of lines, in which the radial dimension is discretized and the resulting system of coupled ODEs is integrated in time with standard methods. The initial static profiles we considered are presented in Table 1.

Profile I is a classical exponential density profile, while Profile II is a deformation of the pressureless star found in Ref. Pani et al. (2011a) (cf. Eq. (48) below). For both class of profiles, is a free parameter related to the slope of the initial density profile at the center.

For testing purposes, we wrote two independent codes that solve the hydrodynamical equations in the Eulerian and in the Lagrangian formulation, respectively. Once convergence is reached, the results of the two codes agree very well. However, the simulations in the Lagrangian formulation generically need a lower resolution, due to the fact that the radial grid evolves along with the fluid collapse. Indeed, the Lagrangian formulation is analog to the comoving frame generically adopted in relativistic collapse simulations. Hence, we shall present results obtained using this formulation. Typically, we used a non-uniform grid with points in the spacial direction, which guarantees convergence of the results for all values of taken into account. As expected, we find that, for a given value of the viscosity constant , convergence and stability of the numerical results is reached provided the mesh is sufficiently fine.

In standard Newtonian gravity, , the hydrodynamics equations can be solved analytically for a constant density profile and they correspond to the relativistic Oppenheimer-Snyder solution Florides (1977). In that case, the dust collapses in a finite time , where and are the initial radius and the total mass of the spherical dust configuration. This analytical results is also valid when non-constant initial density profiles are considered Pani et al. (2011a). Our simulations reproduce this result with very good precision. Furthermore, for any negative value of , we found a qualitatively similar behavior, i.e. the fluid collapses to a singular state in a finite time. For , the time of collapse decreases with .

On the other hand, the behavior for any is drastically different, as shown in Figs. 1 and 2. We present results for and obtained using Profile I in Table (1), but different choices of , different initial profiles or different values of , would give qualitatively similar results.

In Fig. 1 we show the central density as a function of time for two different values of . The central density is always finite and oscillates around a constant value at late time. As we shall show in Sec. (IV.1), the mean value of the oscillations (denoted by an horizontal dashed line in Fig. 1) is precisely the central density of a pressureless star Pani et al. (2011a) (cf. Sec. (IV.1)) with the same mass. In Fig. 2, we show the density and the velocity radial profiles at different instants. The would-be shocks in the velocity profile (the steep region shown in the right panel of Fig. 2) propagates towards the exterior of the fluid without developing a discontinuity, which could not be resolved due to finite grid spacing. The artificial term discussed above smears this discontinuity and ensures stability of the numerical simulations.

Our simulations suggest that, even for an arbitrarily small value of , this oscillatory behavior would continue indefinitely. This is perfectly consistent with the fact that the hydrodynamical equations do not include any dissipative term, so that energy is conserved during the evolution. Nevertheless, our results give strong evidences that, for generic (static) initial profiles, the end-point of the dust collapse is precisely the pressureless star. Indeed, we expect that any dissipative term, which has to be included in realistic situations, would quench the oscillations and drift the system toward a static configuration, which is precisely the pressureless star. This picture is also confirmed by Fig. 3, which shows the period of oscillation, as a function of , compared with the fundamental period of proper oscillation of a pressureless star (cf. Sec. IV.4).

Finally, as shown in Eq. (14), the pressureless collapse in the modified non-relativistic theory is equivalent to the collapse of a polytropic fluid with in Newtonian theory. Hence, our results are consistent with the fact that, in the latter case, the final state is a regular polytropic star with polytropic index (see e.g. Larson (1969) for more realistic simulations of star formation in Newtonian gravity).

iii.3.1 On the analytical method of Ref. Pani et al. (2011a)

In Ref. Pani et al. (2011a), we developed an approximated method to solve the collapse equation analytically for any . Here, we point out a typo in that computation and we discuss why this does not affect the final result, as we have proved in the previous section by solving the collapse equation numerically.

Let us review the procedure adopted in Ref. Pani et al. (2011a). Due to the spherical symmetry, we expect the singularity to initially form at the center . The collapse equations are then solved by expanding the dynamical variables close to the center,

 ρ = ρ0(t)+ρ1(t)r+ρ2(t)r2+O(r3), u = u0(t)+u1(t)r+u2(t)r2+u3(t)r3+O(r4),

and solving a system of ODEs for and . At first order, one finds and

 u1(t)=−˙ρ0(t)3ρ0(t). (39)

At second order, we get

 ¨ρ0(t)ρ0(t)−43˙ρ0(t)2ρ0(t)2−32κρ2(t)−8πGρ0(t) = 0 (40) ˙ρ2(t)ρ2(t)−53˙ρ0(t)ρ0(t)+5ρ0(t)ρ2(t)u3(t) = 0 (41)

The results we presented in Ref. Pani et al. (2011a) originated by erroneously omitting the term proportional to in Eq. (41). Indeed, if , then , which can be substituted in Eq. (40) to obtain a single non-linear ODE for . The solution of that equation is given in the equation below Eq. (7) of Ref. Pani et al. (2011a).

However, in general is non-vanishing and Eqs. (40) and (41) are not sufficient to solve for the three variables , and . In fact, it can be easily proved that, at any order in the series expansion, the number of unknown functions is always larger than the number of differential equations. This is precisely due to the extra derivative term in the Poisson equation (1). Indeed, in the standard case when , Eq. (40) decouples and can be solved in the usual way. As a matter of fact, when , the collapse equations cannot be solved exactly by this simple method, as erroneously stated in Ref. Pani et al. (2011a).

Nevertheless, as we have showed above, our fully numerical simulations not only confirm the results previously reported, i.e. that the collapse does not lead to any singularity, but also that the final state is the pressureless solution reported in Ref. Pani et al. (2011a). Furthermore, the oscillatory behavior found in the numerical simulation qualitatively match the oscillatory behavior presented in Fig. 1 of Ref. Pani et al. (2011a).

Thus, it is relevant to understand whether the analytical method present above, although flawed, can actually reproduce the full solution in some limit. This is certainly true when at any time. The latter condition is indeed satisfied when the velocity and density profiles are close enough to a late-time evolving configuration, for example when the configuration oscillates around the static configuration. Indeed in this case, while , where are the constant central density and second derivative at the center of the static pressureless configuration and is a small number. It follows that in this case, the combination is one order smaller than .

Iv Stars in the non-relativistic limit

In this section, we discuss non-relativistic stellar models in the theory defined by (2). Our results apply to any relativistic theory which reduces to Eq. (1) in the non-relativistic regime. We shall solve the hydrostatic equilibrium equation, supplemented by the standard mass conservation, , and an EOS. Requiring spherical symmetry, the hydrostatic equilibrium equation reads Pani et al. (2011a)

 dPdr=−Gm(r)ρr2−κ4ρρ′. (42)

Clearly, equilibrium configurations in EiBI gravity are different from the standard Newtonian ones in presence of non-trivial density profiles. Note that the equation above can be written in a more evocative form as Casanellas et al. (2011)

 dPdr=−Geff(r)m(r)ρ(r)r2, (43)

where we have defined an “effective” Newton’s constant

 Geff(r)≡G+κ4r2ρ′(r)m(r). (44)

Since inside a star, when . In Ref. Casanellas et al. (2011), the modified hydrostatic equilibrium equation has been used to construct realistic solar models and put the first observational constraint on (see also the recent Ref. Avelino (2012) for more stringent constraints).

Here, we shall mainly focus on polytropic models with EOS of the form . Using the mass conservation equation, we can rewrite Eq. (42) as

 1r2[r2ρρ′dP(ρ)dρ]′+4πGρ+κ4(ρ′′+2ρ′r)=0, (45)

where a prime denotes derivative with respect to . Once the EOS is fixed, Eq. (45) is a second order ODE for , which can be solved by imposing regularity condition at the center, i.e.

 ρ(r)=∞∑n=0ρnrn,r→0. (46)

In the equation above, are constants which can all be expressed in terms of the only free parameter, , by solving Eq. (45) perturbatively close to the center. For example,

 ρ2=−8Gπρ203(4P′(ρ0)+κρ0), (47)

whereas for any odd . It is easy to see that stellar models only exist when , which implies a lower bound, . As we shall see, this limit holds also qualitatively in the relativistic limit and in more realistic models Casanellas et al. (2011).

iv.1 Non-relativistic pressureless stars

In classical Newtonian gravity, a gas of non-interacting particle cannot support self-graviting configurations and it will inevitably collapse due to its self-gravity attraction. However, in theories in which the Poisson equation reads as in Eq. (1), a straightforward solution of the hydrostatic equilibrium equation with and reads Pani et al. (2011a)

 ρ(r)=ρcsin(ϖr)ϖr,ϖ=4√Gπκ. (48)

 R=πϖ,M=4π2ρcϖ3, (49)

respectively, so that for a given , the solution above is uniquely characterized by the central density or, equivalently, by the mass. In particular, the radius of the compact object is uniquely determined by .

The existence of such solutions is remarkable not only because they are interesting models for self-gravitating dark matter Bertone et al. (2005), but also because they provide a non-trivial static and spherically symmetric solution of the hydrodynamics equations governing the collapse of dust. Indeed, as we discussed in Sect. III, our numerical simulation give strong evidences that, as an outcome of the collapse, any initial distribution of matter will eventually relax towards the pressureless solution described by Eq. (48).

Clearly, the generality of this result is only due to the fact that we considered the collapse of pressureless matter, but it also suggests that, when taking into account realistic models in which the matter is described by a specific EOS, the end-state of the collapse will be a modified compact star, rather than a singular state. A detailed analysis in that direction would be certainly interesting. For the time being, it is important to notice that these solutions are also linearly stable, as first shown in Ref. Pani et al. (2011a) and detailed in Sec. IV.4.

We note here that the pressureless star (48) is a sort of solitonic solution, since it solves the Poisson equation for any harmonic function , due to the identically vanishing of the right hand side of Eq. (1), i.e. . In the interior, the Newtonian potential is constant and it matches continuously the vacuum potential at the radius. Indeed, Eq. (1) can be thought as a forced oscillator equation for , where would play the role of the inverse spring constant and is an external force. The pressureless star qualitatively corresponds to the solution of the oscillator with a constant external force.

iv.2 Newtonian polytropic models

For a generic polytropic index , the field equation must be solved numerically, imposing at the center. In this case stellar solutions only exist provided the following condition is satisfied

 κ>−|κc|=−4K(1+1/n)ρ−1+1/nc. (50)

For , condition (50) is always fulfilled. In some cases the Lane-Emden equation obtained from Eq. (42) can be solved analytically Chandrasekhar (1957). For instance if , , and the solution reads as in (48), but with , so that it exists for and reduces to the pressureless case for . This is consistent with the effective pressure term defined in (14). Indeed, the modified pressureless star solution (48) with is exactly the same as the polytropic solution with .

iv.3 Modified Chandrasekhar model for white dwarfs

As a relevant example of non-relativistic polytropic star, let us consider zero temperature non-rotating white dwarfs. The interior of a zero-temperature white dwarf is described by a relativistic EOS, , which is parametrically defined by Chandrasekhar (1935)

 P(x) = πm4ec5μemP3h3[x(2x2−3)√1+x2+3sinh−1x], ρ(x) = 8πm3ec3μemP3h3x3, (51)

where and are the speed of light and the Planck constant, and are the electron and the proton mass, respectively and is the molecular weight. In Fig. 4 we show the mass-radius relation for different values of obtained by integrating numerically the hydrostatic equations.

When we recover the famous Chandrasekhar result: zero temperature white dwarfs have a maximum mass, . However, for any , the behavior is drastically different. White dwarf models can have an arbitrarily large mass, but instead there exists a minimum radius, which approximately scales as .

We can use Landau’s original argument to understand these results. Following Shapiro and Teukolsky’s exposition Shapiro and Teukolsky (1983), let us consider a star of radius composed of fermions, each of mass . In the relativistic regime, the Fermi energy of the system reads

 EF=ℏcN1/3R. (52)

In EiBI theory, the gravitational energy per fermion is approximately

 EG≈−GNm2bR+3κNm2b16πR3, (53)

where we approximated the star’s density as . Thus, the star’s total energy read

 E≡EF+EG=ℏcN1/3−GNm2bR+3κNm2b16πR3, (54)

For small , the total energy is positive, and we can decrease it by increasing . If moreover

 R>Rcrit≡√3κ16πG, (55)

then the total energy can become negative. At some point the electrons become non-relativistic and the purely Newtonian term dominates. It is easy to see that when this happens the energy tends to zero when goes to infinity. Thus, there is a point of minimum energy, and it is delimited by the curve (55). Note that the scaling is consistent with the recent computation of the Jeans length in Ref. Avelino (2012).

On the other hand, even if initially, for sufficiently small condition (55) will be violated: the EiBI term dominates at sufficiently small radii, and collapse is never energetically favored. Relativistic stellar models of white dwarfs are described in Sec. V.4.

iv.4 Stability in the non-relativistic limit

We now discuss stability of the Newtonian configurations against radial perturbations Shapiro and Teukolsky (1983), expanding the discussion in Ref. Pani et al. (2011a). The equations governing the stellar dynamics are the Poisson equation (1) together with the continuity equation and the momentum equation. These are given in Eqs. (22)-(24) by replacing the artificial viscosity term by a physical pressure term, i.e. . In those equations, the hydrostatic equilibrium is recovered when .

We shall focus on spherically symmetric models. The perturbed Euler equation reads

 Δ[dudt+P′ρ+Φ′]=0. (56)

where is the total time derivative, and is the Lagrangian displacement. In order to compute this equation explicitly we use the following Shapiro and Teukolsky (1983)

 Δρ = −ρr2(r2ξ)′, δρ = −(r2ρξ)′r2, (δΦ)′ = −4πρξ+κ4(δρ)′=−4πρξ−κ4[(r2ρξ)′r2]′, ΔP ≡ PγΔρρ,

where the last equation above defines the adiabatic index of the perturbations, . It is then straightforward to obtain the following

 Δ(dudt) = d2ξdt2, Δ(P′ρ) = −Δρρ2P′+Δ(P′)ρ=2ξrρP′+(ΔP)′ρ = 2ξrρP′−1ρ[Pγ(r2ξ)′r2]′, Δ(Φ′) = (ΔΦ)′−ξ′Φ′=(δΦ)′+ξΦ′′ = (δΦ)′+ξ∇2Φ−2rξΦ′ = 2ξρrP′+κ4[2rξρ′−ξ′ρ′−[ρr2(r2ξ)′]′],

where we have used and the background equations. Finally, using the equations above, Eq. (56) reads

 ¨ξ−1ρ[γPr2(r2ξ)′]′+4ρrξP′ +κ4[2rξρ′−ξ′ρ′−[ρr2(r2ξ)′]′]=0. (57)

Assuming a time dependence , the modified eigenvalue equation reads Pani et al. (2011a)

 4ξP′r+κρ4[2rξρ′−ξ′ρ′−[ρr2(r2ξ)′]′]−[γPr2(r2ξ)′]′=ρξω2. (58)

This equation must be solved for , requiring the following boundary conditions Shapiro and Teukolsky (1983)

 ξ(0)=0,ΔP(R)=−γPr2(r2ξ)′=0, (59)

the latter being equivalent to requiring regularity of at the radius. An instability corresponds to an eigenmode with .

iv.4.1 Perturbations of pressureless stars

For and given by Eq. (48), the eigenvalue equation (58) simplifies. In particular it does not depend on the index of perturbations , but only on . Thus, for a given , there is one fundamental mode. By integrating Eq. (58) numerically and imposing Eqs. (59), we find

 ω=αρ1/2c, (60)

where , independently from .

This gives the characteristic period of oscillation

 ΔTM=π5/42α(κM2)3/4, (61)

which is shown in Fig. 3 and compared with the period of oscillation of our numerical solutions around a pressureless star with same mass. The results of the linear analysis perfectly agrees with the simulations, giving further support for the end-state of the collapse.

In addition, we found no unstable modes of the pressureless star, confirming that in the modified Newtonian theory these object are linearly stable Pani et al. (2011a).

iv.4.2 Newtonian polytropic stars

In Newtonian gravity, polytropic models with are marginally stable for any polytropic index  Shapiro and Teukolsky (1983). In our case, these models are stable if and unstable if . A representative example is shown in Fig. 5, where we show the fundamental mode as a function of for the polytropic model with and . For generic values of , positive values of contribute to stabilize the models, while negative values work in the opposite direction.

V Relativistic compact stars

In this section, we study compact objects in the fully relativistic theory defined by the action (2). Let us now consider static and spherically symmetric perfect fluid stars. The metric ansatzen read

where we have used the gauge freedom to fix the function in front of the spherical part of the auxiliary metric .

Notice that the field equations (5) are simply algebraic equations relating and . Inserting the ansatzen above into Eqs. (5), and solving for the coefficients of in terms of , leads to

 F = p(r)[1+κρ]2A3(r), (64) B = h(r)A(r), (65) A = [(1−κP)(1+κρ)]−1/2. (66)

which correctly imply in vacuum. The dynamical field equations (4) read

 F = p−κ(p′(rph′+h(−4p+rp′))−2rhpp′′)4rh2p, (67) B = rFh2−κ(ph′+hp′)rhp, (68) A = 1−κ4r2(4+2rh′h2−4h−2rp′hp), (69)

where , and are given in Eqs. (64)-(66). The equations above, supplied by an EOS in the form , can be solved for , and . Since matter is covariantly coupled with the metric , the standard conservation of the stress-energy tensor follows, , where we recall that the covariant derivative is defined in terms of the physical metric, . We consider perfect-fluid stars with energy density and pressure such that

 Tab≡Tabperfectfluid=[ρ+P]uaub+gabP, (70)

where the fluid four-velocity . Thus, rather than using Eq. (69), it is more convenient to use the conservation of the stress-energy tensor,

 2FP′+F′(P+ρ)=0. (71)

Notice that, although the equation above simply reads as in GR, introduces terms proportional to through Eq. (64). Finally, in the limit , the field equations reduce to that of GR with .

v.1 Integration

We have integrated Eqs. (67), (68) and (71) imposing regularity conditions at the center of the star, where the following expansions hold

 p(r)∼p0+p2r2,h(r)∼1+h2r2, P(r)∼Pc+P2r2,ρ(r)∼ρc+ρ2r2.

We can set by a time reparametrization. Furthermore, using the field equations, the coefficients can be written in terms of and as follows

 P2 = 2(1−Pcκ)(Pc+ρc)(1+κρc)((1+κρc)2A3c−1)3κΔ, Δ = (Pcκ−3κρc−4)(1+κρc)−κ(1−κPc)(Pc+ρc)ρ′c, p2 = (1+