Exact cosmological black hole solutions in Scalar Tensor Vector Gravity

We find an exact solution of Scalar-Tensor-Vector Gravity field equations that represents a black hole embedded in an expanding universe. This is the first solution of the kind found in the theory. We analyze the properties of the apparent horizons as well as the essential singularities of the metric, and compare it with the McVittie spacetime of General Relativity. Depending on the cosmological model adopted and the value of the free parameter $\alpha $ of the theory, the solution describes a cosmological black hole, an inhomogeneity in an expanding universe, or a naked singularity. We use the latter result to set further constraints on the free parameters of the theory.


Introduction
"Probably the most beautiful of all existing theories". These words by Landau and Lifschitz [1] reflect the pleasant aesthetic experience induced on many of us by General Relativity (GR). The theory not only excels in simplicity, symmetry, unification strength, and fundamentality [2], but also has an outstanding predictive and explanatory power. Though Einstein himself remarked GR charm [3], he was quite aware that it was not the ultimate theory of gravitation. He struggled the last decades of his life searching for suitable generalizations of the theory that could accommodate electrodynamics and also include quantum effects.
Besides the inherent deficiencies in the theory, such as the problem of spacetime singularities, GR models do not succeed in reproducing rotation curves of nearby galaxies, mass profiles of galaxies clusters, some gravitational lensing effects, and cosmological data. A possible solution to these problems consist in modifying the right hand side of Einstein equations: a term with a cosmological constant is added and the existence of dark matter is postulated. From an ontological point of view, this approach is quite costly since we are assuming the existence of entities of unknown nature whose properties have never been measured to date [4][5][6].
We can follow a different strategy to explain the astronomical data: modify the theory of gravitation. This is the case of Scalar-Tensor-Vector Gravity (STVG), also dubbed MOdified gravity (MOG) [7]. In STVG, the effects of gravity are not only represented by a metric tensor field but also by a scalar and a vector field. Specifically, the universal constant G along with the massμ of the vector field are the dynamical scalar fields of the theory. When gravity is weak, the equations of the theory reduce to a modified acceleration law characterized by: 1) an enhanced Newtonian constant G = G N (1 + α), and 2) at certain scales, a repulsive Yukawa force term that counteracts the augmented Newtonian acceleration law, in such a way that in the Solar System GR is recovered. The first of the features mentioned above allows to reproduce the rotation curves of many galaxies [8][9][10], the dynamics of galactic clusters [11][12][13]and cosmological observations [14,15], without dark matter 1 .
We can classify the known solutions of the field equations of STVG in two main groups. On the one hand, vacuum and non-vacuum solutions for a given distribution of matter where the spacetime metric is asymptotically flat. This is the case of the Schwarzschild and Kerr STVG black holes 2 found by Moffat [28], and neutron star models constructed by Lopez Armengol and Romero [29]. On the other hand, there are cosmological solutions such as the ones derived by Roshan [30] and Jamali and collaborators [31]. Until now, a third class of solutions remains unexplored in the theory: metrics that represent an inhomogeneity in an expanding universe.
In General Relativity, McVittie [32] was the first to obtain an exact solution of Einstein field equations that corresponds to a central inhomogeneity embedded in a Friedmann-Lemaître-Robertson-Walker (FLRW) background. The McVittie metric and its generalization have been widely studied through the years (see for instance the works by Faraoni and Jacques [33] and Carrera and Giulini [34]). The investigation of such solutions has transcended GR to encompass alternative theories of gravitation [35]. The results of the studies of inhomogeneous spacetimes have direct astrophysical implications: a cosmological force acting on large scales can modify the structure of galaxies and clusters of galaxies, and inhibit accretion processes. The effects of the cosmological expansion, thus, need to be taken into account when modeling the evolution of structure in the universe.
In this work we present exact solutions of STVG that represent an inhomogeneity in an expanding spacetime, and analyze the corresponding properties. We distinguish the metrics that represent cosmological black holes and we compare them with the corresponding solutions in GR.
The paper is organized as follows. We provide a brief introduction to STVG in Section 2. Next, we introduce a solution to the field equations of the theory that represents an inhomogeneity in an expanding universe. In Section 4, we analyze the properties of the metric: singularities and apparent horizons, and in Section 5 we offer a discussion of the results obtained. The last section of the paper is devoted to the conclusions.

STVG action and field equations
The action 3 in STVG theory is [7]: The recent detection of a neutron star merger in gravitational waves [16] (GW170817), and the subsequent observation of the electromagnetic counterpart GRB 170814A [17,18] has been used to show that a large class of alternative theories of gravitation, for instance those in which photons suffer an additional Shapiro time delay, must be discarded [19,20]. As demonstrated by Green and collaborators [21], STVG survives such stringent test: both gravitational and electromagnetic travel on null geodesics in the theory.
2 Different aspects of the STVG black hole solutions have been extensively studied in the literature: accretion disks around Schwarzschild and Kerr STVG black holes [22], shadows cast by near-extremal Kerr STVG black holes [23], black hole superradiance in STVG [24], quasinormal modes of Schwarzschild STVG black holes [25], the process of acceleration and collimation of relativistic jets in Kerr STVG black holes [26], dynamics of neutral and charged particles around a Schwarzschild STVG black hole immersed in a weak magnetic field [27], among others. 3 As suggested by Moffat and Rahvar [9] and Moffat and Toth [36], we dismiss the scalar field ω, and we treat it as a constant, ω = 1. where Here, g µν is the spacetime metric, R denotes the Ricci scalar, and ∇ µ is the covariant derivative; φ µ stands for a Proca-type massive vector field,μ is its mass, and The scalar fields G(x) andμ(x) vary in space and time, and V (G), and V (μ) are the corresponding potentials. We adopt the metric signature η µν = diag(−1, +1, +1, +1). The term S M in the action refers to possible matter sources. The full energy-momentum tensor for the gravitational sources is: (2.9) Following the notation introduced above, T M µν denotes the ordinary matter energy-momentum and T S µν the scalar contributions to the energy-momentum tensor; T φ µν stands for the energymomentum tensor 4 of the field φ µ : The equation of motion for a test particle in coordinates x µ is given by where τ represents the particle proper time, and q is the coupling constant with the vector field. Moffat [28] postulates that the gravitational source charge q of the vector field φ µ is proportional to the mass of the source particle, q = ± αG N m. (2.12) Here, G N denotes Newton's gravitational constant, and α is a free dimensionless parameter. The positive value for the root is chosen (q > 0) to maintain a repulsive, gravitational Yukawalike force when the mass parameterμ is non-zero. We see, then, that in STVG the nature of the gravitational field has been modified with respect to GR in two ways: there is an enhanced gravitational constant G = G N (1 + α), and a vector field φ µ that exerts a gravitational Lorentz-type force on any material object through Eq. (2.11).
3 Solution of STVG field equations

Derivation of the metric
In order to derive the spacetime metric that represents an inhomogeneity in an expanding universe, we make the following assumptions: • The energy-momentum tensor has two components T µν = T M µν + T φ µν , where T M µν stands for the energy-momentum of the cosmological fluid, and T φ µν is the energy-momentum tensor for the vector field φ µ : Here, ρ and p are the density and pressure of the cosmological fluid, respectively, and u µ is the four-velocity of the fluid that in a comoving coordinate system has the form: • Since the effects of the mass of the vector fieldμ manifest on kiloparsec scales from the source, it is neglected when solving the field equations for compact objects such as black holes [28].
• G is a constant that depends on the parameter α [7]: Given these hypotheses, the action (2.1) takes the form, Variation of the latter expression with respect to g µν yields the STVG field equations: where G µν is the Einstein tensor. If we vary the action (3.4) with respect to the vector field φ µ , we obtain the dynamical equation for this field: and We propose the following metric ansatz 5 : where (t, x, θ, φ) are isotropic coordinates. Since the off-diagonal elements of the energymomentum tensor T µν = T M µν + T φ µν are zero, the G tx component of the Einstein tensor yields: If we calculate G tx and substitute in Eq. (3.9), after some algebra we get: .
If we now compare the line element given by Eq. (3.8) with the line element for a Schwarzschild STVG black hole 6 in isotropic coordinates: where dΩ 2 = dθ 2 + sin θ 2 dφ 2 , a possible form for B(t, x) is: The functions l(t) and h(t) are related to the gravitational mass M and gravitational charge Q of the source. Because we assume that both M and Q are not distributed in space but are concentrated in the singularity, l(t) and h(t) depend only on the time coordinate. By substituting Eq. (3.13) into Eq. (3.10), we derive an expression for A(t, x): (3.14) 5 In what follows we work with geometrized units G = c = 1. 6 The coordinate transformation between Schwarzschild coordinates (t, r, θ, φ) and isotropic coordinates (t, x, θ, φ) is: To determine the specific form of the functions f (t), k(t, x), l(t), and h(t), we require that A(t, x) → √ g tt in (3.12) for t = const. This yields: (3.18) In the limit l(t) → 0, h(t) → 0 (the gravitational mass and charge tend to zero), the line element should be that of a FLRW model (for simplicity we assume the spatial curvature κ = 0): Thus, the function k depends only on the temporal coordinate, and from (3.15): . (3.20) Substituting the later expression into Eqs. (3.16) and (3.18) yields: The integration constants M and Q are the gravitational mass and gravitational charge of the central inhomogeneity, respectively, whileã is associated with the scale factor a(t) of the cosmological model asã(t) = a(t).
Finally, by replacing f (t), k(t), l(t), and h(t) into expressions (3.10) and (3.13), the metric (3.8) takes the form: where the corresponding constants have been adequately restored.
In the next subsection, we prove that the metric here obtained does indeed satisfy the field equations of the theory.

Correctness of the metric
The first step in order to show that metric (3.23) satisfies the field equations of STVG given by (3.5), (3.6), and (3.7), is to compute the Einstein tensor. The non-zero components of G µ ν are: where an overdot denotes differentiation with respect to the comoving time t, κ = 8πG/c 4 , and the coefficients a 0 , a 1 and a 2 are given by: We determine the explicit form of the tensor B µν subtracting Eq. (3.25) from Eq. (3.26). After some algebraic manipulation, the non-zero components of B µν are: and B xt = −B tx . Next, we verify that Eq. (3.6) is satisfied. Since B µν is an anti-symmetric tensor: Furthermore, given that the only non-null components of B µν are B tx and B xt , two of the four equations of (3.6) are trivially satisfied. The other two remaining terms read: Thus, Eq. (3.6) holds. On the other hand, it can be easily checked that the tensor B µν with components given by expression (3.31) also satisfies Eq. (3.7).
All these lengthy calculations were necessary to prove that there is an exact solution of STVG field equations that corresponds to an inhomogeneity in an expanding universe. Our next goal is to assess the nature of this spacetime; more specifically, we first analyze whether the metric becomes singular for a certain range of coordinates and values of the parameters. Second, we compute the location of the apparent horizons and determine if they correspond to event or cosmological horizons. These features are essential to obtain a precise characterization of the spacetime and evaluate if cosmological black hole solutions are possible within the theory.

Properties of inhomogeneous expanding spacetimes in STVG
It is convenient to express the line element (3.23) in terms of the parameter α: The limits of this metric are the expected: if a ≡ 1, (4.1) reduces to the line element of a Schwarzschild-STVG black hole written in isotropic coordinates, while in the limit M → 0 Eq. (4.1) tends to the metric of a spatially flat FLRW model. For α → 0, the McVittie metric in GR is recovered.

Singularities
Singularities are a pathological feature of some solutions of the fundamental equations of a theory [37]. In GR and STVG, we can identify singular spacetime models if some physical quantity, for instance density or pressure of the fluid, or some curvature invariant is badly behaved. Thus, we begin computing the Ricci scalar for metric (4.1): Inspection of the latter equation reveals that the Ricci scalar diverges if f (t, x) = 0, that is: According to the classification of spacetime singularities introduced by Ellis and Schmidt [38], the metric possesses a scalar curvature singularity for those values of the coordinate x that satisfy Eq. (4.6). In the limit α → 0, the singular points in McVittie metric are obtained.
The singularities of the metric corresponds to the hypersurface σ(t, r) = 0, where The normal vector n a = ∇ a σ to the hypersurface and its corresponding norm are [39] : (4.8) Since in the limit x → G N (1 + α) 1/2 M/2c 2 a(t), the norm of the normal vector tends to −∞, the surface is spacelike, and consequently the singularity is spacelike. We write the Ricci scalar in terms of the energy density of the fluid and its pressure: we take the trace of Eq. (3.5), and using that T φ = T µ µ φ = 0, we get: On the other hand, we obtain an additional relation between ρ and p by subtracting Eq.
(3.24) from Eq. (3.25). The result is: Equations (4.9) and (4.10) form a system of two equations with two unknowns, ρ and p. The solution is: In the limit α → 0, the corresponding expressions for the energy density and pressure in McVittie spacetime in GR are recovered [40]. Notice that ρ is homogeneous on hypersurfaces of t constant as opposed to the pressure. From expression (4.12), we see that p diverges in the same way as the Ricci scalar. Both the energy density and the pressure have the same qualitative features as in McVittie spacetime in GR [41].

Apparent horizons
We characterize stationary black holes by the presence of event horizons. In dynamical spacetimes, however, to compute the location of the event horizon is an impossible task since we would need to know the entire spacetime manifold to future infinity. Instead, we can resort to the concept of apparent horizon. This is defined as the boundary where the convergence properties of null geodesics congruences change. The apparent horizons are located where: and θ l > 0. (4.14) Here, θ n and θ l are the expansion of the future-directed ingoing and outgoing null geodesics congruences, respectively [40]. Apparent horizons are defined quasi-locally and do not refer to the global causal structure of spacetime [40]. In spherical symmetry, the future-directed ingoing and outgoing null geodesics are radial and their tangent fields are denoted n a and l a , respectively. If the null vectors n a and l a are not affinely-parametrized, their corresponding expansions are calculated as follows: θ n = h ab ∇ a n b = g ab + l a n b + n a l b (−n c l d g cd ) ∇ a n b , (4.15) where in the later equation n b should be substituted by l b in order to calculate θ l . The tensor h ab acts as a projector onto the two-dimensional surface to which n a and l a are normal. The tangent fields to the ingoing and outgoing radial null geodesics of metric (4.1) are: , 0, 0 , (4.16) , 0, 0 . These tangents fields are such that n µ n µ = l µ l µ = 0, and g µν n µ l ν = −2. Given the vectors n a and l b , after some algebraic manipulations, the expansions for θ n and θ l take the form: (4.20) The condition θ n = 0 implies xȧ(t)g(t, x) 2 = cf (t, x), which in terms of the areal radius: can be written as: and H(t) =ȧ(t)/a(t) is the Hubble function. We can gain some insight on the nature of the apparent horizons of metric (4.1) by analyzing the limits of Eq.
The two solutions of the quadratic equation are: These are the outer (+) and inner (−) event horizons in the Schwarzschild STVG black hole [28]. Finally, if we take α → 0, Eq. (4.22) reduces to a cubic equation that locates the apparent horizons in McVittie metric in GR (see for instance Equation (4.25) in [40]). Thus, in the appropriate limits, the apparent horizons become a cosmological or a black hole event horizon, a strong hint that this metric may represent a cosmological black hole. Equation (4.22), nonetheless, has four roots. Using Descartes' rule of sign, we determine that three of them are positive and one is negative. We discard the latter because it has no physical meaning. Let us denote the three positive roots R * , R − , and R + , where R * < R − < R + . We show plots of the roots as a function of the cosmic time from Figures 1 to 4. For Figures 1 and 2, the Hubble function is that of a cosmological dust dominated background model: In Figures 3 and 4, we adopt the scale factor of the Λ Cold Dark Matter model (ΛCDM): Here, H 0 = 2.27 × 10 −18 s −1 ≈ 70 km/s Mpc, and Ω Λ,0 = 0.7 for the Hubble factor and the cosmological constant density parameter, respectively. The value of the parameter α depends on the mass of the gravitational central source. For stellar mass sources, Lopez Armengol and Romero [29] found that α < 0.1. In the case of supermassive black holes (10 7 M ≤ M ≤ 10 9 M ) the range of values are 0.03 < α < 2.47 (see for instance [42], [11], and [22]). Figures 1 and 3 correspond to a stellar mass source (we choose α = 5 × 10 −2 ) while for Figures 2 and 4 the source is a supermassive black hole (we select α = 1, and α = 2.45). In the four plots, we include the apparent horizons in McVittie spacetime in GR (α = 0) for comparison.
There are common features to all these plots: • We distinguish three apparent horizons. Two of them, R − and R + , lay in the causal future of the curvature singularity (see Eq. (4.7)). The innermost apparent horizon R * is bounded by the singularity and disconnected from the exterior geometry. In what follows, we restrict our analysis to the spacetime region that corresponds to the causal future of the curvature singularity.
• The curvature singularity (dashed line in the Figures) is present since t = 0 and its location in terms of the areal radius does not change with cosmic time. Since the surface given by Eq. (4.7) at t = 0 is in the causal past of all the spacetime events of the region of interest, we regard it as a cosmological "Big-Bang" singularity 7 .
• At early values of the cosmic time, we only have a cosmological singularity. Later on, the apparent horizons R − and R + appear together at a specific value of the cosmic time.
The horizon R + becomes larger for growing t, reaching the value of the cosmological apparent horizon in the FLRW model. Conversely, R − gets smaller for increasing values of t, and in the limit t → ∞, it gets closer and closer to the singularity.
• For larger values of the parameter α, the appearance of the apparent horizons occurs at later times. Furthermore, the value of the areal radius of the surface singularity and the horizons is higher for increasing values of α. Since R + expands forever and it tends to the cosmological apparent horizon in the FLRW model, we interpret the surface R + , t finite as a cosmological apparent horizon of the spacetime metric (3.23).
The nature of the apparent horizon R − requires some further analysis. We are particularly interested in the surface R = R − , t = ∞. In the next subsection, we show that independently of the asymptotic form of the Hubble function as long as the null energy condition is satisfied, ingoing null radial geodesics reach the null surface R = R − , t = ∞ in a finite lapse of the affine parameter. In other words, the spacetime metric to null future-oriented ingoing geodesics 8 .

The surface
Consider ingoing null radial geodesics from an initial distance R i > R − and geodesic initial velocity at that point R i < 0. We aim to show that these geodesics arrive in a finite lapse of affine parameter σ to the surface R = R − , t = ∞. More precisely: being ∆σ a finite quantity. Before starting, it is convenient to express the line element (4.1) in terms of the areal radius (see Eq. (4.21)). When doing the coordinate transformation, the algebraic manipulations are considerably simplified if you employ the relation: where we introduce r 0 and r 1 to simplify the notation: After the coordinate transformation, the line element takes the form: For radial null geodesics, ds 2 = 0 and we derive the equation: (4.32) Here, the + (−) sign corresponds to outgoing (ingoing) radial null geodesics, respectively. We focus on ingoing radial geodesics and rewrite Eq. (4.32) as: . (4.33) In the limit R → R − , H 0 is the value of the Hubble function in the limit t → ∞. To leading order, integration of the latter yields: , being γ > 0. Later on, we will use the result given by (4.35).
Next, we compute an additional radial null geodesic equation. After some algebra, we The primes denote the derivative with respect to some affine parameter σ. Notice thatḢ (t) < 0 provided the null energy condition holds. This can be proved by adding Eqs. (4.11) and (4.12): We see that if the null energy condition is satisfied, ρc 2 + p > 0, thenḢ (t) < 0. The following step is to look for an approximated formula forḢ (t): substituting an expression for the energy density of the universe: where s = 3 (1 + w) and w = p/ρ, into Eq. (4.11) and also taking into account the definition of the Hubble function,Ḣ (t) =ȧ(t)/a(t), in the limit t → ∞, and henceḢ (t) ∝ e (−sH0t) . Now, we make use of the approximation given in (4.35): Then, along null ingoing radial geodesics near the surface R = R − , t = ∞: Integration of the latter to leading order yields: (4.44) We denote by R i the initial radial velocity of an ingoing geodesics that begins at the areal radius R = R i > R − ; we also consider that R i < 0. Under this assumption, the constant G is finite and negative: (4.45) 9 The equation for R can be derived from Lagrange equations, defining the Lagrangian L = F+F−, where: After computing the Lagrange equations, it should be set F+ = 0 to derive the equations that corresponds to radial ingoing null geodesics [41].
Using Eq. (4.44), we can straightforward estimate the quantity ∆σ: If αs − 1 < 0, then (4.47) and the integral is convergent. In the case γs − 1 ≥ 0: the integral is also convergent. Consequently, the integral always remains finite and so the quantity ∆σ 10 . Hence, we have proved that ingoing radial null geodesics arrive at the surface R = R − , t = ∞ in a finite lapse of affine parameter.
The event horizon is characterized as a one way membrane: once we have crossed this surface it is physically impossible to cross it back in the opposite sense. This is precisely the case for the surface R = R − , t = ∞ providedH 0 > 0 when t → ∞.
Consider again a radial ingoing null geodesic with initial velocity R < 0. According to Eq. (4.39), and sinceḢ (t) < 0, the acceleration R is negative. This geodesic can never turn back or decrease its speed. Once the geodesic arrives at R = R − , t = ∞ in a finite lapse of affine parameter, it crosses this surface which is perfectly traversable. Recall that R = R − is a null branch of the apparent horizon, and thus constitutes a boundary where the convergence properties of null geodesics change. Right after crossing R = R − , t = ∞, the convergence properties of the geodesic are modified and it is unable to return back. Hence, the surface R = R − , t = ∞ is an event horizon, and the spacetime metric (4.1) represents a cosmological black hole.
The nature of the surface R = R − , t = ∞ whenH 0 = 0 is much more subtle. The equation for the location of the apparent horizons (4.22) can be rewritten as: (4.49) In the limit t → ∞, H(t) → 0, and f ah reduces to: where we employ the equality given by (4.28). The apparent horizons are located where f ah = 0, or equivalently where f (t, x) = 0. As shown in Section 4.1, f (t, x) = 0 identifies a singular surface of the spacetime metric. Consequently, the cosmological solution does not represent a black hole. Further investigation is needed to assess the strength of the singularity, and thus to get a better understanding of the global causal structure of the spacetime. 10 In the case the asymptotic value of H(t) vanishes, the demonstration exists and is quite alike.

Discussion
Until now, we have explored the properties of the solution for a limited range of values of the parameter α. In what follows, we remove such restriction and allow α to freely move in the interval 0 < α < ∞.
As already mentioned in Section 4.1, the location of the singularity is independent of the cosmic time. We write Eq. (4.7) in terms of the areal radius: The function R(α) is strictly increasing: the higher the value of α, the larger the areal radius of the singularity. Also notice that there is no value for α such that the singularity can be avoided, implying that there are no regular cosmological black hole solutions in the theory. The location of the apparent horizons as a function of the parameter α is depicted in Figure 5, for a fixed value of the cosmic time. As before, the dashed line marks the location of the singularity. In the interval 0 < α <α, we identify three apparent horizons: an inner horizon R * that lies beyond the singularity (and hence is not part of the spacetime), a black hole apparent horizon R − , and a cosmological apparent horizon R + . As α gets closer toα, R − increases while R + becomes smaller. For α =α both horizons, R − and R + , become one.
Higher values of α implies an augmented gravitational constant. In STVG the gravitational field is stronger that in GR; the central source drags the cosmological horizon while the black hole apparent horizon enlarges.
If α >α, the apparent horizons disappear and a naked singularity is left behind. Accepting the validity of the cosmic censorship conjecture [43], we see that restrictions can be imposed on the values of the parameter α such that solutions that contain naked singularities are not allowed in the theory. The constraint on α changes for different values of the cosmic time (the coefficients of Eq. (4.22) depend on the Hubble function H(t)). The latter implies that the permitted values of α do not only depend on the mass of the central source but also on the cosmic epoch of the universe.

Conclusions
In this work we derive the first exact solution of STVG field equations that represents an inhomogeneity in an expanding universe. When the Hubble factor is positive at late cosmic times, we prove that the metric describes a black hole immersed in a cosmological background.
The spacetime presents a spacelike singular surface where the pressure of the cosmological fluid diverges, a feature that is common to McVittie metric in GR. We also show that there is no value of the parameter α of the theory such that the singularity can be avoided. This result implies that there are no regular cosmological black hole solutions in STVG.
The metric has two apparent horizons: an inner horizon and an outer horizon that correspond to an event and cosmological horizon for the black hole case. As the value of the parameter α increases, the size of the horizons enlarges as well as the areal radius that locates the singular surface. enlarges. The inner horizon approaches the singularity while the outer one tends to the cosmological horizon in the FLRW model. For a given value of the cosmic time, there is a limited range of values of α such that the solution exhibits an inner and an outer apparent horizon. Beyond this range, both horizons merge and finally disappear leaving behind a naked singularity. If we assume the validity of the cosmic censorship conjecture, we see that only some values of α are allowed. Thus, the value of α is not only dependent on the mass of the central source but on the cosmic epoch. This result should be taken into account when modeling the evolution of the structure and the dynamics of astrophysical systems through cosmic time.
This work is a first step towards a better understanding of cosmological black holes in STVG. There are several issues that remain unexplored; for instance, the strength of the spacelike surface singularity, the nature of the cosmological solution when the Hubble factor is zero at late times, the dynamics of particles in this spacetime, just to mention some. The fact that STVG admits cosmological black hole solutions is yet another positive indicator that the theory offers a suitable classical description of the various manifestations of gravity.