Chaotic Solutions and Black Hole Shadow in $f(R)$ gravity

We discuss the emergence of black hole shadow and photon-sphere in the context of $f(R)$ gravity. It is shown that the shadow is exponentially sensitive to linear instabilities of metric coming from some $f(R)$ solutions. Thus, the instabilities of photon circular trajectories, delimiting the black hole photon-sphere, are double exponentialized. Specifically we individuate two Lyapunov exponents, rather than only one, related to two different sources of chaos in geodesic orbits as a sort of butterfly effect. Such a result violates the black hole chaos bound proposed by Maldacena, Shenker and Stanford for General Relativity. We also explore the impact of the black hole metric instabilities in $f(R)$ gravity on the quasi-normal modes. In the framework of Extended Theories of Gravity, our analysis suggests a new paradigm to deal with black hole shadow and gravitational waves observations coming from black hole merging in the ringdown phase.


INTRODUCTION
One of the most fascinating aspects of black hole (BH) physics is that, below a critical distance, also light gets trapped in circular orbits. In other words, no any signal can escape from the BH gravitational attraction below a radial threshold. The photons orbit around the "BH design" giving rise to the so dubbed BH photon-sphere. Indeed, the BH appears as a black hole shadow for any distant observers, surrounded by a luminous photon ring. In General Relativity (GR), the photon-sphere is limited by a typical radius r C which is related to the Schwarzschild radius r S as 2r C = 3r S [1,2]. Such a circular null-like geodesic is exponentially unstable against arbitrary small geodesic perturbations, as a manifestation of the so-called butterfly effect. In other words, here there is the emergence of chaos: any deviation from the critical photon-sphere trajectory leads to instabilities, exponentially growing in time. Chaos is ubiquitously appearing in many different aspects of Nature, notoriously including atmospheric physics, finance, biology and so on. Thus, it is not surprising to find "chaos footprints" in BH physics. Specifically, a Lyapunov exponent, related to the BH shadow, parametrically controls the exponential chaos growth in time. Maldacena, Shenker, and Stanford (MSS) conjectured that such instabilities cannot grow more than exponentially, with a Lyapunov exponent which cannot exceed a critical value directly proportional the Bekeinstein-Hawking temperature [3] 1 . In general, chaos can be "diagnosed" from analyzing the out-of-time-order (OTO) correlators related to the commutators of time-separated operators [3]. In chaotic systems, the OTO Operators have exponential instabilities controlled by the Lyapunov exponent 2 . It is also interesting to note that the circular unstable photon-sphere is related to the characteristic modes of BH -the so-called quasi-normal modes (QNMs). More precisely, the BH Lyapunov exponent is proportional to the QNM decays after the BH formation. This suggests that BH shadow issues are also related to possible tests of the BH horizon from gravitational waves emitted from the BH merging during the ring-down phase, when a new BH has formed undergoing to relaxation. It is worth stressing how QNMs are crucially important for the gravitational waves physics after LIGO/VIRGO observations [8][9][10][11].
On the other hand, the above considerations may be significantly improved considering Extended Theories of Gravity [12]. In principle, any extensions of GR may lead to revisit our standard conceptions of BH physics, including stability issues and BH thermodynamics. The possibility that the standard Einstein-Hilbert (EH) action may be embedded and corrected extending the gravitational action considering further curvature (and torsion invariants), besides the Ricci curvature scalar R, is one of the most debated issue of contemporary gravitational physics and cosmology. Indeed, within the vast landscape of possible theories of gravity, which are theoretically self-consistent and pass astrophysical and cosmological tests (specifically the Solar System ones), there are several models leading to new predictions and solutions in BH physics.
First of all, the Birkhoff theorem for spherically symmetric BHs can be violated in many gravity theories beyond EH. Let us remind that, in GR, the Birkhoff theorem states that any spherically symmetric solution of Einstein field equations is static and stationary. Despite of this statement, many spherically symmetric solutions of Extended Gravity are unstable under linear metric perturbations [13-21, 23, 24]. Linear metric instabilities were found in many theories, including f (R) gravity, for extremal Schwarzschild-de Sitter and Reissner-Nordstrom BHs [26,27]. Furthermore, several BH solutions have been found and studied in the context of f (R) gravity, see [28][29][30][31][32][33][34][35].
In this paper, we are going to analyze the BH shadow, the chaos photon-sphere instabilities and the QNMs in f (R) gravity. We find that the OTO correlators are double exponentially sensitive for two characteristic Lyapunov exponents rather than the only photon-sphere one. Such a result violates the MSS bound. On the other hand, the scalaron field instabilities, related to f (R), must be fine-tuned against an exponential function of time. If not so, an obvious contradiction with lensing and gravitational waves data would rule out f (R) gravity and similar models that seem retained from observations [25] (see also Ref. [26] for an analysis of Active Galactic Nuclei and lensing data).
This approach may suggest several opportunities to tests f (R) gravity and, in general, Extended Theories of Gravity, by BH physics. First of all, the observations by the the Event Horizon Telescope (EHT) collaboration of the M87 * BH shadow [36][37][38] can be used to probe possible anomalous growth or decrease of the BH area. In general, any supermassive compact object can be adopted to test gravity and forthcoming observations towards our Galactic center Sgr A * , as recently reported in [39], will constitute a formidable laboratory in this perspective.
Specifically, f (R) gravity can predict BHs with an event horizon increasing in radius in time. Furthermore horizonless solutions can be achieved in this context [40]. Finally, the new Lyapunov exponent emerging in f (R) gravity is related to the quasi normal modes (QNMs) of BH in the ringdown phase. Thus, metric instabilities coming from f (R) gravity, beyond the Einstein theory, can be, in principle, tested, and then confirmed or rejected, from gravitational waves observations after the BH or neutron star merging events to a final BH state. This phenomenology could be of high interest for LIGO/VIRGO [8][9][10][11] and KAGRA [41] collaborations.
To fix the ideas, let us consider OTO correlators and BH chaos: the classical chaos growing is related to an exponential instability in correlators of particle trajectory coordinates, i.e. for example, the azimuthal angles ϕ ≡ ϕ(t): with c = C(t = 0), {... , ...} P.B. the Poisson brackets, P ϕ the conjugate momentum, t the time variable and t the initial fixed time. The constant λ introduced above is the Lyapunov exponent, which controls the chaos instability of circular photon trajectories. Eq.1 can be generalized to a quantum version for generic correlators as followŝ whereŴ ,V are two generic Hermitian correlators in the Heisenberg representation, ... = Z −1 Tr[e −βĤ ] is the thermal expectation value at temperature T = β −1 . Z and H are the partition function and the Hamiltonian of the N-body system respectively. Eq.(1) indicates that small deviations from the initial conditions related to trajectories diverge as e λt (conventionally t = 0) with a characteristic time constant τ ∼ λ −1 . The C in Eq.(1) andĈ in Eq.(2) start to be close to one after a critical "scrambling time" t * ∼ λ −1 log c −1 . According to MSS, any correlator in N body systems has a Lyapunov exponent that is below a critical thermal Lyapunov exponent: In the case of a BH, the maximal Lyapunov exponent corresponds to a λ BH = 2πk B T BH / where T BH is the Bekeinstein-Hawking temperature and , k B are the Planck and Boltzmann constants, respectively. Eqs. (1) indicates that the the chaotic instabilties are exponentially sensitive to the Lyapunov exponent, in turn related to the BH radius linear variations: This leads to the relative relations among the variations of the Lyapunov exponent, the BH temperature and the BH radius as follows: while the correlators relative to the linear variations are related to the temperature as follows: In a standard BH, in (semi)classical GR, if infalling matter/energy as well as Hawking radiation are negligible, the Lyapunov BH coefficient is constant in time, i.e. the expectation value is δλ(t) = 0. However, if the Lyapunov exponent has a time dependence, the linear correlator variation has a non-linear dependence inside the exponential function of time. Such a phenomenon can be extremely interesting in Extended Theories of Gravity, where there is a class of spherically symmetric solutions which are unstable under metric perturbations [12, 15-21, 23, 24, 26]. Such a dynamics is a consequence of the energy conditions where the presence of new degrees of freedom, due to Extended Gravity, play an important role into dynamics [21,22]. This fact gives rise to a variation in time of the BH radius, in turn altering the thermal Lyapunov BH exponent in time: where λ g is a new Lyapunov exponent related to the metric instabilities 3 . Thus, the maximally critical correlator has a double exponential form as and O(t) terms in the exponential of Eq.(8) are comparable when the time is t λ −1 g . In the case λ g << λ BH , the double exponentialization is completely negligible compared to the leading MSS exponential and Eq.(1) is re-obtained. However, if λ g λ BH the metric instability cannot be neglected anymore in the chaos correlator after a scrambling time.
In the next sections, we will discuss the instabilties in f (R) gravity and their implications in chaos around the BH shadows and QNMs.

CHAOTIC INSTABILITIES IN f (R) GRAVITY
In this section, we will describe the metric instability phenomena in f (R) gravity (see Refs. [12,[42][43][44] for detailed reviews on the subject).
f (R) gravity is a straightforward extension of the Einstein-Hilbert action where also non-linear terms in the Ricci scalar can be included. It is: where R is the Ricci scalar curvature, g ≡ detg µν is the determinant of the metric tensor. Specifically, f (R) is a generic function of R and S M is the matter action. We set the Newton constant to G N = 1. The field equations read as where T µν is the stress-energy tensor of matter and f R (R) ≡ df /dR. As it is well known, Eq.(10) admits spherically symmetric metric solutions [45] of the form where it is possible to restrict to the sub-class of solutions where M is the BH mass and L 2 dS = 3/Λ is the de Sitter curvature radius. Here Λ is the cosmological constant. It is worth noticing that a particular case of Eq.(11) is the Schwarzschild-de Sitter solution, corresponding to f, g reported in Eq. (12). As discussed in [45], any f (R) = R model in vacuum can be reduced to R + Λ. Thus f, g can be rewritten as where r 0 = −r + − r Λ . Here r + and r Λ are the BH and cosmological horizon respectively. The two radii are related with the dS-length and the BH mass as follows Let us consider a near-extremal limit where the BH radius is almost saturating the Hubble radius: In this limit, it is so that all quantities can be expressed in terms of r + . The null-geodesic circular orbit has a radius associated to a circular orbit angular velocity In the near-extremal limit, one obtains the Nariai BH metric solution, that can be written as follows: being where −∞ < t < ∞ corresponds to −π/2 < τ < π/2. Eqs. (10) in the background (20) read as: The perturbation equations are obtained substituting Eqs. (26) into Eqs. (22), (23), (24), (25) and they correspond to where As it is well known, this set of linear equations contain unstable modes as follows (see Ref. [24]): where t is the cosmological time (not τ ). Here m is the normalized mass in G N = 1 dimensions. The stability is related to the following parameter in turn related to Eqs. (32) and the following linear equation: which simplifies in the regime tanh t 1 for t >> 0. In this situation, 0 < α < 1/2 or α < 0 are the stability regions, while there are two new branches related to the Lyapunov exponents: In the case of solution λ g+ , a critical branch is comparable to the Lyapunov BH exponent. It corresponds to compare the linear perturbation effective mass and the BH temperature, that is Here m 2 and T are both adimensional being G N = M P l = c = = k B = 1. Thus, for macroscopic BHs, the effective (adimensional) mass for the scalar metric perturbation must be extremely suppressed in order to not have a double exponentiation in chaos correlators. In other words, if m 2 > 2πT , one would expect that the MSS bound is violated: there is a dominant Lyapunov exponent not related to the Hawking temperature and corresponding to the linear metric instability. In Fig.2, some BH Lyapunov branches are shown in the cases of evaporating and antievaporating solutions. An important issue is if realistic classes of f (R)gravity models exist exhibiting such a phenomenon or if these aspects are confined within unphysical toy models. Interestingly, the metric instabilities are present in many f (R) gravity models. In this perspective Let us consider three f (R) cases interesting for cosmology and relativistic astrophysics.

The polynomial model
This class of models have been largely investigated starting from the early f (R) = R + γ 2 R 2 . It can be easily generalized as where γ n is a constant with dimensionality dependent on the power n. In this case, the de Sitter solution is obtained for In the dS-branch, Eq.(33) corresponds to It is worth noticing that the stability conditions from Eq.(34) are violated for any n > 2, γ n > 0. Interestingly, n = 2 lives on the edge of the stability condition and not any metric instabilties are predicted, being α = 1/2 and n 2 = 0.

The exponential model
Another interesting case is given by exponential model, see Refs. [47][48][49]. The Lagrangian is Such a model can easily converges to the ΛCDM model while all Solar System tests are recovered in the post-Newtonian limit. In this case, a de-Sitter solution is obtained for R 0 = 4Λ = 3.75Λ ef f , where Λ 0.95Λ ef f . The corresponding α = 0.09 is compatible with the stability bounds; i.e. in the case of the exponential model the metric instabilities are turned off.

The Hu-Sawicki model
This model is particularly interesting to explain the accelerated behavior of cosmic fluid in the framework of f (R) gravity, see [46]. The Lagrangian is A de Sitter solution is achieved for R 0 = 4Λ, where for R 0 = 3.95Λ and Λ = 0.99Λ. The corresponding α = 0.02 lies into the metric stability condition.
To summarize, the metric instabilties propagate in polynomial f (R) models while they are suppressed in exponential and Hu-Sawicki models.

CHAOS AND THE MALDACENA-SHENKER-STANFORD CONJECTURE
In the previous section, we showed the presence of a new Lyapunov exponent, related to the f (R) action, which can be larger than the MSS bound. Here, we wish to explore and clarify how the MSS violation can be possible in f (R) gravity by analyzing generic Out of Thermal Equilibrium (OTO) correlators.
As mentioned in Introduction, the MSS conjecture states that, considering a generic N-particle system, the influence of chaos in commutators of two (Hermitian) operators can evolve in time no faster than exponentially and with a Lyapunov exponent λ ≤ 2πT (k B = = 1), where T is the temperature of the system. In the case of BHs, the temperature coincides with the Hawking temperature. Eq.(2) is an example of operator which provides a measure of chaos propagation in the N-field system. In particular, if such a correlator can be effectively factorized in form C(t) 2 V (0)V (0) W (t)W (t) for large t, then this will be a clear manifestation of the butterfly effect 4 .
In the time behavior of C, there are two important characteristic time scales: the dissipation time or collision time t d , which is typically t d ∼ β = 1/T in strongly coupled systems, which controls the exponential decay of the two point function V (0)V (0) and for C(t) ∼ 2 e λt , it corresponds to t * ∼ λ −1 log −1 . Typically, in the macroscopic limit, t d << t * , their ratio is proportional to the Planck constant (→ 0 in the classical limit).
While the C-operator can be safely computed in condensed matter lattice systems, it has regularization problems in Quantum Field Theory as reported in [3]. Alternatively, a more controllable prescription consists in treating the following correlator:C Such a correlator is related to which is the one investigated in the main MSS argument [3]. Indeed, the relation between Eq.(42) and Eq. (43) is with V ≡ V (0). The first two terms can be re-absorbed in thermal state normalizations, while only the last F -terms contribute to the time growth [3]. As mentioned above, for t d << t << t * , the F (t) is expected to be approximately factorable to as a product of disconnected correlators of V and W . In the large N-field limit, the factorization has the form where the first two terms exponentially decay for t >> t d . In a chaotic system, the F (t) departs from the constant value F d after a critical scrambling time t * . MSS conjectured that such a growth rate is bounded as follows: where is an initial small quantity, and The conjecture is supported by a simple argument that we review in Appendix A. The main assumption of the MSS argument consists into the treatment of the error variable as where the error ≡ t depends on the reference time considered. In other words, Eq.(49) converges to Eq.(47) only under proper assumptions for the time t that imply t << 1. In case GR, it is sufficient to consider a time scale much below the BH information scrambling time for ensuring that Eq.(49) converges to Eq.(47), i.e. t<<t * << 1. Such an apparently innocuous assumption can be violated in f (R) gravity where metric instabilities with a new Lyapunov exponent λ g ≥ λ are present. Indeed, the correlators can violate the Eq.(49) bound for t < t * , if the characteristic time scale related to the metric instability is below the thermal MSS Lyapunov exponent. In this case, a new characteristic scrambling time, related to the metric instabilities, appears as t g * ∼ λ −1 g log(λ 0 BH ) −1 . In particular, for BH growing solutions, Eq.(49) is replaced by leading to a double exponentialization of the F -function. Behind the MSS violation coming from f (R) gravity metric instabilities, there is an important difference between f (R) and GR related to their BH thermodynamical properties. Indeed, an important assumption of the MSS conjecture is that the N-system has a fixed temperature T and thus it is in thermal equilibrium. In other words, the chaos effects of correlators are thought as a non-equilibrium transient before reaching an equilibrium to the final thermalization. Indeed the Hawking radiation is computed from semiclassical saddle approximations of the Euclidean quantum gravity action, coinciding with the partition function of the BH emitting Hawking radiation in thermal equilibrium with the external environment. On the other hand, if the BH metric has a horizon area growing without reaching any constant plateau in time, as an effect of growing metric instabilities of f (R) gravity, then a perfect thermal equilibrium will be never reached. As pointed out in Ref. [17], a dynamical BH horizon leads to revisit the standard BH thermodynamical considerations, including the BH temperature. The time exponential growing of the BH radius leads to a continuous increasing of chaos in the system without reaching any asymptotic thermalization.
In Ref. [3], MSS comment on the possible effect of higher-derivative corrections to Einstein gravity, arguing that the bound does not receive any corrections. Such a statement is based on a result holographically connected to Einstein gravity related to the AdS/CFT correspondence [51,52], that is: where f 0,1 ∼ O(1) are positive constants which depend on the particular V, W considered. Such a result saturates the MSS bound, suggesting a t d ∼ T −1 BH and t * = T −1 BH log N 2 . However, the apparent disagreement of our results with CFT computations in the case of f (R) metric instability can be understood from several prospectives. First of all, f (R) gravity cannot be fully identified with a GR with quantum higher-derivative corrections. Indeed, f (R) gravity can be rewritten as GR with a scalar field after a conformal transformation. In this sense, the MSS bound stability with respect to higher derivatives does not apply on f (R) gravity. On the other hand, the dS/CFT holographic equivalence is notoriously a not so solid conjecture as (specific) AdS/CFT cases. Indeed, the holographic principle works in physical systems where no any metric perturbations destabilized the fixed (A)dS background. In a dynamical background, the holographic principle may loose any controllability and reasonable sense.
The case of a BH horizon saturating a de Sitter Hubble radius, the extremal Schwarzschild-de Sitter solution, was never studied from the perspective of dS/CFT and, in principle, such a correspondence may be lost in f (R) gravity or just from quantum gravity higher-derivative operators beyond specific special and integrable cases such as 5 AdS 5 × S 5 /CF T 4 , AdS 2 /CF T 1 , ....

THE GEODESIC INSTABILITY
In classical chaotic systems, the Lyapunov parameters control the average rate of convergence or divergence of body trajectories within the phase space {X i , P i }, with i the space-coordinate index.
Let us consider the generic classical equation of motion where P j ≡ P j (X i ). Let us consider a certain trajectoryX i (t) which is a solution of the Eq. (52). In order to perform the geodesic stability analysis, we can consider the linearized equation of motion around the trajectoryX i , as follows: where K ij (t) = ∂P j /∂X i | Xj (t) . The integral solution of the linearized equation can be rewritten in the formal form where the operator U satisfies the following equation: The principal Lyapunov parameter λ is related to the following asymptotic limit: However, in our case, the particle orbit has a radius that is thought, in first approximation, as slowly growing in time: In the case of λ g > 0 and t << λ −1 g , such a dynamics would correspond to a spiral trajectory slowly increasing the radial coordinate in time. In the opposite regime, where t >> λ −1 g , the radial instabilities would dominate on any other circular instabilities.
In the case of circular geodetic orbits, the K matrix and the Lyapunov parameter can be related to the metric tensor components: where In our case, in the approximation of a slow radial time-variation, a pure circular trajectory will be substituted by a spiral trajectory with a slowly increasing radius in time. As a specific case of Eq.(57), we obtain an exponential growth in time of the Lyapunov exponent: where which can be reduced to a point-like Lagrangian and, from the Euler-Lagrange equations, we get the effective potential The Lyapunov exponent can be related to the effective potential as follows: In Eqs.(62)-(64) a slow time variation has been considered, neglecting all time derivative terms of the BH radius, i.e. r BH 0,r BH 0, .... With these onsiderations in mind, let us discuss now the quasi normal modes.

THE QUASI NORMAL MODES
Let us consider a probe scalar field around the BH: it is described by the Klein-Gordon (K.G.) equation in the curved background as where ∂ µ is the partial four-derivative. In a spherical space-time background, the K.G. equation has a solution that can be decomposed as where Y l ≡ Y l (θ, φ) are the spherical harmonic function with l the angular momentum eigenvalues satisfying the following inserting Eq.(66) in the K.G. equation, one obtains the following equation for the radial part of the scalar field: known as the Regge-Wheeler wave equation [50]. Let us introduce now the tortoise radial coordinate The V s (r) appearing in Eq.(68) is the effective potential reading V s (r) = f (r) l(l + 1) In the case of a generic integer spin field such a potential can be generalized as V s (r) = f (r) l(l + 1) However, in case of metric instabilities, the background has a radius changing in time; thus the potential V s (r) would have an effective time dependence through the time variation of the BH radius. Assuming that the time-variation of the BH radius is much slower than the characteristic field frequency, we can neglect anyṀ ,M , ... derivatives and we can treat the time dependence inside the potential redefining f (r) and M (r) as V s (r, t) f (r, r BH (t)) l(l + 1) where time derivative terms of r BH (t) are neglected. The usual approach for solving the Regge-Wheeler equation is a separation of the wave function assuming R(r, t) = Ψ(r)ζ(t), with ζ(t) ∼ e −iωt , obtaining where ω corresponds to the QNM frequency, having a real and an imaginary parts as with ω R the real oscillation frequency and ω I the imaginary part related to damping or exponential instabilities of the field modes. However, the separation of radial and time parts cannot be performed, in our case, if the background time variation is comparable with the inverse of the characteristic real and imaginary parts.
Assuming that m << ω R , ω I , the time derivative can be neglected and we obtain where ψ has a slow varying time variation and R(r, t) = ψ(r, t)e −iω(t)t , where ω(t) derivatives are neglected, i.e. d n ω/dt n = 0 for any n > 0. Under this assumption, we can search for QNM as saddle solutions of the following WKB condition: where the only time dependence in the potential enters through the r 0 (t). This implies that ω has a slow time dependence as ω ≡ ω(t). The V (r 0 (t)) is the second derivative of the effective potential with respect to r at the maximum point r 0 related to dV /dr * | r * =r0 = 0. In our case, the maximum is slowly changing in time. In the deep WKB regime, l >> 1 and q(r, t) and Such a condition corresponds to the following QNM frequency: which can be rewriten as ω QN M (t) = Ω c (r 0 (t))l − i(n + 1/2)|λ(t)| .
In the case of the near Schwarzschild-de Sitter BH, λ = Ω c where Ω c (r 0 (t)) = r Λ (t) − r + (t) 2r 2 + (t) In the case of a perfectly extremal Schwarzschild-de Sitter BH, the horizons coincide as soon as r Λ = r + and thus Ω c = 0; i.e. no QNM is found. From Eq.(84), one expects that Ω c decreases with the cosmological time if a + > 0 (antievaporation) while increases for other cases (evaporation). This corresponds to the case where the QNMs can be altered by the slow dynamical variation of the BH geometry sourced by gravitational actions like f (R) gravity.

THE ANOMALOUS BLACK HOLE SHADOW
As mentioned in the Introduction, the BH shadow is defined as the area delimited by the critical photon-trajectory 6 , r C = 3r BH /2 = 3M BH [56]. Nevertheless, the critical impact parameter, where light is trapped, is higher than the r C by a factor √ 3, i.e. b C √ 3r C 5.2 r BH [57]. The photon-sphere appears as a luminous ring around the BH shadow, notoriously interpreted as the evidence of the BH horizon. Indeed, EHT observed the M 87 * center with a BH size resolution, measuring the photon-ring around the BH shadow. The interpretation favors the presence of a BH horizon and it seems to rule out horizonless exotic compact objects [36][37][38] 7 .
Here, we wish to remark that a test for the BH horizon existence does not coincide with a general test for Extended Theories of Gravity, including f (R) gravity. Indeed, there are several BH solutions, having an external null-like event horizon, obtained in large classes of alternative theories of gravity and, in particular, in f (R) gravity. In this category can be included the Schwarzschild-de Sitter solutions discussed above. Nevertheless, the metric instabilities of f (R) gravity, related to the new Lyapunov exponent λ g , can give rise to anomalous variations of the BH radius and, thus, of the the BH shadow with respect to similar features in GR.
Specifically, the ratios among the critical photon-orbit, the impact parameter and the photon-ring are fixed and controlled by the time-varying BH radius: where R ± are the maximal and minimal radii delimiting the photon-ring and ∆R is the photon-ring thick length; k is a parameter giving the numerical hierarchy between b C and R − that will be discussed soon. Eq.(85) implies that the photon-ring is time-dependent and it is controlled by the metric Lyapunov exponent λ g , while ∆R remains constant in time.
For fixing the ideas, let us consider light emission with an intensity peak I 0 around the critical photon-orbit, as displayed in the left side of Fig.3. It is worth to remark that the emitted intensity profile significantly differs from the observed intensity on terrestrial experiments. Indeed, in the right side of Fig.3, we show how the profile of emitted intensity would appear for an asymptotic observer. In this case, the observed intensity peak appears as a narrow distribution around r/M 5.5, while it is highly suppressed for r/M < 4 and r/M > 6 [58][59][60]. The lensing ring emission is within a thin ring because of the BH demagnification effect at r/M > 6 and the efficient gravitational captures within the BH shadow from r/M < b C .
Let us now consider the case of a positive geometric Lyapunov exponent: in this case, while the emitted intensity spectrum retains the same profile in time, the photon-ring dynamically grows in time preserving its geometrical profile (see Fig.3). The photon-ring, at a fixed time t 0 , is absorbed by the growing BH shadow at t 1 > t 0 , but it is replaced by higher distance photons within R ± (t 1 ) > R ± (t 0 ).
Such a phenomenon could be tested from BH shadow measurements, searching for anomalous growing or shrinking of the photon-ring in time. Clearly, due to the very fine measurements required for this pupose, it could be feasible for the Next Generation Event Horizon Telescope Design Program [66].

r/M(t)
Ie/I0 rc(t)/M(t) Figure 2: An example of emitted intensity profile Ie, normalized with the maximal peak I0, as a function of the radius, normalized by the time-dependent BH mass M (t), is displayed on the left. On the right side, we show its corresponding photon-ring intensity seen by an asymptotically far observer. In case of λg > 0, the photon-ring dynamically increases with time (here the case of t > λ −1 g is qualitatively displayed).

PERSPECTIVES AND CONCLUSIONS
In this paper, we analyzed the emergence of chaos for circular trajectories around spherically symmetric BHs in the framework of f (R) gravity. In particular, we focalized on the spherically symmetric solutions exhibiting an anomalous growing or shrinking of the BH event horizon area sourced by metric instabilities. We found that the The new Lyapunov exponent of f (R) gravity dominates on the MSS one if the effective scalaron mass is heavier than a critical scale proportional to the BH temperature. The BH temperature for supermassive BHs of the order of 10 9 M is around 10 −17 Kelvin corresponding to 10 −21 eV or so. Thus if the scalaron has a mass higher than 10 −21 eV , the MSS can be violated for a Supermassive BH in the context of f (R) gravity. We also showed that such metric instabilities lead to a distortion of the Quasi-Normal-Modes, with a potential impact on BH ringdown phase after the BH merging emitting gravitational waves. In conclusion, these features could be, from one hand, signatures to test alternative theories of gravity. On the other hand, they could constitute straightforward explanations for classifying anomalies of compact objects, if revealed by next generation experiments.
Let us define the function where is an infinitesimal quantity and t 0 a reference time. At this point the proof of the MSS bound will be obtained if the f (t) satisfies the hypothesis of the following theorem.
Theorem. Let us consider a function f (t) satisfying the following properties: i) f (t + iτ ) is an analytic function in t > 0 and −β/4 ≤ τ ≤ β/4, where t, τ are the real and the imaginary parts of the complex variable z = t + iτ . Let us assume that f (t) is real for τ = 0.
ii) |f (t + iτ )| ≤ 1 in t > 0 and −β/4 ≤ τ ≤ β/4. Then, the (i)-(ii) hypothesis imply: The Theorem can be proved considering a map of the half strip {t > 0, −β/4 ≤ τ ≤ β/4} to the unit circle, in the complex plane. Such a transformation can be obtained using the following map: The f (Z) is an analytic function in the unit circle as guaranteed by the (ii)-hypothesis. From the Schwarz-Pick theorem, the function does not increase in the hyperbolic metric ds 2 = 4dZZ/(1 − ZZ) 2 , i.e.
Such an inequality corresponds to Eq.(88) for τ = 0 (having use the (i) hypothesis): The rest of the MSS argument reduces to show that Eq.(87) satisfies |f (t + iτ )| ≤ 1 on three boundaries of the half strip {t > 0, −β/4 ≤ τ < β} and that f is constant in every points inside the half strip. Indeed, if the previous properties are proved, then the Phragmén-Lindelöf principle will guarantee that |f | will be bounded by 1 in every points of the half-strip interior. On the |τ | = β/4 edges, the F reduces to Such an equation can be considered as a contraction of two matrix elements [ζV W ζ] ij and [ζW V ζ] ij ; thus, we can apply the Cauchy-Schwarz inequality For times larger than the dissipation timescale, the operators factorize with a certain error where the error depends on the reference time t and the first term on the right side is nothing but the definition of F d . If we consider t much minor than the scrambling time, then the will be small and the the |f | ≤ 1 condition will be guaranteed. Let us now consider the boundary t = 0, corresponding to F (iτ ) with −π/4 ≤ τ ≤ π/4. Following similar arguments as above, assuming t sufficiently smaller than the scrambling time, it is easy to show that |f | ≤ 1 is satisfied.
Finally, to complete the proof, we show that |f (z)| ≤ C with C as a constant with C ≥ 0. Using the chaos factorization and Hermiticity of the V and W operators, we obtain where α = 4τ /β. Thus, the theorem can be applied on Eq.(87) and we obtain where << 1 if t << t * , with t * the scrambling time. As we discussed in Sec.3, this last assumption should be relaxed in case of f (R) gravity with metric instabilities, evading the MSS argument.