Observational constraints on Yukawa cosmology and connection with black hole shadows

We confront Yukawa modified cosmology, proposed in arXiv:2304.11492 [Jusufi et al. arXiv:2304.11492], with data from Supernovae Type Ia (SNe Ia) and Hubble parameter (OHD) observations. Yukawa cosmology is obtained from a Yukawa-like gravitational potential, with coupling parameter $\alpha$ and wavelength parameter $\lambda$, which gives rise to modified Friedmann equations. We show that the agreement with observations is very efficient, and within $1\sigma$ confidence level we find the best-fit parameters $\lambda=\left(2693_{-1262}^{+1191}\right)\, \rm Mpc$, $\alpha=0.416_{-0.326}^{+1.137}$, and a graviton mass of $m_{g}=\left(2.374_{-0.728}^{+2.095}\right)\times 10^{-42}\, \text{GeV}$. Additionally, we establish a connection between the effective dark matter and dark energy density parameters and the angular radius of the black hole shadow of the SgrA and M87 black holes in the low-redshift limit, consistent with the Event Horizon Telescope findings.


I. INTRODUCTION
According to modern cosmology, the universe's largescale structure is homogeneous and isotropic. Additionally, it is believed that cold dark matter, a type of matter that is not visible and only interacts through gravity, exists [1][2][3][4][5]. However, despite numerous efforts, there have not been any direct detection of dark matter particles and their existence is only inferred from its gravitational effects on galaxies and larger structures. On the other hand, dark energy is also introduced to explain the universe accelerated expansion [6], supported by numerous observations [7][8][9].
On the other hand, black holes are intriguing astronomical objects and can potentially test the theories of gravity in strong gravitational fields. One of the most fascinating aspects of black holes is their shadow image. The black hole's silhouette is a dark region that results from the immense gravitational pull of the black hole, which bends the path of light rays near it. Specifically, photons emitted from a bright source close to a black hole can either be drawn into the black hole or scattered away from it and into infinity. Additionally, critical geodesics separate the first two sets, known as unstable spherical orbits. By observing the critical photon geodesic trajectories in the sky, we can obtain the black hole shadow . With the recent results, the black hole shadow for the supermassive black holes M87 and Sgr A was confirmed by the Event Horizon Telescope (EHT) collaboration [77][78][79][80][81][82].
In the present paper, we follow an approach motivated by cosmology and quantum field theories; we aim to study the dark sector by introducing the Yukawa potential [83][84][85][86][87]. We adopt the viewpoint of Verlinde and con-sider gravity as an entropic force caused by the changes in the system's information [88]. Verlinde further argued that dark matter is an apparent effect, i.e., a consequence of the baryonic matter [89]. Furthermore, the corresponding entropic force was used in deriving the corrected Friedmann equations due to the minimal length, as recently studied in [90][91][92][93].
In particular, as it was recently shown in [93], dark matter can be explained by the coupling between baryonic matter through a long-range force via the Yukawa gravitational potential. This coupling is characterized by the coupling parameter α, the wavelength parameter λ, and the Planck length l 0 . The modified Friedmann equations are derived using Verlinde's entropic force interpretation of gravity based on the holographic scenario and the equipartition law of energy. An equation connects the dark matter density, dark energy density, and baryonic matter density. It is worth noting that dark matter is not associated with a particle but is an apparent effect. Dark energy is related to graviton mass and α, indicating that the cosmological constant can be viewed as a self-interaction effect between gravitons. The model parameters were estimated as λ ≃ 10 3 [Mpc] and α ∈ (0.0385, 0.0450).
In this work, we are interested in performing detailed observational tests of the cosmological scenario based on Yukawa potential. In particular, we wish to constrain the parameters of the model using Supernovae Type Ia (SNe Ia) and Hubble parameter (OHD) observations. Additionally, we are interested in investigating the connection to black hole physics. In particular, at low redshifts, one can use the angular radius of the black hole shadow to constrain the Hubble constant independently. The manuscript is organized as follows. In Sect. II we review the Yukawa cosmological scenario, while in Sect. III we extract observational constraints. In Sect. IV we study the relation between the modified Yukawa cosmology and black hole shadows, and in Sect. V we comment on our findings.

II. YUKAWA MODIFIED COSMOLOGY
In this section, we shall review the model that was recently studied in [93]. The gravitational potential considered is modified via the non-singular Yukawa-type gravitational potential with l 0 being a small quantity of Planck length order, i.e. l 0 ∼ 10 −34 cm and α > 0. Note that the wavelength of massive graviton we have λ = ℏ mgc > 10 20 m, that leads to m g < 10 −64 kg for the graviton mass [94]. If we use the relation F = −∇Φ(r)| r=R , and by neglecting the term αl 2 0 /R 2 → 0 we get the modified Newton's law of (2) We can further elaborate that such a force can be obtained if we modify the expression for the entropy by means of Verlinde's entropic force interpretation. In this theory, when a test particle or excitation moves apart from the holographic screen, the magnitude of the entropic force on this body has the form [88] F △x = T △S. Specifically, one can get Newton's law of gravitation via the entropy-area relationship S = A/4, however in general, one can modify the total expression for the entropy as S = A/4 + S(A) [93]. One uses the holographic scenario and the equipartition law of energy to get the expression for force. The entropy of the surface changes by dS = [1/4 + ∂S/∂A] dA; at the same time, we can relate the area of the surface to the number of bytes according to A = QN, where Q is a fundamental constant, and N is the number of bytes. We can now use the equipartition law of energy, we get the total energy on the surface via E = N k B T /2, and further taking △N = 1, and △A = Q, we get the modified gravity force This means that, by changing the entropy, we end up with a modified law of gravity. Let us define η = 1/8πk B ; we get Q = 1, Since for large-scale distances l 0 is unimportant, we can set it to zero, i.e. l 0 = 0. We now see that we can get the corrections to the entropy if we compare 1 + 1 2πR Solving for entropy, we obtain [93] We have therefore shown that the Yukawa modified force follows from the modification of the entropy. This method can be viewed as a generic result; if we modify the entropy, we modify Newton's law of gravity. It is worth noting that the correction of the entropy can be viewed as a volume law entanglement to the entropy due to the gravitons. We see from the last equation that α appears only in the second term; hence we can say that α is a result of the entanglement due to the volume law entropy contribution. We proceed by studying the implications of the above-modified law of gravity in cosmology. First, we assume the background spacetime to be spatially homogeneous and isotropic, described by the Friedmann-Robertson-Walker (FRW) metric ds 2 = −dt 2 + a 2 (t) dr 2 1 − kr 2 + r 2 (dθ 2 + sin 2 θdϕ 2 ) , (7) with R = a(t)r, x 0 = t, x 1 = r, and the two dimensional metric h µν , and where k is the spatial curvature (k = 0, 1, −1 corresponding to flat, closed, and open universes, respectively). In addition, we have an apparent dynamic horizon, which the following relation can determine h µν (∂ µ R) (∂ ν R) = 0. It is easy to show that the apparent horizon radius for the FRW universe reads as On the other hand, we have a matter source which can be assumed to be a perfect fluid described by the stressenergy tensor along with the continuity equationρ + 3H(ρ + p) = 0, with H =ȧ/a being the Hubble parameter. Let us consider a compact spatial region V with a compact boundary S, corresponding to a sphere of radius R = a(t)r, where r is a dimensionless quantity. Through Newton's law, we can write the gravitational force on a test particle m near the surface [93] as In Newtonian cosmology, we can take ρ = M/V inside the volume V = 4 3 πa 3 r 3 , hence we can rewritten the above equation as [93] which is the dynamical equation for Newtonian cosmology. To obtain the Friedmann equations in general relativity, we must use the active gravitational mass M rather than the total mass M . By replacing M with M, we obtain where the active gravitational mass can also be computed via Using these equations, we obtain the modified acceleration equation for the dynamical evolution of the FRW universe [93] Furthermore, we can simplify the work since l 0 is a very small number; we can consider a series expansion around x = 1/λ via (15) provided that αR 2 /λ 2 ≪ 1. In general, we expect α < 1, and λ to be some large number of magnitude comparable to the radius of the observable Universe R ∼ 10 26 m. In summary, the corresponding Friedmann equation where we have included several matter fluids with a constant equation of state parameters ω i along with the continuity equationρ i + 3H(1 + ω i )ρ i = 0, that yield an expression for densities ρ i = ρ i0 a −3(1+ωi) . Inserting these into (16) and integrating we obtain [93] Using the fact that R[a] = ra, we further, geṫ with r nearly a constant. Considering the equations of state, ω i / ∈ {−1, 1/3}, we havė implying that at leading order terms (l 2 0 /λ 2 → 0), where the Newton's constant is shifted as G eff = G(1+α), along with the definitions [93] If, for example, we assume only a matter source, at leading order terms we can write Focusing on the flat case (k = 0), we have R 2 = 1/H 2 , yielding Finally, by expanding around l 0 , making use of , and neglecting the terms ∼ O(l 0 α 2 /λ 2 ), we obtain

A. Late time universe
Let us now study the modified Friedmann equation's phenomenological aspects. In particular, we are interested in studying the late universe, which implies we can neglect the quantum effects by setting l 0 → 0 [Γ 1 = 0]. This gives and using ρ crit = 3 8πG H 2 0 we acquire two solutions: where Ω i = Ω i0 (1 + z) 3(1+ωi) , Ω i0 = 8πGρ i0 /(3H 2 0 ), with E = H/H 0 . In addition, we point out that the total quantity Ω 2 in the square root should be viewed as the root-mean-square density energy, i.e Ω ≡ ⟨Ω 2 ⟩ along with Ω 2 = Ω 2 B + Ω 2 Λ . As explained in [93], the most interesting implication of the last equation relies on the physical interpretation of the term 2Γ 2 (ω i )Ω i (1 + α)/H 2 0 . In particular, it was shown that this term precisely mimics the effect of cold dark matter of the ΛCDM model. Taking the term Γ 2 Ω i and set ω i = 0 we define the quantity (here we shall add the constant term c to make the equation consistent) [93] Thus, we can obtain an equation for dark matter as From this equation, we can deduce that dark matter can be viewed as an effective sector, a manifestation of modified Newton's law, quantified by ∼ α Ω B . Additionally, we define the quantity Finally, comparing the last expression with ρ Λ = Λc 2 8π G , we can estimate the effective cosmological constant to be Note that one can combine the above expressions and relate baryonic matter with the effective dark matter and dark energy, acquiring In summary, Eq. (27) can be re-written as [93] We can now introduce the split where

III. OBSERVATIONAL CONSTRAINTS
In the previous section, we presented Yukawa-modified cosmology Hence, in this section, we can proceed to observational confrontation with Hubble parameter data (OHD) and Type Ia supernovae (SN Ia) data, in order to extract constraints of the free parameters. For this purpose, we compute the best-fit of the free parameters and their corresponding confidence regions at 1σ (68.3%) of confidence level (CL) with the affine-invariant Markov chain Monte Carlo (MCMC) method [95], implemented in the pure-Python code emcee [96]. In particular, we have considered 100 chains or "walkers", using the autocorrelation time provided by the emcee module as a convergence test. In this sense, we computed at every 50 step the autocorrelation time, τ corr , of the chains. If the current step is larger than 50τ corr and the values of τ corr changed by less than 1%, then we will consider that the chains are converged, and the constraint is stopped. We discard the first 5τ corr steps as "burn-in" steps. Finally, we compute the mean acceptance fraction of the chains, which must have a value between 0.2 and 0.5 [96] and can be modified by the stretch move provided by the emcee module.
For this Bayesian statistical analysis, we need to construct the following Gaussian likelihoods: where χ 2 OHD and χ 2 SNe are the merit function of the OHD and SNe Ia data, respectively. The Gaussian likelihood for the joint analysis SNe Ia+OHD is constructed as In the following subsections, we will briefly describe the construction of the merit function of each data set.

A. Observational Hubble parameter data
For the OHD, we consider the sample compiled by Magaña et al. [97], which consists of 51 data points in the redshift range 0.07 ≤ z ≤ 2.36, and for which we construct their merit functions as where H i is the observational Hubble parameter at redshift z i with an associated error σ H,i , all of them provided by the OHD sample, H th is the theoretical Hubble parameter at the same redshift, and θ encompasses the free parameters of the model under study. The theoretical Hubble parameter is obtained from Eq. (33), which we conveniently rewrite as with Ω Λ,0 given by Eq. (30), while the condition H(z = 0) = H 0 , leads to and the free parameters of the Yukawa-modified cosmology are θ = {H 0 ; Ω B,0 ; Ω Λ,0 }. Note that one has the relations and Eq. (37) becomes which is related to the Hubble parameter for the stan- It is important to mention that the OHD, as well as the SNe Ia data, are not able to independently constraint Ω ΛCDM ; Ω m,0 }. Therefore, we consider the following priors for our MCMC analysis: H 0 = 100 km/s M pc h, with 0.55 < h < 0.85, 0 < Ω m,0 < 1, 0 < Ω B,0 < 0.2, 0 < Ω Λ,0 < 1, and the condition α > 0 implies 0 < Ω B,0 + 2Ω B,0 Ω Λ,0 + Ω Λ,0 < 1.

B. Type Ia supernovae data
For the SNe Ia data, we consider the Pantheon+ sample [98], which is the successor of the original Pantheon sample [99] and consist of 1701 data points in the redshift range 0.001 ≤ z ≤ 2.26. In this case, the merit function can be conveniently constructed in matrix notation (denoted by bold symbols) as where [∆D(z, θ, M )] i = m B,i − M − µ th (z i , θ) and C = C stat + C sys is the total uncertainty covariance matrix, where the matrices C stat and C sys accounts for the statistical and systematic uncertainties, respectively. In this expression, µ i = m B,i − M is the observational distance modulus of Pantheon+ sample, obtained by a modified version of the Trip's formula [100], with three nuisance parameters calibrated to zero with the BBC (BEAMS with Bias Corrections) approach [101]. Therefore, the Pantheon+ sample provides directly the corrected apparent B-band magnitude m B,i of a fiducial SNe Ia at redshift z i , with M the fiducial magnitude of an SNe Ia, which must be jointly estimated with the free parameters of the model under study. The theoretical distance modulus for a spatially flat FLRW spacetime is given by with d L (z i , θ) the luminosity distance given by where c is the speed of light given in units of km/s. Note that the luminosity distance depends on the theoretical Hubble parameter, which is given by Eq. (37) for the Yukawa cosmology and Eq. (42) for the ΛCDM model. Therefore, we only add to the free parameters the nuisance parameter M , for which we consider the following prior to our MCMC analysis: −20 < M < −18.
Similarly to the Pantheon sample, there is a degeneration between the nuisance parameter M and H 0 . Hence, to constraint the free parameter H 0 using SNe Ia data alone, it is necessary to include the SH0ES (Supernovae and H 0 for the Equation of State of the dark energy program) Cepheid host distance anchors of the form the Cepheid calibrated host-galaxy distance obtained by SH0ES [102]. Hence, for the correspondence, we use the Cepheid distances as the "theory" instead of using the model under study to calibrate M , considering that the difference µ i (M ) − µ Cepheid i is sensitive to M and H 0 and also is largely insensitive to other parameters like Ω m,0 . In this sense, the Pantheon+ sample provides µ Cepheid i , and the total uncertainty covariance matrix for Cepheid is contained in the total uncertainty covariance matrix C. Therefore, we can define the merit function for the SNe Ia data as where (48) It is essential to mention that from now on, we will omit the free parameter M and we will focus our analysis only on the free parameters of each model. Besides, considering that the best-fit parameters minimize the merit function, we can use the evaluation of the best-fit parameters in the merit function, χ 2 min , as an indicator of the goodness of the fit: the smaller the value of χ 2 min is, the better is the fit.

C. Results and discussions
In Table I, we present the total number of steps, the values used for the stretch move, the mean acceptance fraction, and the autocorrelation time for the free parameters h and Ω m,0 of the ΛCDM model, and h, Ω B,0 , and Ω Λ,0 of the Yukawa modified cosmology. Additionally, in Table II, we present their respective best-fit values at 1σ CL with the corresponding χ 2 min criteria. In Figs. 1 and 2, we depict the posterior 1D distribution and the joint marginalized regions of the free parameters space of the ΛCDM model and the Yukawa-modified cosmology. The admissible joint regions presented correspond to 1σ, 2σ (95.5%), and 3σ (99.7%) CL, respectively. These results were obtained by the MCMC analysis described in Section III for the SNe Ia data, OHD, and their joint analysis.
As we can see from Table II, the values obtained for the χ 2 min criteria show that the Yukawa modified cosmology can fit the observational data of SNe Ia, OHD and SNe Ia+OHD as accurately as the ΛCDM model. Even more, the value of the Hubble constant is the same in both models, which agrees with our previous identification in Eqs. (37) and (42), where H 0 = H ΛCDM 0 . The only difference between these models relies on the rescaling of energy densities due to the contribution of the α parameter. On physical grounds, the main difference between these models are that in Yukawa cosmology, dark matter is effective and precisely mimics the cold dark matter of the ΛCDM scenario.
To establish the last point, we use the results of our MCMC analysis for SNe Ia+OHD to calculate the values of Ω ΛCDM Posterior 1D distribution and joint marginalized regions of the free parameters space of the ΛCDM model, obtained by the MCMC analysis described in Section III. The admissible joint regions correspond to 1σ, 2σ, and 3σ CL, respectively. The best-fit values for each model free parameter are shown in Table II.

FIG. 2:
Posterior 1D distribution and joint marginalized regions of the free parameters space of the Yukawa modified cosmology, obtained by the MCMC analysis described in Section III. The admissible joint regions correspond to 1σ, 2σ, and 3σ CL, respectively. The best-fit values for each model free parameter are shown in Table II    ogy can mimic the late-time ΛCDM model, as we can see from Fig. 3, where we depict the Hubble parameter as a function of redshift z for the Yukawa and ΛCDM cosmologies, given respectively by Eqs. (37) and (42). Furthermore, in Fig. 4 we depict the matter density pa-rameter Ω m as a function of redshift for both models, were Ω m = Ω m,0 (1 + z) 3 /E 2 . Both figures were obtained with the results of our MCMC analysis described in Section III for the joint analysis.
The joint analysis has revealed the best-fits values for Yukawa parameters at 1σ CL. These are λ We have compared our best-fit values for α, λ, and m g at 1σ CL with previous references in Table III. Our results are consistent with more stringent Yukawa constraints on m g [105]. For instance, our findings show that the graviton mass is less than 6 × 10 −41 GeV and λ is greater than 97.22 Mpc obtained for Weak Lensing measurements [105], and α is less than the value of 0.581 obtained for the Local Group estimations [107]. Also, we provide Diagram 5 for upper bounds on the graviton mass (m g ) for the Yukawa potential, which we compared with our best-fit values using SNe Ia+OHD. Our research has revealed that the graviton's mass is lower than previously estimated. This discovery significantly contributes to the scientific community and complements other significant findings related to Yukawa constraints from various systems [103][104][105][106][107][108][109][110][111].

IV. RELATING YUKAWA COSMOLOGY AND BLACK HOLE SHADOWS
In this section, we will present a connection between the dark matter/dark energy densities and the angular radius of the black hole shadow. We will closely follow the approach developed in [74,75], in which one can employ the standard definition of the luminosity distance for a flat ΛCDM model as [74,75] where c is the speed of light, H 0 is the Hubble constant. The quantity I(z) is given in terms of the integral (1 +z) 3 + Ω ΛCDM with Ω ΛCDM m,0 = Ω ΛCDM B,0 + Ω ΛCDM DM,0 , that provides the present values of the critical density parameters for matter and a dark energy component, respectively. On the other hand, one can define the luminosity distance, which is related to the angular diameter distance d A (z) in terms of the equation [74,75]: To examine the connection with black holes shadows, let us recall that, by definition, the observed angular diameter of an object (say, a black hole) is given by θ = R/d A , with R is some proper diameter of the object. Hence, we can extract information about the cosmological parameters, having measured one of these distances at a particular redshift z. Assuming Schwarzschild black holes, as a first approximation, an observer located far away from the black hole can construct the shadow image in the center with an angular radiuŝ where R SH is the shadow radius and M BH is the mass of the supermassive black hole. As it was pointed out in [75], the above equation for the angular radius of the black hole shadow is valid only when the radial coordinate is large enough in comparison with the size of the black hole shadow radius 3 √ 3GM BH /c 2 . If we now combine Eqs. (49)-(51), we can obtain Tsupko et al. [74], Escamilla-Rivera and Torres Castillejos [75] In the present work, we are interested in the lowredshift limit, hence utilizing Eqs. (50) and (53) we can obtain [74,75] implying that we can estimate the angular radius of the black hole if we know the Hubble constant, the redshift z of the black hole, along with its mass. Since we have explicit expressions of the dark matter/dark energy parameters in terms of H 0 (see Eqs. (29) and (30)), we can solve for H 0 in Eq. (30) and directly relate the angular radius with the dark energy density parameter aŝ or using the notation of (42): Similarly, we can express this relation in terms of the effective dark matter mass aŝ  or in terms of the notation of (39), (40) and (41), aŝ In summary, we can obtain the angular radius of black holes using dark matter and dark energy densities. Note that in this modified cosmological scenario, the shadow radius depends on the distance from the black hole to the observer.

A. Black hole solution
Let us see how the spacetime geometry around the black hole us modified in this theory. The general solution in case of a static, spherically symmetric source reads ds 2 = −f (r)dt 2 + dr 2 f (r) + r 2 (dθ 2 + sin 2 θdϕ 2 ). (59) The energy density of the modified matter can be computed from ρ(r) = 1 4π ∆Φ(r). In astrophysical scales we can set l 0 /r → 0, in that case using (1) we acquire The negative sign reflects that the energy conditions are violated inside the black hole. On the other hand, we assume that Einstein field equation with a cosmological constant holds in the sense that the effect of effective dark matter is encoded in the total energy-momentum part; namely, G µν + Λg µν = 8πT µν . Then, from the gravitational field equations for the t − t component, we obtain yielding the solution The third term is due to the apparent dark matter effect, while the last term is the contribution due to the cosmological constant. We can perform a series expansion around x = 1/λ, yielding It is, therefore, natural to define the true or the physical mass of the black hole to be M = M (1 + α) and write the solution in terms of the physical mass At this point, it is interesting to notice that the last equation can be linked to the spherical black hole solutions in de Rham, Gabadadze and Tolley (dRGT) massive gravity [112] for ζ = 0 and γ = Mα/λ 2 . Further, we can neglect the term O(α/λ 2 ) in (64) and we obtain the Kottler spacetime, i.e., Schwarzschild black hole with a cosmological constant Following [71], it is easy to show that the shadow radius for the metric (64) can be written as In other words, such a change is outside the scope of the present technology. We can approximate the shadow radius in both cases to be R SH ≃ 3 √ 3GM BH /c 2 . Now, it is well known that the real part of quasinormal modes ω ℜ in the eikonal limit is related to the shadow radius of BHs via ω ℜ = lim l≫1 1 RSH (l + 1/2) [113][114][115][116][117][118] where l is the angular node number. This correspondence is achieved based on the geometric-optics correspondence between the parameters of a quasinormal mode and the conserved quantities along geodesics. This connection allows the testing of gravitational waves with the nextgeneration Event Horizon Telescope [119,120]. However, here we see that the frequency of the quasinormal modes emitted by a perturbed black hole in the eikonal limit will depend on the effect of cosmological constant and apparent dark matter. In particular, we obtain the following relationα SH = lim l≫1 , respectively. These relations are valid in specific conditions, i.e. the eikonal regime and the low-redshift limit.
In what follows, we will apply Eq. (57) to compute the angular radius, assuming known black hole mass and Yukawa parameters.

V. CONCLUSIONS
In the present work, we extracted observational constraints on the Yukawa cosmological model. In this scenario, dark matter appears effectively, and a relation exists between dark matter, dark energy, and baryonic matter. In particular, the effective dark matter is attributed to the long-range force between the baryonic matter particles. Such a Yukawa-like gravitational potential modifies Newton's law of gravity in large-scale structures. It is characterized by the coupling parameter α and has a graviton with non-zero mass (which is inversely related to the wavelength parameter λ). We used SNe Ia and OHD observational data, and we found within 1σ CL the best-fit parameter λ = 2693 +1191 −1262 Mpc and α = 0.416 +1.137 −0.326 , respectively. With these values, we acquire the value m g = 4.233 +3.735 −1.298 × 10 −69 kg, or, equivalently, m g = 2.374 +2.095 −0.728 ×10 −42 GeV, for the graviton mass, complementing other significant findings related to Yukawa constraints from different systems [103][104][105][106][107][108][109][110][111]. Additionally, we found a black hole solution and a relation between the dark matter/dark energy density parameters and the angular radius of black hole shadows. These equations allow us to constrain the graviton mass directly from the EHT results for Sgr A and M87 supermassive black holes. We can further use the following Gaussian likelihoods L Shadow ∝ exp (−χ 2 Shadow /2), where χ 2 shadow = i=1 [α observed SH −α theory SH /σα ,i ] 2 , and the modify the Gaussian likelihood for the joint analysis SNe Ia+OHD+Shadow as L joint = L SNe + L OHD + L shadow . We will consider and explore this possibility in a separate project.
Our research combines two systems from varying scales, successfully integrating cosmological constraints on black hole shadows.
That establishes a multimessenger constraint on gravity, modelled explicitly as the Yukawa potential. Our application mainly focuses on constraints related to late-time cosmology. However, constraints from higher redshifts, such as those derived from Cosmic Microwave Background (CMB) and Big Bang Nucleosynthesis (BBN), could effectively restrict the Yukawa parameters in future research.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.