Embedding black holes and other inhomogeneities in the universe in various theories of gravity: a short review

The classic black hole mechanics and thermodynamics are formulated for stationary black holes with event horizons. Alternative theories of gravity of interest for cosmology contain a built-in time-dependent cosmological"constant"and black holes are not stationary. Realistic black holes are anyway dynamical because they interact with astrophysical environments or, at a more fundamental level, because of backreaction by Hawking radiation. In these situations the teleological concept of event horizon fails and apparent or trapping horizons are used instead. Even as toy models, black holes embedded in cosmological"backgrounds"and other inhomogeneous universes constitute an interesting class of solutions of various theories of gravity. We discuss the known phenomenology of apparent and trapping horizons in these geometries, focusing on spherically symmetric inhomogeneous universes.


Introduction
Black holes are a fundamental prediction of General Relativity (GR) and have been the subject of a celebrated theory of dynamics and thermodynamics developed in the 1970s [133]. The black holes considered in these studies are stationary solutions of GR. Real black holes, however, are non-stationary due to several possible processes: • A cosmological background: real black holes are embedded in the universe and, although this feature may be irrelevant for their astrophysics, their asymptotics are quite relevant in problems of principle and in mathematical physics. Alternative theories of gravity advocated to explain the present acceleration of the universe without an ad hoc dark energy contain a built-in, time-dependent effective cosmological "constant" and black holes in these theories are asymptotically Friedmann-Lemaître-Robertson-Walker (FLRW).
• An astrophysical environment, which could include a companion in a binary system (this is the situation which led to the detection of gravitational waves in a black hole merger by LIGO [1,2,3]), an accretion disk, or spherical accretion.
• Hawking radiation and evaporation affect all black holes and are, therefore, important and unavoidable in all problems of principle. They become important for hypothetical small black holes in their final stages or perhaps even for primordial black holes in the early universe.

A problem of principle
The non-stationarity of black holes constitutes a problem that may be safely ignored in most astrophysical contexts but it becomes important, or even crucial, in a fundamental physics context. Let us recall a few motivations for the study of dynamical black holes (a more extensive discussion can be found in [122,45]). The first problem to note is the fact that black holes are embedded in a universe, that is, they are not asymptotically flat. Asymptotic flatness is an essential ingredient in the proof of the no-hair theorems of scalar-tensor gravity, which provide a fundamental understanding of black holes in these theories. While there is literature on asymptotically de Sitter or anti-de Sitter black holes, it is much more realistic to consider asymptotically FLRW black holes. In the astrophysical context, the cosmological asymptotics cannot be neglected for primordial black holes.
Theories of gravity alternative to GR, aiming at explaining the cosmic acceleration without dark energy (e.g., f (R) gravity [24,26]) contain a built-in, time-dependent cosmological "constant". Therefore, black holes in theories of this class which are relevant for present-day cosmology are not asymptotically flat. Understanding black holes is an important part of understanding a theory of gravity, and the FLRW asymptotics are essential here.
Dynamical inhomogeneous spherically symmetric universes are theoretical testbeds which allow one to probe ideas about the backreaction of inhomogeneities, or the idea that we live in a giant void, which have been proposed as alternatives to dark energy to explain the cosmic acceleration (see the review [19]).
Black hole mechanics and thermodynamics were developed for stationary black holes with event horizons, which are null surfaces. Realistic black holes are instead dynamical and have apparent horizons, which are timelike or spacelike surfaces and they may change their causal nature during the evolution.
If primordial black holes were formed in the early universe, at some time they could have had a size non-negligible in comparison with the Hubble scale ∼ H −1 and very dynamical horizons. A relevant question is how fast they could have accreted material and have grown, even though primordial black holes are unlikely to be a dominant component of dark matter today. Another interesting problem of principle is the accretion of dark energy, especially phantom energy, onto black holes [9,28,73,39,87,57,62,125,126,59,68,8,105,27]. The simplest model is, of course, spherical accretion but this is already a very non-trivial problem requiring the introduction of some toy model.
Another motivation of principle is the study of the spatial variation of the fundamental constants of physics. This idea goes back to Dirac and is partially implemented in scalartensor gravity, in which the gravitational coupling strength G ef f becomes a dynamical field (roughly speaking, the inverse of the Brans-Dicke scalar field φ). While the time variation of G ef f has been studied extensively in scalar-tensor cosmology, its spatial variation requires the analysis of inhomogeneous universes, which have been the subject of a much smaller literature.

A practical problem
There is a much more practical and urgent problem in gravitational physics: highly dynamical black holes are very important in certain astrophysical applications. The first direct detection of gravitational waves by the LIGO interferometers involved a black hole merger [1]. In the final approach and during coalescence of these black holes, the spacetime is extremely dynamical. The signal to noise ratio in LIGO interferometers is typically low and one needs to match the signal to templates calculated theoretically. Numerical simulations of black hole mergers are used to generate banks of templates of gravitational waveforms for LIGO detection. These simulations do not rely on event horizons, but they use instead apparent and trapping horizons. Event horizons are essentially useless for practical purposes such as the generation of these banks of templates. Even though an unusually loud gravitational wave signal from a binary black hole merger can be spotted without templates, the determination of the orbital parameters of the black hole binary requires fitting the signal with templates computed by using apparent horizons. In the calculations producing these templates, "black holes" are identified with outermost marginally trapped surfaces and apparent horizons (e.g., [129,14,29,33]). It is fair to say that, in practice, in astrophysics we use apparent horizons and not event horizons.
It is impossible to discuss dynamical black holes without speaking of apparent horizons. Therefore, in the next section we introduce apparent and trapping horizons and their problems. Sec. 3 presents a selection of exact solutions of various theories of gravity describing black holes embedded in cosmological "backgrounds". 1 The last section outlines open problems.

Apparent horizons and their problems
In the words of Rindler [113], a horizon is a frontier between things observable and things unobservable. The horizon, a product of strong gravity, characterizes a black hole.
In black hole thermodynamics, if the "background" is not Minkowskian, the internal energy appearing in the first law of thermodynamics must be defined carefully. This is identified with a quasi-local energy. In spherical symmetry, this is the Misner-Sharp-Hernandez energy [95,70], to which the more general Hawking-Hayward quasilocal energy [64,65] reduces. In spherical symmetry, the Hawking-Hayward/Misner-Sharp-Hernandez quasilocal energy is in practice closely related to the notion of apparent horizon.
Let us begin by reviewing the concepts of null geodesic congruence and trapped surface necessary to define apparent horizons [133,112]. Consider a congruence of null geodesics, each with tangent l a = dx a /dλ, where λ is an affine parameter along the geodesic. (The assumption of an affine parameter is not mandatory, but it simplifies several of the following equations.) With affine parametrization, the tangent satisfies l a l a = 0 and l c ∇ c l a = 0. The metric h ab in the 2-space orthogonal to l a is determined as follows. Pick another null vector field s a such that s c s c = 0 and l c s c = −1, then (2.1) The tensor h ab is purely spatial and h a b is the projection operator on the 2-space orthogonal to l a . The choice of s a is not unique, but the geometric quantities of interest to us do not depend on it once l a is fixed [112]. Let η a be the geodesic deviation vector and define the tensor orthogonal to the null geodesics. The transverse part of the deviation vector is and the orthogonal component of l c ∇ c η a , denoted by a tilde, is Now decompose the transverse tensor B ab as where the expansion θ = ∇ c l c propagates according to the Raychaudhuri equation [133,112] If the null geodesic congruence is not affinely parametrized, the geodesic equation is (where κ is sometimes used as a possible notion of surface gravity) and The difference between Eq. (2.7) and our previous definition θ ≡ ∇ c l c is due to the fact that the geodesics in Eq. (2.7) are not required to be affinely parametrized, while the definition θ ≡ ∇ c l c does require an affine parameter and κ measures precisely the failure of geodesics to be affinely parameterized. Since κ arises purely from the choice of parametrization, its significance as a surface gravity is questionable and, indeed, several alternative definitions of "surface gravity" exist in the literature (see Refs. [99,111]) for reviews).
With non-affine parametrization, the Raychaudhuri equation then becomes (2.10) A compact and orientable surface has two independent directions orthogonal to it, corresponding to ingoing and outgoing null geodesics with tangents l a and n a , respectively. Let us give some basic definitions for closed 2-surfaces.
Definition: A normal surface corresponds to θ l > 0 and θ n < 0.
Definition: A trapped surface corresponds to θ l < 0 and θ n < 0.
The outgoing, in addition to the ingoing, future-directed null rays converge here instead of diverging and outward-propagating light is dragged back by strong gravity.
Definition: A marginally outer trapped (or marginal) surface (MOTS) corresponds to θ l = 0 (where l a is the outgoing null normal to the surface) and θ n < 0.
Definition: An untrapped surface is one with θ l θ n < 0.
Definition: A marginally outer trapped tube (MOTT) is a 3-dimensional surface which can be foliated entirely by marginally outer trapped (2-dimensional) surfaces.
A 1965 result by Penrose [110] states that if a spacetime contains a trapped surface, the null energy condition holds, and there is a non-compact Cauchy surface for the spacetime, then this spacetime contains a singularity. Timelike naked singularities are believed to be unphysical since they spoil the initial value problem and cause the loss of predictability. Therefore, black hole horizons assume a crucial role in physics, which is hiding these singularities from the world outside the horizon.
Trapped surfaces are essential features in the concept of black hole. "Horizons" of practical utility are identified with boundaries of spacetime regions containing trapped surfaces. 2

Event horizons
Definition: An event horizon is a connected component of the boundary ∂ (J − (I + )) of the causal past J − (I + ) of future null infinity I + .
An event horizon is a causal boundary separating a region from which nothing can come out to reach a distant observer from a region in which signals can be sent out and eventually arrive to this observer. An event horizon is generated by the null geodesics which fail to reach infinity. Provided that it is smooth, an event horizon is a null hypersurface.
In order to define and locate an event horizon, one must know all the future history of spacetime: the event horizon is a globally defined concept and has teleological nature. As an example, consider a collapsing Vaidya spacetime described by the geometry [131] where v is a retarded time, r is the areal radius, dΩ 2 (2) ≡ dθ 2 + sin 2 θ dϕ 2 is the line element on the unit 2-sphere, and m(v) is a time-and radius-dependent mass function. The apparent horizon is located by the root of the equation ∇ c r∇ c r = g rr = 0, which yields r AH = 2m(v). The apparent horizon is distinct from the event horizon, which is given approximately by the expression [10,20,15,97,137] whereṁ ≡ dm/dv. The phenomenology of the event horizons is reported in several studies (e.g., [10,20,15,97,137] and we will not repeat the analysis here. To summarize the situation, an event horizon forms and grows starting from the centre and an observer can cross it and be unaware of it even though his or her causal past consists entirely of a portion of Minkowski space. This event horizon cannot be detected by this observer with a physical experiment and it "knows" about events belonging to a spacetime region very far away and in its future, but not causally connected to it (a property called "clarvoyance") [7,15,16]. Properly speaking, an event horizon H is a tube in spacetime. There is some language abuse in the literature, which consists of referring to the intersections of H with surfaces of constant time (which produce 2-surfaces) as "event horizons".

Apparent and trapping horizons
Definition: A future apparent horizon is the closure of a 3-surface which is foliated by marginal surfaces.
An apparent horizon is defined by the conditions on the time slicings [66] θ l = 0 , (2.13) where θ l and θ n are the expansions of the future-directed outgoing and ingoing null geodesic congruences, respectively, which have tangent 3 fields l a and n a (outgoing null rays momentarily stop expanding and turn around at the horizon). In simple words, an event horizon is a surface from the interior of which nothing will ever emerge, while an apparent horizon is a surface from the interior of which nothing can escape now. The inequality (2.14) distinguishes between black holes and white holes. Apparent horizons are defined quasi-locally and they do not require the knowledge of the global structure of spacetime, including future null infinity. In this sense, they are much more practical entities than event horizons. However, apparent horizons have a serious limitation: they depend on the choice of the foliation of the 3-surface with marginal surfaces. This foliation-dependence is exemplified dramatically by the fact that non-symmetric slicings of the Schwarzschild spacetime exist for which there is no apparent horizon [134,120]. In nonstationary situations, apparent horizons and event horizons do not coincide. This foliationdependence problem is overlooked: it implies that the very existence of a dynamical black hole depends on the observer. This problem is alleviated (or, for purely practical purposes, largely eliminated), but not solved, in spherical symmetry. In this case, the apparent horizons coincide in all spherical foliations [51].
In GR, a black hole apparent horizon lies inside the event horizon provided that the null curvature condition R ab l a l b ≥ 0 for all null vectors l a is satisfied. However, Hawking radiation itself violates the weak and the null energy conditions, as do quantum matter and non-minimally coupled scalars, hence apparent horizons cannot always be expected to lie inside event horizons in the presence of quantum fields or in theories of gravity alternative to GR.
Definition: A future outer trapping horizon (FOTH) is the closure of a surface (usually a 3-surface) foliated by marginal surfaces such that on its 2-dimensional "time slicings" [66] θ l = 0 , The last condition distinguishes between inner and outer horizons and between apparent and trapping horizons (the sign distinguishes between future and past horizons).
Definition: A past inner trapping horizon (PITH) is defined by exchanging l a with n a and reversing the signs of the inequalities, The PITH identifies a white hole or a cosmological horizon. As one moves just inside an outer trapping horizon, one encounters trapped surfaces, while trapped surfaces are encountered as as one moves just outside an inner trapping horizon. As an example, consider the Reissner-Nordström black hole with the natural spherically symmetric foliation. The event horizon r = r + is a future outer trapping horizon (FOTH), the inner (Cauchy) horizon r = r − is a future inner trapping horizon (FITH), while the white hole horizons are past trapping horizons (PTHs).
Black hole trapping horizons have been associated with thermodynamics. It is claimed that it is the trapping horizon area, and not the area of the event horizon, which should be associated with entropy in black hole thermodynamics [63,71,32,96]. This claim is somewhat controversial [121,34,98]. The Parikh-Wilczek [108] "tunneling" approach, originally applied to study Hawking radiation from stationary black holes, is in principle applicable also to apparent and trapping horizons.
In general, trapping horizons do not coincide with event horizons. Dramatic examples are spacetimes which possess trapping horizons but not event horizons [115,67].

Spherical symmetry
As usual, assuming a symmetry greatly simplifies a physical problem, and this is even more so for the assumption of spherical symmetry when discussing apparent horizons and dynamical black holes. The notion of internal energy in black hole thermodynamics is crucial for the first law and, in spherical symmetry, it is universally identified with the Misner-Sharp-Hernandez mass notion [95,70].
The Misner-Sharp-Hernandez mass [95,70] is defined in GR 4 and in spherical symmetry. The more general Hawking-Hayward quasilocal energy [64,65] reduces to the Misner-Sharp-Hernandez notion in spherical symmetry. A spherically symmetric line element can be written as where R is the areal radius, a geometrical quantity defined in a coordinate-invariant way using the area A of 2-spheres of symmetry, R = A/(4π), and h ab is the 2-metric on the submanifold orthogonal to the time and radial directions. Then the Misner-Sharp-Hernandez mass M M SH is defined in an invariant way by the scalar equation [95,70] It is often useful to use the formalism of Nielsen and Visser [100], in which a general spherical line element is written as where M (t, R), a posteriori, turns out to be the Misner-Sharp-Hernandez mass. This notation interpolates between a gauge originally used to describe wormholes and the more familiar Schwarzschild-like coordinates. The line element (2.23) is recast in Painlevé-Gullstrand coordinates as where the functions φ(τ, R) and M (τ, R) are defined implicitly. By introducing the quantities In this gauge, the outgoing radial null geodesics have tangents while the ingoing radial null geodesics have tangents 30) and the expansions of the radial null geodesic congruences are The apparent horizons are located by The inverse metric is In practice, the condition g RR = 0 is a very convenient recipe to locate the apparent horizons in spherical symmetry.
the horizon is outer. The condition for the apparent horizon to be also a trapping horizon iṡ If matter satisfies the null energy condition, and assuming the Einstein equations, the area of the apparent horizon cannot decrease. The Nielsen-Visser surface gravity is

A selection of exact solutions in various theories of gravity
Having introduced the basic quantities, we can now proceed to see a little bestiary of analytical solutions of the field equations of various theories of gravity describing dynamical black holes embedded in cosmological spacetimes at least part of the time. Not many such exact solutions are known and some of them have physical problems such as negative energy densities, usually in spatial regions close to a black hole apparent horizon. Nevertheless, a complete review would be too long 5 and a somewhat arbitrary selection needs to be made, but the geometries appearing most often in the literature are reported below. For extra discussion of cosmological black holes in alternative gravity see [122,69,45,130].

Schwarzschild-de Sitter/Kottler spacetime
The oldest known cosmological black hole geometry is the Schwarzschild-de Sitter/Kottler solution of GR and of many other theories [77]. It solves the vacuum Einstein equations with positive cosmological constant Λ where R ab is the Ricci tensor and R ≡ g ab R ab is its trace. The Jebsen-Birkhoff theorem familiar from GR textbooks and stating the uniqueness of the Schwarzschild solution under the assumptions of vacuum, spherical symmetry, and asymptotic flatness can easily be generalized to show that the Schwarzschild-de Sitter/Kottler geometry is the unique spherically symmetric solution of Eq. (3.37) with Λ > 0 [127,119,50]. The Schwarzschild-de Sitter/Kottler line element is in static coordinates, where H = Λ/3. This geometry is locally static in the region comprised between the black hole and the cosmological horizons R 1 < R < R 2 (with R 1,2 defined below), where a timelike Killing vector exists. The apparent horizons are located by the usual equation The formal roots of this cubic equation are with sin(3ψ) = 3 √ 3 mH. Here m and H are positive, which implies that R 3 < 0 and there are at most two apparent horizons. When R 1 and R 2 are real, R 1 is a black hole apparent horizon while R 2 is a cosmological apparent horizon. Both apparent horizons are null surfaces, due to the (locally) static character of this geometry. Three situations are possible: • Two apparent horizons exist if 0 < sin(3ψ) < 1.
• For sin(3ψ) > 1 there is a naked singularity. The physical interpretation of this situation is that the black hole horizon becomes larger than the cosmological one and effectively disappears so, in the region below the cosmological horizon, the central singularity is not screened by a black hole horizon.

McVittie solution
The The original motivation for the McVittie solution was the study of the cosmological expansion (which in those days was still a fairly recent discovery) on local systems. This problem was later investigated also by Einstein, who was unaware of McVittie's work and led to the Einstein-Straus vacuole, or Swiss-cheese model [40,41]. McVittie made a crucial simplifying assumption in the derivation of his solution using spherical coordinates [93]: the "no-accretion condition" G 1 0 = 0. This implies T 1 0 = 0 due to the Einstein equations and, therefore, zero radial energy flow onto, or out of, the central inhomogeneity. The McVittie line element in isotropic coordinates is [93]  which leads to The cross-term in dtdR can be eliminated defining the new time coordinate T (t, R) by where F (t, R) is an integrating factor. The cross-term in the new coordinate system vanishes by imposing The line element is then recast as wherer = m/2 corresponds to R = 2m a(t) = 2m 0 , that is, the finite radius singularity does not expand with the rest of the universe in which it is embedded. The McVittie metric admits arbitrary FLRW "backgrounds" generated by fluids with any constant equation of state. Let us restrict to a fluid which reduces to dust at spatial infinity (w = 0), for the sake of simplicity. The fluid pressure is The apparent horizons are located at the roots of the equation This equation is the same cubic already seen in the Schwarzschild-de Sitter/Kottler case, but now it has a time-dependent coefficient H(t). Its roots R 1,2 (t) are given again by the expression already seen, but now with time-dependent coefficient H(t). This means that the location of the apparent horizons depends on time (i.e., on the comoving time t of the FLRW "background"). Both apparent horizons exist if mH(t) < 1/(3 √ 3), an inequality which is satisfied only if t > t * . The critical time at which mH(t) = 1/(3 √ 3) for a dust "background" is There are three distinct epochs in the history of this inhomogeneous universe ( Fig. 1): 3 H(t) and both R 1 (t) and R 2 (t) are complex. There are no apparent horizons.  3 H(t) . A black hole apparent horizon cannot be accommodated in this small universe and, at t < t * there is a naked singularity at R = 2m 0 . At t * an instantaneous black hole apparent horizon and a cosmological apparent horizon appear together at areal radius R 1 (t * ) = R 2 (t * ) = 1 √ 3 H(t * ) , in analogy with the extremal Nariai black hole. As t > t * , this single horizon splits into an evolving black hole apparent horizon surrounded by an evolving cosmological horizon. The black hole apparent horizon shrinks, asymptoting to the 2m 0 singularity as t → +∞.
Motivated by the discovery of dark energy, and by recurrent claims by astronomers that this dark energy may have a phantom equation of state with w ≡ P/ρ < −1 (although such a form of dark energy seems to be in extreme conflict with known theory), one can consider a phantom "background" universe in the McVittie solution. In this case the FLRW cosmic scale factor is where t rip is the time of the Big Rip and A is a positive constant. In this situation, the behaviour of the apparent horizons (Fig. 2) appears to be the time reversal of the behaviour already seen in a non-phantom universe.  An idealized interior solution for the McVittie exterior metric has been found [101]. It describes a relativistic star of uniform density in a FLRW "background" [101]. The equivalent of the Tolman-Oppenheimer-Volkov equation in the interior of a McVittie star is derived in [47].
Recent works on the McVittie spacetime have studied its conformal structure [76,80,81,37], which entails integrating numerically the null geodesics or deriving general analytical results upon making some assumptions about the cosmic expansion. Lake & Abdelqader [80] have found that null geodesics asymptote to the singularity without entering it. Depending on the form of the scale factor, a bifurcation surface may appear which splits the spacetime boundary into a black hole horizon in the future and a white hole horizon in the past (it is not clear whether this is due directly to the McVittie no-accretion condition). da Silva, Fontanini, and Guariento [37] find that the presence of this white hole horizon depends crucially on the expansion history of the universe.
It was found recently that the McVittie geometry is also a solution of cuscuton theory, a special Hoȓava-Lifschitz theory, and of shape dynamics [4,58,61].

Generalized McVittie solutions
It is interesting to remove the no-accretion restriction from the McVittie solutions [47]. In principle, this becomes the "Synge approach" in which one postulates a metric and runs the Einstein equations R ab − Rg ab /2 = 8πT ab from the left to the right to derive the form of an effective stress-energy tensor T ab that generates this metric. Usually this approach is fruitless because the effective T ab thus generated is physically pathological and violates all the energy conditions. However, reasonable matter sources exist in this case. Assume the line element to be where m(t) ≥ 0 , (3.59) and where now the new function m(t) is not determined by the McVittie no-accretion condition. The mixed Einstein tensor is For the special subclass of solutions which has m = m 0 =const., the quantity this is the non-rotating Thakurta solution [128], which will re-appear later in our discussion. Atr = m/2, the quantity C reduces to We must now discuss the matter source. If this is a single perfect fluid, only the McVittie solutions are possible. However, imperfect fluids can be matter sources for generalized McVittie spacetimes. Two cases have been studied [47] and are reported below. In addition, since the McVittie condition has been removed and we now have a radial energy flow, we must model it somehow and an imperfect fluid is the simplest model for this purpose.

Imperfect fluid and no radial mass flow
In this case the right hand side of the Einstein equations (3.43) contains the matter stressenergy tensor T ab = (P + ρ) u a u b + P g ab + q a u b + q b u a , (3.70) where the purely spatial vector q c describes a radial energy flow, and u c u c = −1. The Einstein equations yielḋ The radial energy flow, the area A of a 2-sphere, and the accretion rate are related by [47] M (t) = −aB 2 Aq . (3.73) In the case of inflow (q < 0), on a sphere of radiusr ≫ m this accretion rate becomesṀ ≃ aA |q|. The mass M changes due to inflow of matter and to the evolution of the cosmological fluid in it. The energy density and the pressure of the fluid are The evolution of the quantity C is regulated by the generalized Raychaudhuri equation [47] (3.76)

Imperfect fluid and radial mass flow
A more general situation is the one in which both an imperfect fluid and a radial flow of material and energy are allowed, as described by the stress-energy tensor (3.70) with (i.e., with non-vanishing radial velocity) and with spacelike flux density In this case the accretion rate onto the central object iṡ 79) and the energy density is found to be [47] 8πρ = (3.80) The generalized McVittie geometry is also a solution of Horndeski theory [5].

The non-rotating Thakurta solution
The choice M (t) = m 0 a(t) selects a special subclass of generalized McVittie solutions of the Einstein equations, the non-rotating Thakurta solution [128]. The "comoving mass" attractor family coincides with the non-rotating Thakurta solution [128] of GR discovered in 1981. This is the special case, corresponding to zero angular momentum, of the Thakurta solution describing a Kerr black hole embedded in a FLRW background [128]. Recently, the nonrotating Thakurta solution has been discussed, with different purposes, in Refs. [36,94]. The apparent horizons exhibit the same phenomenology of apparent horizons described by a C- The non-rotating Thakurta subclass is not a mere curiosity, but it is important because it is a late-time attractor of generalized McVittie solutions [57]. Given that the Jebsen-Birkhoff theorem [75,18] fails in the presence of matter and/or when the asymptotics are not Minkowskian, and then no "general solutions" are known, this result is significant because it provides a general solution in the generalized McVittie class, provided that one waits long enough (either in a universe expanding forever or in a phantom-dominated universe approaching a Big Rip singularity at a finite future). The non-rotating Thakurta solution is generic under certain assumptions, in the sense that all other generalized McVittie solutions approach it at late times.
Another peculiarity of this "comoving mass" Thakurta attractor resides in the structure of its apparent horizons. In general, in spherical symmetry, the radii of the apparent horizons can only be found numerically because Eq. (2.22) locating them is trascendental. It is very rare to be able to obtain analytical explicit expressions for apparent horizon radii, which would be useful, for example, for quantum field theory calculations in these curved spacetimes to determine Hawking radiation fluxes and so on. This task is possible for the non-rotating Thakurta family of solutions, for which the areal radii of the cosmological and black hole apparent horizons are [57]

The Husain-Martinez-Nuñez solution
The 1994 Husain-Martinez-Nuñez solution of the Einstein equations brings in new phenomenology of the apparent horizons with respect to what we have seen thus far. The Husain-Martinez-Nuñez spacetime describes an inhomogeneous universe with a spatially flat FLRW "background" sourced by a free, minimally coupled, scalar field. The Einstein equations (3.43) assume the form while the free scalar field φ obeys the curved space Klein-Gordon equation The Husain-Martinez-Nuñez line element is where A 0 , B 0 , C, D ≥ 0 are constants, α = ± √ 3/2, and η > 0. The additive constant B 0 becomes irrelevant and can be dropped if A 0 = 0. When A 0 = 0, the Husain-Martinez-Nuñez line element degenerates into the static Fisher one [55] where V (r) = 1 − 2µ/r, µ and ν are parameters, and the Fisher scalar field is The Fisher solution also goes by the names Buchdahl-Janis-Newman-Winicour-Wyman solution because it was rediscovered many times by these authors [17,74,22,136]. It always hosts a naked singularity at r = 2C and is asymptotically flat. The general Husain-Martinez-Nuñez metric is conformal to the Fisher metric with conformal factor Ω(η) = √ A 0 η + B 0 equal to the scale factor a(η) of the "background" FLRW space and with only two possible values of the parameter ν. In the following we set B 0 = 0.
The Husain-Martinez-Nuñez geometry is asymptotically FLRW as r → +∞ and is exactly FLRW if C = 0, in which case the constant A 0 can be eliminated by rescaling the time η. The Ricci scalar is This Ricci scalar shows that there is a spacetime singularity at r = 2C (for both values of the parameter α). The matter scalar field φ also diverges there, and there is a Big Bang singularity at η = 0. We have 2C < r < +∞ and the value r = 2C corresponds to zero areal radius We can transform the time coordinate and use the comoving time t of the "background" FLRW space instead of the conformal time η related to t by dt = adη = √ A 0 ηdη [45]. Then we have and The Husain-Martinez-Nuñez solution in comoving time reads The areal radius increases with r for r > 2C. In terms of the areal radius R, and using the notation we have R(t, r) = a(t)rA A time-radius cross-term in the line element is eliminated by introducing a new time coordinate The choice The apparent horizons are located by the roots of the equation g RR = 0, which reads (3.100) As r → +∞, corresponding to R → +∞, Eq. (3.100) reduces to R ≃ H −1 , the radius of the cosmological apparent horizon in FLRW. Let x ≡ C/r, then the equation locating the apparent horizons is The left hand side of this equation is while the right hand side is therefore the apparent horizons radii are expressed in parametric form by with parameter x. The numerical solution for the radii of these apparent horizons is shown in Fig. 3.
If α = √ 3/2, between the Big Bang and a critical time t * there is only one expanding apparent horizon, then two other apparent horizons are created at the critical comoving time t * . One of them is a cosmological apparent horizon which expands forever, while the other is a black hole apparent horizon which contracts until it meets the first (expanding) black hole apparent horizon [72]. When they meet, these two apparent horizons "annihilate" and a naked singularity appears at R = 0 in a FLRW universe [72,45]. This phenomenology of apparent horizons is dubbed "S-curve" behaviour from the shape of the curve representing the evolution of the apparent horizon radii in Fig. 3. This "S-curve" phenomenology appears also in Lemaître-Tolman-Bondi spacetimes, in which a dust fluid curves the cosmological "background" [21]. Multiple "S"s are also possible (e.g., five of them may appear in these Lemaître-Tolman-Bondi solutions). The scalar field is regular on the apparent horizons.
For the parameter value α = − √ 3/2, there is only one cosmological apparent horizon and the universe contains a naked singularity at R = 0.
The apparent horizons are spacelike: the normal vector to these surfaces always lies inside the light cone in an (η, r) diagram [72], in agreement with a general result of Booth, Brits, and Gonzalez [21] stating that a trapping horizon created by a massless scalar field must be spacelike.
The singularity at R = 0 is timelike for both values of the parameter α. The central black hole in this spacetime is created with the universe in the Big Bang and not in a collapse process (although sometimes the literature states otherwise).
Clifton's 2006 solution of f (R) = R n gravity [30] (where R is the Ricci scalar) and some perfect fluid solutions of Brans-Dicke gravity found by Clifton, Mota, and Barrow [31] exhibit the same S-curve phenomenology of apparent horizons originally discovered in the Husain-Martinez-Nuñez spacetime [44,53]. f (R) gravity is described by the fourth order field equations where f ′ (R) ≡ df /dR. It can be shown that f (R) theories are equivalent to a subclass of Brans-Dicke theories with effective Brans-Dicke field φ = f ′ (R), coupling parameter ω = 0, and the special (implicit) potential (see Refs. [123,38,106] for reviews).

The Fonarev and phantom Fonarev solutions
Another inhomogeneous universe representing a central object embedded in a FLRW "background" is the Fonarev solution of GR [56], which generalizes the Husain-Martinez-Nuñez spacetime (the latter has a free scalar field as the matter source) to the case in which the scalar field acquires the exponential potential where V 0 and λ are positive constants. The coupled Einstein-Klein-Gordon equations are now The Fonarev line element and scalar field are [56] 115) and where w and A 0 are constants, while η is the conformal time of the FLRW substrate. Choosing A 0 = 1 and restricting to the parameter value w = 0, the Fonarev solution reduces to spatially flat FLRW. If a ≡ 1 and α = 1, it reduces to the Schwarzschild geometry, while if λ = ± √ 6 and V 0 = 0, it reproduces the Husain-Martinez-Nuñez spacetime. A generalized Fonarev solution corresponding to a phantom scalar field solution of the Einstein equations was presented in [57] and is obtained from the Fonarev solution by means of the transformation The equation locating the apparent horizons is a quartic which has only two real positive roots, corresponding to a cosmological apparent horizon of radius R c and to a black hole apparent horizon of radius R b with the same qualitative behaviour of the apparent horizons of the (generalized) McVittie solution.
One can regard the Fonarev solution of GR as the Einstein frame version of a solution of Brans-Dicke gravity, which can then be mapped back to the Jordan frame 6 [46]. The conformal cousin of the Fonarev geometry is a 4-parameter family of solutions which are spherical, non-static, and asymptotically FLRW. The solutions of this family which have been studied explicitly contain only wormholes and naked singularities. Special cases of this 4parameter family provide solutions of vacuum f (R) = R n gravity [46].

Other GR solutions
Other better-known solutions of GR which may be interpreted as black holes embedded in cosmological substrata include Swiss-cheese models [40,41], Lemaître-Tolman-Bondi black holes (with a dust-dominated FLRW "background"), members of the large Barnes family of solutions [13], and the Sultana-Dyer solution [124]. Several other analytical solutions of GR exist which describe central inhomogeneities in FLRW "backgrounds" [56,132,109,114,23,11,35,83,84]. Many do not have reasonable matter sources and the energy density becomes negative in certain spacetime regions. Often these solutions are obtained by conformally transforming a stationary black hole solution or by performing a Kerr-Schild transformation of a stationary black hole metric g ab , g ab →ḡ ab = g ab + λk a k b , (3.117) where λ is a function of the spacetime position and the vector k a is null and geodesic with respect to both g ab andḡ ab [78,79,89,90,91,92].

Perfect fluid solutions of Brans-Dicke gravity
Perfect fluid solutions of Brans-Dicke gravity describing dynamical cosmological black holes in certain regions of the parameter space were found by Clifton, Mota, and Barrow [31]. The Brans-Dicke field equations are where ω = −3/2 is the Brans-Dicke coupling parameter and T is the trace of the matter stress-energy tensor T ab = (P + ρ) u a u b + P g ab (3.120) describing a perfect fluid with constant equation of state while the Brans-Dicke field φ is free and massless. The line element of this family of solutions is spherical, inhomogeneous, and asymptotically FLRW: where while the Brans-Dicke scalar field and the fluid energy density are The areal radii of the apparent horizons of this class of solutions were studied in [53], disclosing a rich variety of behaviours as the three parameters vary. Some possibilities, corresponding to selected regions of parameter space, are shown in Figs. 4-7 (we refer the reader to [53] for details).

Is there a relation between S-curve and C-curve?
Thus far, we have encountered two main phenomenological classes describing the time behaviour of apparent horizons in spherical inhomogeneus spacetimes: the "C-curve" behaviour exemplified by McVittie solutions and the "S-curve" phenomenology originally discovered in the Husain-Martinez-Nuñez solution of GR. These two phenomenologies are encountered in Figure 4: The areal radii of the apparent horizons versus the FLRW comoving time (both in units of (ma 0 ) 1/(1−β) ) for ω = −17/12. The dashed curve corresponds to dust (γ = 1) and the solid curve corresponds to both radiation (γ = 4/3) and stiff matter (γ = 2). For dust, a single apparent horizon expands to a maximum size and then shrinks.  there is a single horizon for all three values of γ. As time goes by, two more apparent horizons appear. Two of these horizons eventually merge and disappear, leaving a naked singularity in a FLRW universe, which has its own cosmological horizon. Fig. 7 zooms in on the third curve, which here appears flattened along the time axis.  GR and in other theories of gravity. We know that other behaviours are in principle possible (cf. Sec. 3.6) and that the catalog of analytic solutions with these physical properties is scarce in GR and even leaner in alternative theories of gravity. However, the tentative assumption that these two behaviours have some generality has some basis for the moment. One is then led to wonder whether the "S-curve" and the "C-curve" behaviours are completely disconnected and mutually exclusive, or whether there may be some relation between them. As the Brans-Dicke coupling ω tends to infinity, the Clifton-Mota-Barrow solution of Brans-Dicke gravity discussed in the previous subsection asymptotes to the comoving mass non-rotating Thakurta solution of GR [45]. While the ω → ∞ limit of Brans-Dicke theory is expected to reproduce a corresponding GR solution [135], exceptions are known [88,116,117,118,107,6] and often explained [12,42,43]. In this case the non-rotating Thakurta solution is indeed the GR limit of the Clifton-Mota-Barrow solution, which was not realized nor expected a priori. But the most interesting thing is the behaviour of the apparent horizons of the Clifton-Mota-Barrow geometry in this limit: the S-curve describing the apparent horizons areal radii in the (t, R) plane reduces to a C-curve as ω → ∞ because the lower bend of the S-curve is pushed to infinity [48], as shown in Fig. 8. Is the C-curve described by the apparent horizons radii of a spherical, asymptotically FLRW spacetime always a limit of an S-curve? Does the relation between these two phenomenologies have any physical meaning or are we just looking at a coincidence? These questions are currently under investigation.

Conclusions
To conclude this brief excursion on the subject of dynamical cosmological black holes and their relations with modified gravity, one should first realize that, while it is auspicable to find new inhomogeneous solutions amenable to this physical interpretation, it is usually harder to find whether apparent horizons exist, and to determine their location, nature, and dynamical behaviour, and some effort should be devoted to these goals.
There are many open questions, some of which are more general than the subject of dynamical cosmological black holes, and are related to the nature of horizons, which ultimately define the very concept of black hole. The main question seem to be whether apparent horizons are the "right" surfaces to characterize black holes. It is a fact that apparent and trapping horizons are used in the numerical prediction of the waveforms of gravitational waves emitted by binary systems and used in their interferometric detection. The recent successes of the LIGO interferometers [1,2,3] in detecting the elusive gravitational waves from black hole mergers seem to answer the question affirmatively, but it is in principle possible that a different notion of horizon could be as successful as that of apparent/trapping horizon. In the meantime, the gravitational wave community is not apologetic in using apparent and trapping horizons and in disregarding the teleological event horizon.
A different but related question, which has not been addressed here, is whether apparent horizons are the "right" surfaces to use in the study of the thermodynamics of dynamical black holes. There are strong claims that they are, but no convincing definitive proof has been provided. It is difficult to perform calculations of quantum field theory in curved space even on a prescribed time-dependent black hole background (that is, neglecting backreaction). The situation at the moment is summarized in the fact that the tunneling method provides a definite and general result for the Hawking temperature while other methods, at best, have produced results only for specific and very particular background metrics and a general computation cannot be completed as is done in the tunneling method. However, at present there is no guarantee that the result and procedure provided by the tunneling method are actually correct. Progress in this direction is very slow. At the very least, apparent horizons would require some adiabatic approximation if they are to make sense for dynamical black holes, otherwise one is probably discussing non-equilibrium thermodynamics, which is a tall order in the context of a problem which is proving to be already very difficult in equilibrium thermodynamics.
Assuming that apparent horizons are the "correct" notion of horizon, their biggest problem is no doubt their foliation-dependence, which amounts to the fact that the very existence of a black hole seems to depend on the observer. This problem is exemplified dramatically by the realization that there exist slicings of the Schwarzschild spacetime with no apparent horizons [134,120]. The fact that these slicings are rather contrived is not of much consolation. From the practical point of view, the problem is mitigated by the fact that, in spherical symmetry, the apparent horizons coincide in all spherically symmetric foliations [51], but this fact does not really solve the problem and, moreover, no analogous result has been proved beyond spherical symmetry.
In our examples of apparent horizon phenomenology, we have restricted ourselves to spherical symmetry. Evolving apparent horizons exhibit a rather rich phenomenology and dynamics, but there seem to be two main phenomenological classes for the behaviour of the areal radii R AH (t) of apparent horizons versus the comoving time of the FLRW substratum in which the spherical inhomogeneities are embedded. These are the "C-curve" first encountered in the McVittie spacetime of GR [93] and the "S-curve" first discovered in the Husain-Martinez-Nuñez solution of the same theory [72] and then rediscovered in f (R) and Brans-Dicke gravity [44,53]. There are hints that some relation may exist between these two classes, with the C-curve being some limit of the S-curve, but the situation is not clear at present.
A final question is provided by the long-standing issue of cosmic expansion versus local dynamics which motivated the introduction of the McVittie metric in 1933 and of the Einstein-Straus Swiss-cheese model a decade later. This problem of principle, which at some point was brushed off as being irrelevant, keeps being discussed in the literature (see [25] for a 2010 review). We have seen already in the few examples discussed here that sometimes black hole apparent horizons expand (in rare cases they are even comoving with the cosmic substratum), whereas some other times they resist the cosmic expansion or even contract. No general rule appears here and the apparent horizons of dynamical cosmological black holes provide no definite answer to the puzzle. All these theoretical questions of classical gravity are worth understanding and investigating in the future. The recent detection of gravitational waves and the future development of gravitational wave astronomy seem to reformulate the issue of apparent horizons in more practical terms and to make its understanding a more pressing issue than it was before.