Equatorial orbits and imaging of hairy cubic Galileon black holes

Null and timelike equatorial orbits are investigated in a family of hairy black holes in the cubic Galileon theory. These include rotating generalizations of static black hole metrics supporting a time-dependent scalar field. Depending on the coupling and rotation, the properties of the geodesics expectedly deviate from general relativity. In particular, it is found that stable circular geodesics only exist below a critical coupling, which is related to the existence of an outermost stable circular orbit. Focusing on the strong-field region, images of an accretion disk are also produced to highlight tendencies that would constrain the model given further accurate observations of supermassive black holes.

It must be noted though that unequivocally excluding a given model is a complex endeavour in most cases. Current instrumental limitations and extremely large parameter spaces (describing e.g. the emission flow or sources of noise) require to rely on simplifying prescriptions [32] to realize and interpret observations, such as the EHT reconstructed image. This leaves too much uncertainty to draw definitive conclusions today, and extensive studies are needed to explore significant parts of the parameter spaces and guess which observations could be decisive in discriminating some given models. But the understanding of astrophysical black holes gradually progresses by improving current instruments and analysis tools, and developing ideas for future enhanced observations.
In this context, the present paper focuses on the characteristic features and preliminary constraints arising from the strong field observables of a family of black hole spacetimes within the cubic Galileon theory. The ("covariant generalized") Galileons are scalar-tensor theories which coincide with Horndeski theories in four dimensions, meaning that they are the most general scalar-tensor theories leading to second-order field equations [33,34,35,36,37,38]. Galileons thus provide a relevant framework to search for observable deviations from general relativity (GR), as most signatures of alternative theories are expected to be well described by scalar-tensor theories at least in some effective range. In particular, the cubic Galileon emerges from effective formulations of higher dimensional theories, either in the decoupling limit of braneworld models such as the popular 5-dimensional Dvali-Gabadadze-Porrati (DGP) model [39,40,34,41,42], or from Kaluza-Klein compactification of Lovelock theory [43,44]. Explicitly, its vacuum action writes where (∂φ) 2 ≡ ∇ µ φ∇ µ φ and ζ, η and γ are coupling constants. It is the simplest of Galileons with higher order derivatives, and it is compatible with the observed speed of gravitational waves [45,46,47,48]. As a well-motivated, consistent theory, the cubic Galileon has been investigated in various contexts, from laboratory tests [49] to cosmology [50,51,52,53,54]. Besides current observations of supermassive black holes, further interest in the characteristics of cubic Galileon black holes comes from the fact that, together with most shift-symmetric ‡ Horndeski theories, the cubic Galileon is subject to a "no-scalar-hair" theorem: in the asymptotically flat framewok, the only static, spherically symmetric black hole metric and scalar field are the Schwarzschild metric together with a trivial scalar field [55,56,57,58] (see also [59] for an extension to slow rotation and [60] for stars). Consequently, modified gravity effects can only occur in systems breaking one of these hypotheses. Indeed, one does obtain non-GR metrics coupled to non-trivial scalar hair when enforcing all hypotheses except the stationarity of the scalar field. Such minimal violation of the hypotheses is possible when the scalar field features a linear time-dependence [61] where q is a non-zero constant and Ψ is time-independent. ‡ Shifting the scalar field by a constant (φ → φ + constant) preserves action (1).
The shift-symmetry is what makes a linear time-dependence compatible with a static and spherically symmetric metric. Such configuration was considered in various contexts such as cosmology, and the linear time-dependence was physically interpreted as a first order approximation to a slowly evolving scalar field [62,63,53,54]. Besides, using ansatz (2) (with Ψ depending not only on the radial but also the angular coordinate), previous numerical work [64] produced rotating generalizations of hairy static and spherically symmetric solutions derived in the cubic Galileon theory [65]. At the level of the metric, these rotating black holes significantly deviate from Kerr spacetime, implying possibly observable modified gravity effects in black hole environments.
Such hairy configurations were constructed as circular spacetimes, meaning that they were assumed to admit a quasi-isotropic coordinate system with respect to which the line element writes Circular spacetimes represent a large subclass of stationary and axisymmetric spacetimes, and their quasi-isotropic coordinates are well suited to study geodesics. Yet compatibility of a circular metric (3) with a linear time dependence (2) is exact only in the non-rotating case, while errors arise as rotation increases and could become significant at high rotation. This is why these spacetimes are only considered at low and moderate rotation such that the errors on the solutions are negligible [64]. Despite such restriction, these configurations allow to observe non-perturbative effects of rotation (examples of which are still not so abundant in modified gravity, although relevant for astrophysical black holes which are expected to rotate).
Besides, the metrics were constructed to be asymptotically flat, as it is a natural hypothesis of the no-scalar-hair theorem. Yet this is realized by taking η = Λ = 0 in the cubic Galileon model (1), leaving the scalar field ruled by the non-standard "DGP term" (∂φ) 2 φ. This induces a non-standard asymptotic behaviour of the metric which yields a vanishing Komar mass at infinity. This theoretical fact is unusual, but astrophysically, it does not forbid the hole to feature an effective mass on finite distances. More precisely, we will see that there do exist stable orbits, yet only up to an OSCO (outermost stable circular orbit). OSCO's are not so unusual even in GR in presence of a positive cosmological constant, and they can emerge on short enough scales to be astrophysically relevant (for instance, using the Schwarzschild-de Sitter metric outside of a galaxy, the OSCO is of the order of the inter-galactic distance [66]). In the cubic Galileon case, studying the geodesics of the asymptotically flat hairy black holes will show that the OSCO actually emerges on an even shorter scale when the remaining parameters of the model are chosen so as to generate non-negligible metric deviations in the strong-field region. The asymptotically flat model itself is then strongly constrained by requiring stability of all the orbits that would be dominantly ruled by a central supermassive black hole.
Note however that non-zero couplings η and Λ lead to Schwarzschild-de Sitter asymptotics in the static case [65], which generically restore a large OSCO compatible with observations. Despite different asymptotic properties, models with non-zero η and Λ, and asymptotically flat solutions, might present analogous dependencies on the DGP coupling and on rotation, and share common observable characteristics in the strong field region. This is why images of an accretion disk orbiting the asymptotically flat black holes are also presented below, while future work will construct asymptotically de Sitter solutions to assess such similarities (and more generally to identify degeneracies of the observables within the cubic Galileon model).
Besides, if some strong field characteristics of the asymptotically flat configurations turned out particularly interesting without being preserved for any non-zero η and Λ, one might deal with the unusual large distance properties of the asymptotically flat model by relying on convenient mechanisms or fields, such as a second, "screened" scalar field χ. Screened scalar fields are commonly invoked in cosmology, in different realizations of massive gravity and hence in Galileon theories [41] to recover the successful predictions of GR on short scales, e.g. solar system scales, while providing new relevant cosmological phenomenology in regard of dark energy. In standard screening mechanisms (Chameleon [67,68], Symmetron [69,70], Vainshtein [71,72,73,42]), the mass of χ or its coupling to matter effectively decrease in regions of high matter density. Such processes are realized through non-standard kinetic terms and/or appropriate couplings to standard model matter, such as Dirac spinors, Higgs scalars or other fundamental fields [69]. Although the idea has not been considered elsewhere, these mechanisms might be adapted to the present asymptotically flat case, e.g. by introducing analogous interactions that suitably couple χ to the black hole hair φ (rather than to standard matter): χ would be screened where φ adopts its strong field profile, i.e. the immediate vicinity of the black hole, while it would modify the geometry on larger distances, and in particular the location of the OSCO. Thus, the asymptotically flat case may provide an effective description of the strong field region of more complete models that feature more standard weak field properties thanks to a standard kinetic term for φ, a cosmological constant Λ, or alternative fields or mechanisms.
The plan of the article is as follows. The properties of the equatorial timelike circular geodesics and photon rings (location, stability, deviation from Kerr spacetime) are studied in the static, spherically symmetric case in section 2.1, and the rotating case in section 2.2. This leads to to compute the images of an accretion disk orbiting the static, spherically symmetric black holes in section 3.3, and the rotating black holes in section 3.3. Notations and general results on equatorial geodesics in quasi-isotropic coordinates are summarized in appendix Appendix A.

Static and spherically symmetric case
To study the geodesics of the cubic Galileon static and spherically symmetric black holes obtained in [64,65], the procedure is to first characterize the circular geodesics.     [64]), so that D > 0 implies that the denominator in (A.11) is positive. Therefore, V + > 0 (prograde orbits) and V − = −V + < 0 (retrograde orbits).
As detailed in appendix Appendix A, regions of positive discriminant D (A.10) are first checked on figure 1a for different values ofγ = q 3 r H γ/ζ, where r H is the radial coordinate of the event horizon §. Function D appears positive everywhere down to the horizon, where it diverges because of division by the lapse N which cancels at the horizon. Therefore, circular geodesics a priori exist everywhere for all couplings, but they necessarily become superluminal near the horizon according to (A.11). This is what figure 1b confirms: for each couplingγ, velocity diverges at the horizon so that there exists a photon ring (marked with a vertical line from 0 to 1), only beyond which timelike circular geodesics exist. Although D does not vary much with coupling on figure 1a, the velocities more strongly depend onγ because function B in the denominator of (A.11) does vary (see figure 2(b) of [64], knowing that B = A everywhere in spherical symmetry). More precisely, at fixed radius, the velocity of the circular geodesic decreases with increasing coupling. As a consequence, the photon ring gets closer to the horizon asγ increases.
These results are related to the following facts mentioned in the introduction. The metric functions N and A = B converge faster to Minkowski at infinity asγ increases [64]. Therefore, at fixed radius away from the strong-field region, gravitation § Recall thatγ is the only remaining dimensionless coupling parametrizing deviations because η and Λ are set to 0; the values ofγ considered in figure 1a are picked in the range leading to non-negligible metric deviations from GR [64].    gets naively weaker asγ increases, so that the velocity of the circular geodesic must be smaller. In addition, for anyγ = 0, convergence to Minkowski spacetime is always much faster than that of Schwarzschild spacetime: N and A = B converge to 1 as 1/r 4 rather than 1/r, yielding a vanishing Komar mass at infinity [64]. As a result, velocities given by (A.11) converge to zero like r −α/2 with α = 1 in Schwarzschild spacetime and α = 4 in Galileon spacetimes.
Such asymptotic behaviours are highlighted in figures 2. In all cases, the Lorentz factor displayed in figure 2a logically converges to 1. However, according to (A.12), the Killing angular momentum per unit massL displayed in figure 2b behaves like rV r 1−α/2 , hence the divergence in Schwarzschild spacetime and convergence to zero for anyγ = 0 . Finally,Ē = ΓN converges to 1 in all cases on figure 2c, which will The numerical solutions contain information at infinity confirming this fact even for small couplings 20 40 60 80r hold true in the rotating case since function ω will converge to 0 like 1/r 3 regardless of whetherγ is zero or not. At the photon ring, all the kinematic quantities displayed in figure 2 naturally diverge. Figure 3 assesses the stability of circular orbits for various couplings, based on the functions V ± given by (A.14). As explained in appendix Appendix A, their sign at a given radial coordinate r 0 is the same as V r 0 , 1,Ē ± (r 0 ),L ± (r 0 ) respectively. It appears that for any non-zeroγ, both an innermost and an outermost stable circular orbit (ISCO and OSCO) exist. They respectively correspond to the smallest and greatest r 0 such that V ± (r 0 ) = 0. Since V ± globally decreases asγ increases, the ISCO radius increases while the OSCO decreases from infinity (where it is formally located in the Schwarzschild caseγ = 0). In particular, the photon ring is always unstable because its location decreases whenγ increases (figure 1b), so that it remains below the ISCO as in Schwarzschild spacetime.
The ISCO and OSCO eventually merge for a critical couplingγ c 2.2 × 10 −2 , beyond which no stable circular orbit exists anywhere. Therefore, the mere existence of stars orbiting black holes such as Sgr A* in a seemingly stable way requiresγ γ c for the present static Galileon black hole to be viable. The existence of a particularly close OSCO further constrains the model since e.g. the well-known star S2 lies beyond 2500M , where M 4×10 6 M is the observed mass of Sgr A*. Although its orbit is non-circular, it indicates that a stable circular orbit exists between its apsides, and hence the OSCO must lie beyond. For instance, this further requiresγ 10 −3 . Note though that a perturbation δ away from a circular orbit at some radius r 0 > r OSCO obeys equation likeγ = 10 −3 whose convergence to zero becomes apparent very far from the horizon.
so that the instability timescale is Since V asymptotically goes to zero as r −6 , τ diverges and rapidly becomes larger than the age of the universe (e.g. around 30r OSCO for r H of the order of the Schwarzschild radius of Sgr A* and forγ ∼ 10 −2 ). Although this does not render the model viable ¶, such instability is weaker than in the standard cases featuring an OSCO (in Schwarzschild-de Sitter spacetime, τ converges to a finite value, fixed by the cosmological constant). This goes along with the fact that this close OSCO is a singular artefact of the marginal combination η = Λ = 0, and corroborates the idea that a much larger stability region would be restored with any additional Lagrangian terms or mechanism mentioned in the introduction. With this in mind, the next section focuses on the effects of rotation on the orbits.

Rotating case
Rotation breaks spherical symmetry so that the "+" and "-" quantities are no longer equal or opposite, as shown in figures 4 (in which solid lines correspond to the analogue quantities in Kerr spacetime). Yet all these quantities have the same behaviour at the boundaries as in the static, spherically symmetric case. Besides, figure 4a shows that V + > 0 and V − < 0 still hold everywhere. However V + = −V − no longer does, so that there exist a prograde photon ring and a distinct retrograde one for each angular velocityΩ H = r H Ω H . As in Kerr spacetime, prograde and retrograde velocities decrease asΩ H increases, so that the prograde (resp. retrograde) ring radius decreases (resp. increases). The dependence onΩ H yet seems stronger in Kerr spacetime, meaning e.g. that the prograde ring radius decreases faster than for any non-zeroγ. Since the photon ring of the static, spherically symmetric Galileon spacetime is below that of Schwarzschild spacetime, the relative positions of the Kerr and Galileon prograde rings are inverted for someΩ H ( 0.03 forγ = 10 −2 ). On the contrary, the Kerr retrograde ring grows away from its Galileon counterpart.
The fact that the dependence onΩ H is qualitatively the same in Kerr and Galileon spacetimes, but stronger in the former, also applies toL (figure 4b),Ē (figures 4c and 4d) and V ± (figures 4e and 4f). Besides, V + (resp. V − ) globally increases (resp. decreases) asΩ H increases. Therefore, both the prograde ISCO and retrograde OSCO (resp. prograde OSCO and retrograde ISCO) radii decrease (resp. increase) with rotation. Interestingly, since the ISCO radius of static, spherically symmetric Galileon black holes is beyond Schwarzschild's ISCO, and Kerr's retrograde ISCO increases faster withΩ H , the relative positions of the Kerr and Galileon retrograde ISCO's are inverted for someΩ H ( 0.06 forγ = 10 −2 ). Furthermore, sufficiently highΩ H makes it possible for V + to become positive even forγ greater than the critical couplingγ c 2.2 × 10 −2 ; this is illustrated in figure 5a          forγ = 1 >γ c . Therefore, for each couplingγ >γ c , there is a minimal angular velocitiy beyond which stable orbits reappear, yet only prograde ones. On the contrary, asΩ H increases, the Galileon retrograde ISCO and OSCO eventually merge for a critical angular velocityΩ c H ( 0.09 forγ = 10 −2 ), beyond which no stable retrograde orbit exists anywhere. Therefore, the fact that stars stably orbit Sgr A* in both directions + leads to even tighter constraints thanγ γ c if Sgr A* is modeled as a rotating black hole.

Principle of ray-tracing
In the present section, images of an accretion disk orbiting the black holes are computed numerically. This is again motivated by the idea that the present model may capture the significant characteristics of more complete, better-behaved models. Computations are performed by the free, open-source ray-tracing code GYOTO [74], which features an efficient approach to integrate the geodesic equations from the knowledge of the 3+1 quantities decomposing a numerical metric [75]. In our case, the shift β and spatial metric g 3 corresponding to the quasi-isotropic metric (3) are Images are computed in the following way. An explicit model of accretion flow is set around the black hole * . A telescope set in the numerical metric mimicks the + The spin direction of Sgr A* is unknown. * Rough estimates, confirmed by simple exact models of accretion disks, show that the gravitational influence of an accretion disk is usually completely negligible with respect to the black hole. Thus, observing wavelength (1.3 mm), the distance (16.9 Mpc), field of view (120 µas) and orientation of the Event Horizon Telescope with respect to M87* (the black hole being set at the origin and the disk lying in the equatorial plane θ = 90 • , the colatitude of the Earth is θ = 160 • , while the vertical axis of the screen of the EHT is rotated by 110 • clockwise from the projection on the screen of the spin axis of the disk). Each pixel of its focal screen corresponds to a spatial direction, which uniquely defines the initial tangent vector of a null affinely parametrized geodesic. The latter is integrated backwards in time until a stopping condition is met, e.g. the photon gets too close to the event horizon, or definitely leaves the strong field region. Otherwise, every time the geodesic crosses the accretion disk, the radiative transfer equations ruling the specific intensity are integrated along the segment lying within the disk. The cumulated specific intensity is eventually plotted on the initial pixel.
Yet, determining the nature and properties of a compact object based on the image of its accretion flow is a very degenerate inverse problem [19,26,76]. This is for instance evidenced in reference [13] in which the same model of accretion disk is set around a boson star and a black hole: the deviations between the resulting images are very subtle although the natures of the accreting objects are very different. Furthermore, the resolution of present and future instruments like the EHT is limited, making it even harder to distinguish subtle features .
Then, the purpose of numerical images is not systematically to check whether the image constructed by the EHT [8] can be reproduced for different accreting compact objects. This indeed requires costly general relativistic magnetohydrodynamics (GRMHD) simulations, together with a model of the EHT itself. Instead, strong efforts are made to propose fairly simple and yet realistic models of accretion disks [77,78,79,9]. In particular, such models are assumed to be good approximations of stable steady solutions of the GRMHD equations. Comparing the resulting images for different compact objects provides a more efficient and still relevant method to evaluate how degenerate the problem is. The hope is that such an approach should help isolating the causes of differences between images, e.g. being able to guess the nature and amplitudes of the modifications that result from changing the accretion model or the theory used to describe the whole system.
As a result, a simple model of accretion disk, recently introduced in [9], is used in the sections below. Like Sgr A*, supermassive black hole M87* features a very low-luminosity accretion flow, revealing an inefficient radiative cooling and hence a high temperature. It is consistently modeled as a low accretion rate, geometrically thick, optically thin disk † [80]. Besides these properties, only the thermal synchrotron the vacuum black hole metrics are still valid in presence of an accretion disk. See section 6.5 of [1] for quantitative arguments.
Regarding the particular problem of distinguishing boson stars from black holes, see however [14] for possibly detectable deviations arising from dynamical effects revealed by 3D general relativistic magnetohydrodynamics simulations. † An accretion disk is geometrically (resp. optically) thin when the opening angle (resp. optical depth) is smaller than 1. It is geometrically (resp. optically) thick otherwise. emission is computed, following a method exposed in [81]. In the end, the complete model is described by very few input parameters: the opening angle and inner edge of the disk (which is set at the ISCO in our case), the magnetization parameter (which determines the ambient magnetic field strength), the electron density and temperature at the inner edge (which determine the density and temperature profiles).

Static and spherically symmetric case
To connect with section 2, the black holes on figure 6 all have the same radius, namely the Schwarzschild radius of M87* (in quasi-isotropic coordinates, i.e. r H = M M87* /2 where M M87* 6.5 × 10 9 M ). It is actually more natural to manipulate the horizon radius than any notion of mass: because of the non-Schwarzschild asymptotics of the Galileon black holes recalled in section 2.1, there is no relevant mass parameter that can be extracted from star trajectories (through Kepler's law) or surface integrals (such as the Komar or Arnowitt-Deser-Misner masses) in the weak field region, since its value could not be compared in any meaningful way to that of a Schwarzschild black hole (recall for instance that these Galileon black holes feature a vanishing Komar mass at spatial infinity, and yet admit stable orbits). It is thus equally legitimate to make comparisons based on parameters specific to the strong field region such as the horizon, lightring and ISCO radii, and the characteristics of the images described below.
In the non-rotating case, images 6a and 6b compare the Schwarzschild limitγ = 0 to a static and spherically symmetric Galileon black hole below the critical coupling. In both cases, the ISCO radius is read from figure 3. The asymmetries within each image come from the configuration of the EHT with respect to the accretion disk. More precisely, the 110 • clockwise rotation of the vertical axis of the screen from the projected spin axis of the disk explains the position on the images of the brighter spot that results from the relativistic beaming and blueshift affecting the part of the disk rotating towards the screen. Besides, the inside luminous ring corresponds to the secondary and higher order images, which aymptotically accumulate in the direction of the lightring. Rather than being centered within the primary image of the disk, this ring appears shifted towards the top of the spin axis because of the inclination angle θ = 160 • of the EHT observer (the disk is almost seen from below).
This ring gets smaller with respect to the primary image asγ increases because the ISCO radius increases (figure 3) while the lightring decreases ( figure 1b). Explicitly, in Schwarzschild case, the internal diameter of the primary image is D 50 µas while the diameter of the secondary ring is d 42 µas, which is compatible with the EHT image [8]. In the Galileon case, D 46 µas while d 34 µas. Besides, the image is globally brighter asγ increases. This is consistent with the fact that asymptotic convergence to Minkowski is faster asγ increases (as recalled in section 2.1). In other terms, the strong field region shrinks asγ increases, so that most light rays undergo a weaker gravitational redshift. Therefore, if one considered a Galileon black hole with e.g. a 1.12 times   greater radius, so as to fit the internal and secondary diameters close enough to the Schwarzschild values (the deviation on their ratio being possibly undetectable by the resolution of the EHT ‡), the corresponding image would get even brighter than figure 6b. One could accordingly reduce the density parameter to avoid tension with the observed luminosity. Yet the fact would remain that the ISCO radius would be even greater than it already is with respect to its Schwarzschild analogue when the horizon radii are the same (figure 3). It is hardly conceivable that the ISCO radius of M87* could be estimated by other means than EHT-type observations, so that no incompatiblity related to this parameter could be exhibited. On the other hand, if images of Sgr A* were obtained, tensions about the ISCO location should arise based on the high precision astrometric observations made by the instrument GRAVITY which detected "flares" [5] close to Schwarzschild ISCO, knowing that such bright spots are expected to materialize near the inner edge of the accretion disk of Sgr A*.

Rotating case
Images 6c and 6d compare a Kerr black hole to a rotating Galileon black hole sharing the same angular velocity. In both cases, the inner edge of the disk is set at the prograde ‡ Actually, in order to reproduce even more closely the diameters, and hence their ratio, one could fine tune the dimension of the internal diameter of the primary image by considering a more realistic scenario involving accreting matter on non-circular trajectories below the ISCO. In this sense, the secondary ring is more reliable as an observable of the gravitational field as it weakly depends on the boundaries and physical properties of the accretion process.
ISCO, which is read from figure 4e. Besides, the inside ring intersects the primary image, meaning that the ISCO is close enough to the light ring to allow photons to cross the disk at least a second time. This is consistent with the fact that the prograde ISCO decreases (figure 4e) faster than the prograde lightring (figure 4a) asΩ H increases. As in the static case, the image is globally brighter asγ increases. But at a fixedγ, the images get darker asΩ H increases because the strong field region expectedly expands with rotation. As illustrated by figure 7, such tendencies leave room for a black hole with moderately low rotation, and radius greater than M M87* /2, to closely mimick Schwarzschild image 6a in terms of dimensions and brightness. Yet regarding the ISCO radius, the fit is expectedly not so good with respect to Schwarzschild (9% and 50% deviation on the prograde and retrograde ISCO respectively). Of course, clearing away degenerate tendencies and exhibiting significant incompatibilities would only be possible if images and more precise spin measurements of Sgr A* were available.

Conclusion
Investigating the geodesics of asymptotically flat cubic Galileon black holes exhibits non-viable characteristics when the couplingγ is in a range leading to non-negligible metric deviations from GR. As such, the model would be constrainted to a negligible couplingγ for the OSCO to lie far enough to be compatible with distant stars observed around Sgr A*. Yet it is thought that constraints only based on the OSCO might be dismissed by restoring terms in Lagrangian (1), or invoking other appropriate fields, which would preserve the tendencies and strong-field characteristics of the asymptotically flat model. This is why images of accretion disks have also been computed, as these allow to probe the close environment of a supermassive black hole. The dependence of this observable on coupling and rotation may lead to deviations from GR in terms of global brightness and relative dimensions of the luminous structures. Yet these may often be compensated by adjusting density and the horizon radius.
Proper constraints will hopefully arise from future images and precise measurements of the spin and ISCO location of Sgr A*. The latter being one of the main targets of the EHT, knowledge about this system will keep increasing in the coming years. Besides probing the inner accretion flow [82,83], the observations from the EHT will contribute to evaluating its spin (while former estimations are based e.g. on quasi-peridodic oscillations of the radio emission [84]). Indeed, EHT data was already used to establish bounds on the spin of M87* [85,86], and this will improve with further observations from the EHT [31]. Information on the spin may also be extracted from instruments of the Very Large Telescope such as the interferometer GRAVITY which can monitor faint stars very close to Sgr A* [87]. From a broader point of view, GRAVITY also contributes to better understanding the astrophysics of massive black holes [88], and such general knowledge will also improve in future decades thanks to space-based gravitational-wave observatory LISA [89,90], which will extract precious and precise information from stellar-mass compact objects spiralling around massive black holes. On an even longer term, reference [91] proposes a space-based very long baseline interferometry experiment which would characterize the ring-shaped structures in images of black hole accretion flows with extremely high precision, measuring the diameter size to 0.04% accuracy, while the EHT is limited to approximately 10% accuracy. This highlights again the fact that constraining theories of gravity with massive black holes is a long-term enterprise, as determining the parameters of black holes and their accretion flow is a degenerate problem conditioned by instrumental accuracy and modeling limitations [92,93,94]. It is in particular delicate to predict when and which opportune combination of observations will provide inflexible constraints or could at least strongly disfavour a model.
The models considered in the present paper are not free of simplifications either, both in the metric and the model of disk, which could be addressed in later work. As explained in the introduction, the circular metric (3) is only exact in the static case. Reproducing the rotating solutions in a noncircular framework would allow to reach accurate rapidly rotating solutions (possibly up to an extremal case). Such general framework would also be used to construct the asymptotically de Sitter configurations corresponding to non-zero η and Λ, and perform similar study and comparisons. Furthermore, in order to make precise quantitative comparisons between actual images and numerical predictions in a consistent way, more realistic models of disk such as ion tori could be considered. The latter are also geometrically thick and optically thin structures, yet featuring more complex density and temperature profiles derived from first principles, as well as isotropic [77] or toroidal [79] magnetic fields, hence allowing to explore other parts of the parameter space describing accretion disks. Thus, these four conservation equations are four necessary conditions for a curve C to describe an equatorial trajectory of a free particle. For non-circular orbits (i.e.ṙ = 0 almost everywhere), they are sufficient and have a unique solution provided that the initial sign ofṙ is fixed and the initial coordinates (t 0 , r 0 , π/2, φ 0 ), the Killing energy E, the Killing angular momentum L and the mass m satisfy V(r 0 , m, E, L) < 0 and E − ω(r 0 )L > 0 (to guarantee that the trajectory is initially causal, i.e. future-oriented). For a circular geodesic at radial coordinate r 0 , one hasṙ = 0 so that V(r 0 , m, E, L) = 0. Yet, as detailed in [97], the additional condition where E is the energy measured by the ZAMO: for a massive particle, E = Γm where Γ is the Lorentz factor of the particle with respect to the ZAMO, while for a massless particle, E = hν where ν is the frequency measured by the ZAMO. Finally, the radial equation (A.6) expectedly provides a stability criteria for circular geodesics based on convexity: for any perturbation to be bounded in some neighbourhood of a geodesic at r, V (r, m, E, L) must be positive. Since the values of E and L for a circular geodesic at r are necessarily given by relations (A.12) and (A.13), one only has to study the sign of the two functions V ± : r → V (r, m, E ± (r), L ± (r)) (A.14) on the set on which the discriminant D is non-negative. Actually, the expressions V ± (r) are homogeneous with respect to E, so that their sign do not depend on Γ ± m in the massive case nor on hν in the massless case. Therefore, the stability of the causal circular geodesics only depends on the sign of the two functions (A.14) and corresponds to massive particles where V ± (r) ∈ (−1, 1) and massless ones where V ± (r) = ± 1, regardless of whether the expressions used for E ± and L ± apply to a massive or a massless particle. These are therefore the two functions that are plotted in section 2 to study the stability of circular geodesics in the cubic Galileon spacetimes.