Anisotropic stars as ultracompact objects in General Relativity

Anisotropic stresses are ubiquitous in nature, but their modeling in General Relativity is poorly understood and frame-dependent. We introduce the first study on the dynamical properties of anisotropic self-gravitating fluids in a covariant framework. Our description is particularly useful in the context of tests of the black hole paradigm, wherein ultracompact objects are used as black hole mimickers but otherwise lack a proper theoretical framework. We show that: (i) anisotropic stars can be as compact and as massive as black holes, even for very small anisotropies; (ii) the nonlinear dynamics of the 1+1 system is in good agreement with linearized calculations, and shows that configurations below the maximum mass are nonlinearly stable; (iii) strongly anisotropic stars have vanishing tidal Love numbers in the black-hole limit; (iv) their formation will usually be accompanied by gravitational-wave echoes at late times.


Introduction.
A foundational result in General Relativity (GR) states that the maximum compactness of a self-gravitating, isotropic, spherically-symmetric object of mass M and radius R is M/R = 4/9, if the object is composed of a perfect fluid [1] (we use G = c = 1 units). As a corollary, under the assumptions above, the existence of exotic compact objects (ECOs) of compactness arbitrarily close to that of a Schwarzschild black hole (BH, with M/R = 1/2) is ruled out. Thus, tests of the BH paradigm -are the dark and massive objects that we see really BHs? -are challenging to devise, impacting our ability to quantify statements about evidence for BHs [2][3][4][5], or to discover new species of compact objects.
It has been realized that Buchdahl's bound above relies strongly on the hypothesis of isotropy. Anisotropies in matter fields arise naturally at high densities [6][7][8] and may play an important role in the interior of compact objects. The simplest known example is that of a scalar field minimally coupled to gravity, which indeed gives rise to anisotropic pressure in boson stars [9]. Other examples include electromagnetic fields, fermionic fields, pion condensed phase configurations in neutron stars [10], superfluidity [11], solid cores [6], etc. In the real world, anisotropic pressures are the rule rather than the exception.
Surprisingly, anisotropic stars in GR are poorly studied. While various solutions have been obtained, both in closed form [12][13][14][15][16][17][18][19] and numerically [20][21][22][23][24][25], none arises from a consistent covariant model (see Ref. [26] for some progress). The lack of a proper framework prevents the exploration of outstanding questions associated to these objects, such as their stability, dynamical formation, and phenomenology. This is in stark contrast with the excellent knowledge on the dynamics of BHs and neutron stars, and is also the most important limitation in the study of ECOs [2,3] (the only exception being boson stars which, however, are even less compact than the Buchdahl's bound [9] and do not belong to the ClePhOs category introduced in Refs. [2,3]).
Here we introduce a covariant and self-consistent model for spherically-symmetric anisotropic fluids in GR, admitting stable and well-behaved ultracompact solutions which we term C-stars. Covariant approach to anisotropies. Consider an anisotropic fluid with radial pressure P r , tangential pressure P t , and total energy density ρ, described by the stress-energy tensor where u µ is the fluid four-velocity and k µ is a unit spacelike vector orthogonal to u µ , i.e. k µ k µ = 1 = −u µ u µ , u µ k µ = 0. Here, Π µν = g µν + u µ u ν − k µ k ν is a projection operator onto a two-surface orthogonal to u µ and k µ , i.e., u µ Π µν V ν = k µ Π µν V ν = 0 for any vector V µ . At the center of symmetry of the fluid the anisotropy P r − P t must vanish [12]. There is a certain degree of arbitrariness to satisfy this condition in a covariant fashion, the simplest possibility is where f (ρ) is a generic function of the density and the free constant C is a (generically dimensionful) parameter that measures the deviation from isotropy. By construction, P t = P r at the center of static and sphericallysymmetric objects, since ∂ r P r | r=0 = 0. By defining σ := f (ρ)k µ ∇ µ P r , we can write Eq. (1) as the stress-energy tensor of an isotropic perfect fluid plus an anisotropic contribution, In the spherically-symmetric case, u µ = (u 0 , u 1 , 0, 0), k µ = (k 0 , k 1 , 0, 0), and all dynamical variables are func-arXiv:1811.07917v1 [gr-qc] 19 Nov 2018 tions of (t, r) only. The orthogonality conditions provide two constraints on k µ , which is therefore completely fixed in terms of u µ . It is straightforward to show that Π µ ν = diag(0, 0, 1, 1), which simplifies some of the computations presented below. C-stars: equilibrium configurations. For static solutions, the metric can be written as ds 2 = −e ν(r) dt 2 + (1 − 2m(r)/r) −1 dr 2 + r 2 dθ 2 + sin 2 θdϕ 2 . The metric variables satisfy the standard relation, m (r) = 4πr 2 ρ and ν (r) = 2(m + 4πr 3 P r )/(r(r − 2m)), whereas the radial pressure satisfies a modified Tolman-Oppenheimer-Volkoff equation which reduces to the isotropic case when C = 0. Two equations of state for P r and f (ρ) are necessary to close the system. The simplest choice for the function f would be f (ρ) = 1, but in this model P t is discontinuous at the stellar radius since P t (R) = 0; see Eq. (2). The simplest model that ensures continuity of P t and its derivative at the radius is f (ρ) = ρ. We focus on this case here, although other models (e.g. f (ρ) = ρ n ) give similar results. With this choice, Eqs. (4) and (2) guarantee that P r = P r = P t = P t = 0 at r = R.
Remarkably, Eq. (4) can be solved in closed form for a toy model of incompressible fluid (ρ(r) = const), although the solution is cumbersome. We focus instead on a standard polytropic equation of state P r = Kρ γ 0 with adiabatic index γ = 2, where ρ 0 = ρ − P r /(γ − 1) is the rest-mass density. Our results are qualitatively the same for other standard equations of state.
The mass-radius diagram of C-stars is shown in Fig. 1 for different values of C. There are several features worth highlighting: (i) overall, C-stars can be much more compact and massive than isotropic stars and their maximum compactness always approaches that of a Schwarzschild BH in some region; (ii) More importantly, C-stars exist across a wide range of masses, evading one of the outstanding issues with BH mimickers: in most theories giving rise to ECOs, these approach the BH compactness in a very limited range of masses, thus being unable to describe both stellar-mass and supermassive BH candidates across several orders of magnitude in mass. On the other hand, C-stars can do so when C/M 3 1. (iii) As shown in the inset of Fig. 1, compact configurations exceed the Buchdahl's limit and can even classify as ClePhOs [2,3]. (iv) When C ≥ 0 the fluid has P t > 0 everywhere inside the star, and satisfies the weak and the strong energy conditions [27] (ρ+P r +2P t ≥ 0, ρ+P r ≥ 0, and ρ+P t ≥ 0), whereas very compact configurations violate the dominant energy condition (ρ ≥ P r and ρ ≥ P t ) near the radius, where P t attains a maximum and P t > ρ. This maximum is very sharp and approaches the radius of the star for very compact configurations, in a way reminis- cent of gravastars with a thin layer of strongly anisotropic pressure [28][29][30][31]. Radial stability. For any C > 0 the compactness and the anisotropy grow in the high-density region, eventually reaching the BH compactness (see inset of Fig. 1). Thus, even a vanishingly small value of anisotropy parameter C can give rise to strongly-anisotropic quasi-Schwarzschild equilibrium solutions. When C is small standard analysis of the turning points in the mass-radius diagram [32] suggests that these configurations are unstable. On the other hand, in the strong-anisotropy regime, the massradius relation of a C-star approaches that of a BH already on the stable branch. To test these issues, we perform a linear stability analysis of C-stars under radial perturbations. The spacetime metric is written as µν is the metric of a background C-star solution and h µν = diag (H 0 (r), H 2 (r), 0, 0) e −iωt is a small perturbation in Fourier space. Likewise, we expand the fluid density, pressure, and vector components u 0,1 and k 0,1 as X = X 0 + δXe −iωt , where X 0 collectively denotes the background quantities and δX is the corresponding radial perturbation. The orthogonality conditions on u µ and k µ can be used to relate δu 0 , δk 0 and δk 1 to the remaining functions. Radial fluid perturbations propagate at the speed c s = ∂P r /∂ρ, which is always real and subluminal for these configurations. On the other hand, the tangential speed of sound can only be computed for nonspherical perturbations. The linear system can be reduced to a second-order differential equation for the fluid displacement, ξ(r) = i u r ω e ν/2 . The eigenvalue problem is solved by requiring ξ(0) = 0 and ∆P r (R) = 0, where ∆P r = δp + ξ∂ r P r is the Lagrangian variation of the pressure [33]. This selects a discrete set of frequencies ω 2 , with ω 2 > 0 (ω 2 < 0) defining stable (unstable) modes.
Our results are summarized in Fig. 2, where we show the fundamental modes as a function of the compactness for representative values of C. All the expectations based on the mass-radius diagrams are confirmed: configurations with central density below (above) that corresponding to the maximum mass are linearly stable (unstable). Strongly-anisotropic configurations are linearly stable for M/R 0.42, while they become linearly unstable for higher values of the compactness. Configurations on the left of the cusps (corresponding to the zero crossing of ω 2 ) are linearly stable, whereas those on the right are linearly unstable. The threshold corresponds to the maximum mass of the object shown in Fig. 1. We also show the echo delay time (6) for these configurations (dashed curves). The markers refer to the time scale extracted from the nonlinear evolutions, which systematically predict more stable configurations.
Evolution equations of anisotropic fluids. We now discuss the full nonlinear theory. The fluid dynamics is governed by the conservation of energy-momentum and baryonic number. Notice that, since Π µ ν = diag(0, 0, 1, 1), the anisotropies only appear in the spatial projection of the stress-energy tensor. One of our main results is that the systems of partial differential equations describing the fluid and the dynamical spacetime, which are de- scribed in detail in the Supplement Material, is well behaved, and fully nonlinear simulations can be performed. Our simulations confirm the stability properties of the equilibrium configurations found in the previous section. Figure 3 displays the evolution of the central value of the rest-mass density of both stable and unstable equilibrium configurations, for different values of the parameter C, as a function of time. Clearly, small numerical perturbations drive unstable solutions away from their original configuration, whereas they remain bound for solutions in the stable branch. The full nonlinear results for the timescale τ are compared with the linear analysis in Fig. 2. We find good agreement for stable configurations at moderately small compactness. For large compactness and large values of C, the nonlinear evolution shows that the threshold for stability is pushed to larger compactness, as compared to linear analysis. Furthermore, although not shown in Fig. 2, we find indication that unstable configurations typically have a lifetime longer than the one predicted solely by linearized studies. We do not have a solid explanation for this discrepancy; it could be due to nonlinearities driving the energy to higher modes, or to other effects. We postpone a detailed analysis for the future, but we point out that such nonlinear results are potentially exciting: the merger of two C-stars might form an ultracompact configuration which lies on the unstable branch, but with long lifetimes and therefore with the potential to impart unique signatures on the postmerger GW signal, some described below. In addition, C-stars with compactness M/R ≈ 0.4 are nonlinearly stable. To the best of our knowledge, this is the first model of ultracompact objects featuring a clear photon-sphere (at R = 3M ) and which are dynamically well-behaved. Tidal Love numbers. The tidal Love numbers (TLNs) define the deformability of a star immersed in an external field, such as the one produced by a companion in a binary [34]. These quantities are particularly useful for GW astronomy, since they affect the late-inspiral GW signal from a coalescence and contain information about the nature of the merging objects [35]. The prime motivation to measure TLNs is to constrain the neutron-star equation of state [35,36] and to convey information on the nature of compact objects [37][38][39][40]: in GR, the TLNs of a BH are precisely zero [41][42][43][44][45][46], but are nonvanishing for ECOs [37,[47][48][49], being thus a smoking gun for ultracompact horizonless objects.
TLNs can be computed with standard techniques [34,[37][38][39]50]. As a proof of principle, we focus on the quadrupolar scalar TLNs [37], which are expected to be qualitatively similar to the gravitational case but avoid the study of nonspherical gravitational perturbations. In the large-compactness limit, our results are consistent with the relation where k scalar 2 is the scalar TLN as defined in Ref. [37], ∆ is the proper distance [51] between R and the Schwarzschild radius 2M , and a ∼ O(1), p ≈ 1.2, and n ≈ (3 − 3.5) mildly depend on C. Remarkably, this behavior is markedly different from that of other ECO models, for which the TLNs vanish logarithmically, k 2 ∼ log(∆/M ), in the BH limit [37,39], and shows that the TLNs of Cstars are very small as M/R → 1/2. As reference, for a neutron star k scalar 2 ≈ k gravitational 2 ≈ 200 [36].
Smoking guns: echoes. GW echoes in the post-merger GW signal from a binary coalescence are a smoking gun for structure at the horizon scale [52][53][54][55]. Our scope here is to simply show that perturbed C-stars produce echoes when sufficiently compact, a more detailed analysis is left for the future. We consider a test free scalar field on the background of a C-star. Figure 4 shows the linear response of a C-star with initial condition ∂ t ψ(x, 0) = exp(−(x − x 0 ) 2 /σ 2 ) and ψ(x, 0) = 0. The time delay between echoes roughly corresponds to twice the light crossing time from the center of the star to the photon sphere [52,53,55], Interestingly, this delay time is typically dominated by the the travel time within the star, not by the the Shapiro delay factor ∼ log(1 − 2M/R) near the surface [53]. This property is akin to (isotropic) ultracompact stars near the Buchdahl's limit [3,55] and to certain phenomenological models [57] considered in the past. Discussion. To summarize, we introduced a covariant framework to study anisotropic stars in GR, and have shown that the resulting system of equations is wellposed. Thus, relativistic anisotropic fluids can be explored in full-blown nonlinear evolutions, as we hope to do in the near future. Our results are exciting: C-stars provide a prototypical model for ultracompact objects, of immense utility in the quest to quantify the evidence for BHs and in identifying possible smoking guns for new physics [2,3,5]. Such configurations can be metastable and display the whole phenomenology recently predicted for ultracompact horizonless objects. In particular, they can be as massive and compact as BHs, have vanishingly small TLNs [37], and produce GW echoes [52,53] when perturbed (evading recent results [58], due to anisotropy).
It is intriguing to notice that C-stars share many key properties with gravastars [28,30], although the dynamics of the latter lacks a solid theoretical framework (see [59] for recent progress). Thus, C-stars might serve as an effective model for semiclassical corrections near the horizon, as predicted in other contexts [31,60]. Some of these models are nonperturbative in the Planck length P , and they would predict C/M 3 ∼ M/ P ∼ 10 38 for M = M , which motivates the strong-anisotropy regime explored here. It is also likely that the magnitude of anisotropies grows with the compactness of the object. Anisotropic effects might become stronger during the merger and an ordinary neutron star might "anisotropize" dynamically.
We have worked with a very crude toy model. Generalizations include, for example, models with f (ρ) = ρ n ; preliminary results for n > 1 show that stars in the stable branch are even more compact than the models presented here. We are tempted to conclude that there are very generic models of anisotropy which lead to the same phenomenology as the one we report.
Finally, we focused on the nonlinear dynamics in the spherically symmetric case; extensions of our covariant formalism to less symmetric configurations and on simulations of binary C-stars are ongoing. This is particularly interesting in light of our results: the mass-radius diagram of C-stars suggests that, for (say)C = 10 8 , two merging C-stars with equal mass M ≈ 11 M (compactness M/R ≈ 0.18) might give rise to a stable C-star near the maximum mass, M f ≈ 21 M and with compactness M f /R ≈ 0.43, being thus a viable candidate for an ECO + ECO → ECO coalescence. verse" and during the "Fundamental Physics with LISA" workshop, respectively. PP acknowledges financial support provided under the European Union's H2020 ERC, Starting Grant agreement no. DarkGRA-757480 and the kind hospitality of the Universitat de les Illes Balears, where part of this work has been performed. CP and MB acknowledge support from the Spanish Ministry of Economy and Competitiveness grants AYA2016-80289-P (AEI/FEDER, UE) and AYA2017-82089-ERC. CP also acknowledges support from the Spanish Ministry of Education and Science through a Ramon y Cajal grant. MB would like to thank CONICYT Becas Chile (Concurso Becas de Doctorado en el Extranjero) for financial support. MB also acknowledge the support of the PHAROS COST Action (CA16214). V.C. acknowledges financial support provided under the European Union's H2020 ERC Consolidator Grant "Matter and strong-field gravity: New frontiers in Einstein's theory" grant agreement no. MaGRaTh-646597. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 690904. We acknowledge financial support provided by FCT/Portugal through grant PTDC/MAT-APL/30043/2017. We acknowledge the SDSC Comet and TACC Stampede2 clusters through NSF-XSEDE Award Nos. PHY-090003, as well as MareNostrum and the technical support provided by Barcelona Supercomputing Center (AECT-2018-1-0003). The authors would like to acknowledge networking support by the GWverse COST Action CA16104, "Black holes, gravitational waves and fundamental physics." We acknowledge support from the Amaldi Research Center funded by the MIUR program "Dipartimento di Eccellenza" (CUP: B81I18001170001).

Appendix A: Supplement Material: Evolution equations of anisotropic fluids
The conservation of stress-energy momentum and the conservation of baryonic current implies, respectively, the conservation of energy and momentum density and the conservation of mass that govern the fluid equations.
To write this covariant conservation law as an evolution system in spherical symmetric, one needs to split the spacetime tensors and equations into their space and time components by means of the 1 + 1 decomposition. The line element can be decomposed as where α is the lapse function, g rr and g θθ are positive metric functions, and dΩ 2 = dθ 2 + sin 2 θdϕ 2 the solid angle element. These quantities are defined on each spatial foliation Σ t with normal n a = (−α, 0) and extrinsic curvature K ij ≡ − 1 2 L n γ ij , where L n is the Lie derivative along n a .
Notice that, since Π µ ν = diag(0, 0, 1, 1) in spherical symmetry, the anisotropy function σ enters only in T θ θ and T φ φ , being the rest of T µ ν formally the same as for an isotropic fluid. The projections of this tensor and the baryonic current, in spherical symmetry, are given by D = ρ 0 W , U = hW 2 − P r , S r = hW 2 v r , (A3) S r r = hW 2 v r v r + P r , S θ θ = P r − Cσ, where we have defined the enthalpy h = ρ + p = ρ 0 (1 + )+p in terms of the rest mass density ρ 0 and the internal energy . Furthermore, we have defined ∂ t P r = f (α, g rr , u r , ∂u r , ∂ r P r , ∂ r ρ, σ;C) , It is straightforward to obtain generic evolution equations, in the sense that they do not depend on the specific form of the stress-energy tensor, for these projected quantities by projecting the conserved equation A1. The evolved conserved quantities {D, U, S r } are not modified by the anisotropies. Therefore, the algorithm to convert from conserved to primitive or physical fields {ρ 0 , , P r , v r }, given an equation of state P r = P r (ρ 0 , ), is the same as for isotropic fluids.
Einstein's equations can be written by using the Z3 formulation in spherical symmetry [61]. This formulation introduces independent variables in order to form a first order evolution system. The final system must be complemented with gauge conditions for the lapse. We use the harmonic slicing condition ∂ t ln α = −α trK, where trK = K r r + 2K θ θ . Finally, the evolution system is written in balance law form [62] where U = {α, g rr , g θθ , K r r , K θ θ , A r , D r rr , D r rθ , Z r , D, U, S r }, is a vector containing the final set of evolution field.