Yukawa Gravity: Exponential Screening

• Yukawa gravity is a gravitational model in which the classical 1/r potential is modified by an exponential term, effectively introducing a scale-dependent force.
• It emerges as the weak-field limit of analytic f(R) gravity and is applied to explain phenomena from finite graviton mass effects to galaxy rotation curve anomalies.
• Observational and experimental studies reveal its impact on orbital precession, strong-field metrics, and laboratory tests, while parameter degeneracies and environmental factors complicate its interpretation.

Searching arXiv for recent and foundational papers on Yukawa gravity relevant to the provided corpus. Yukawa gravity denotes a class of gravitational models in which the Newtonian $1/r$ potential is modified by an exponentially screened term, typically written in the form

$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$

or equivalently with alternative strength parameters such as $\alpha$, $\beta$, or $\kappa$, depending on convention. In this framework, the parameter controlling the amplitude of the correction is dimensionless, while the length scale $\lambda$ or $\Lambda$ sets the range over which the modification is important. Across the literature, Yukawa gravity appears in several distinct roles: as the weak-field limit of analytic $f(R)$ gravity, as a phenomenological fifth-force parameterization, as a proxy for finite graviton mass, as a cosmological screening law in an expanding Universe, and as a laboratory-scale non-Newtonian interaction. Its observable consequences include orbital precession, modified rotation curves, altered timing signals, and changes in strong-field photon dynamics (Laurentis et al., 2018).

1. Definitions and parameterizations

The most widely used phenomenological form of the Yukawa-corrected potential is

$$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$

with $M$ the central mass, $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$0 a dimensionless coupling, and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$1 the interaction range. The same structure appears in several notational variants, for example with $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$2, $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$3, or $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$4 replacing $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$5, and with $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$6 replacing $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$7. In all such cases, $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$8 or the equivalent vanishing-coupling limit recovers the Newtonian potential $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$9. The correction becomes important when the orbital or geometric scale is comparable to the Yukawa range, whereas for $\alpha$0 the exponential is suppressed and for $\alpha$1 the modification can often be expanded perturbatively (Borka et al., 2013).

A recurrent subtlety is that the Yukawa potential is not merely a Newtonian potential plus a small short-range perturbation. In the $\alpha$2 weak-field form used in relativistic orbital studies, the prefactor $\alpha$3 rescales the $\alpha$4 term itself, so the potential behaves as if generated by an effective mass $\alpha$5 even before the exponential piece becomes numerically important. This implies that Yukawa gravity generically changes both the normalization and the radial profile of the force law. A plausible implication is that observational constraints on the Yukawa amplitude can be entangled with independent uncertainties in source mass or mass modeling.

The associated force law for a central potential $\alpha$6 is

$\alpha$7

showing explicitly that the Yukawa correction modifies the inverse-square law by a scale-dependent factor. In some studies the same range parameter is identified with a graviton Compton wavelength through

$\alpha$8

with the convention depending on the paper. This convention dependence is important because the two expressions differ by a factor of $\alpha$9; the literature discussed here uses both forms, and each must be interpreted in its own stated normalization.

2. Relation to analytic $\beta$0 gravity and weak-field relativity

A major theoretical realization of Yukawa gravity arises in analytic $\beta$1 gravity. Starting from the action

$\beta$2

and expanding

$\beta$3

around a Minkowski background, the weak-field solution yields a Yukawa-like metric potential with a characteristic range

$\beta$4

together with the identification

$\beta$5

In that construction, the Yukawa scale is generated by the extra $\beta$6 degree of freedom, and the gravitational potential becomes

$\beta$7

Thus the Yukawa form is not inserted ad hoc, but emerges as the weak-field limit of a higher-derivative gravitational action (Laurentis et al., 2018).

The same connection was used earlier in galaxy modeling, where the Yukawa-like $\beta$8 potential was written as

$\beta$9

with $\kappa$0 and $\kappa$1 in the Taylor expansion $\kappa$2. That study emphasized that $\kappa$3 can be interpreted as the effective length scale of an extra scalar degree of freedom and remarked that $\kappa$4 may also be seen as a Compton length. It also noted that the condition for a real Yukawa exponent in the phenomenological parametrization corresponds to $\kappa$5, precisely the interval that can flatten galaxy circular-speed curves in their framework (Napolitano et al., 2012).

Relativistic formulations of Yukawa gravity are not unique. Some analyses retain distinct weak-field metric potentials $\kappa$6 and $\kappa$7, while others approximate $\kappa$8 when the fractional difference is only a few percent for $\kappa$9. More recent post-Newtonian work for S-stars argues that replacing the full Yukawa metric by a Newtonian Yukawa potential or by the simplified $\lambda$0 ansatz can be inaccurate at the precision needed for Galactic-center orbital tests. In that treatment, the weak-field metric takes the form

$\lambda$1

and the resulting post-Newtonian Yukawa equation of motion contains both Schwarzschild-like terms and mixed PN–Yukawa terms (Tan et al., 2024).

3. Orbital dynamics, precession, and binary systems

The cleanest classical and relativistic signature of Yukawa gravity in bound systems is apsidal precession. In the relativistic $\lambda$2-motivated treatment, geodesic motion in the Yukawa-deformed metric yields a modified periastron advance that reduces exactly to the general relativistic expression

$\lambda$3

when $\lambda$4. Under the short-orbit approximation $\lambda$5, the Yukawa-modified advance is given in closed form by

$\lambda$6

This shows that the precession depends on the central mass, semimajor axis, eccentricity, and Yukawa parameters, with tighter and more eccentric orbits amplifying the effect. The same study stresses, however, that this analytic formula is valid only for systems whose semi-major axis is much smaller than the Yukawa length; when orbital scales approach $\lambda$7, the full geodesics must be integrated numerically (Laurentis et al., 2018).

At the purely perturbative level, a central Yukawa correction produces no secular drift in $\lambda$8, $\lambda$9, $\Lambda$0, or $\Lambda$1, but it does generate a secular advance of periastron. For a drag-free Earth satellite with perturbing potential

$\Lambda$2

the orbit-averaged disturbing function is

$\Lambda$3

leading to

$\Lambda$4

That analysis explicitly notes that because the Yukawa perturbation is central, the angular momentum direction is conserved, so any claim of Yukawa-driven inclination change in such a setup would contradict central-force dynamics (Deng et al., 2013).

Galactic-center stellar orbits provide a second major arena. Simulations of the S2 orbit under the same phenomenological potential conclude that the orbit is generally non-closing and that the precession direction depends on the sign of the coupling: it is in the same direction as general relativity for $\Lambda$5 and $\Lambda$6, and in the opposite direction for $\Lambda$7. This sign dependence makes S2 useful as a discriminator between relativistic precession, extended-mass retrograde precession, and Yukawa deviations (Borka et al., 2013).

A further complication is environmental degeneracy. When an extended mass distribution near Sgr A* is included, its retrograde contribution competes with the prograde Yukawa contribution. In that case, matching the total precession to the general relativistic value requires smaller $\Lambda$8 as the assumed extended density increases. This was quantified for bulk densities $\Lambda$9, for which the inferred $f(R)$0 at fixed $f(R)$1 drops sharply relative to the point-mass case. This suggests that orbital precession constraints on Yukawa gravity are only as robust as the modeling of extended matter near the central source (Jovanović et al., 2021).

Binary systems introduce an additional layer of degeneracy because the Yukawa correction enters together with the total mass $f(R)$2. In the $f(R)$3 binary extension, the periastron advance rate retains the structure of the GR expression multiplied by a Yukawa correction factor, but the extra Yukawa parameters are degenerate with the masses. The practical consequence is that, unlike in GR where two post-Keplerian parameters can determine two masses, Yukawa or $f(R)$4 binary timing generally requires either independent mass information or multiple post-Keplerian observables to break the degeneracy (Laurentis et al., 2018).

4. Galactic dynamics and the dark-sector interpretation

Yukawa gravity has been repeatedly studied as an alternative to particle dark matter in galaxies. In elliptical galaxies, a Jeans analysis using the $f(R)$5-motivated Yukawa potential

$f(R)$6

was applied to NGC 3379, NGC 4374, and NGC 4494. The best-fit values were $f(R)$7 kpc and $f(R)$8, with fit quality comparable to that of dark-halo models and orbital anisotropies similar to those inferred in $f(R)$9CDM analyses. That work concluded that the Yukawa-like correction can reproduce the observed stellar kinematics of the three ellipticals without adding an explicit dark halo, though it also emphasized the small sample size and modeling simplifications (Napolitano et al., 2012).

For spiral galaxies, two distinct Yukawa programs appear. One adopts a direct modification of the baryonic gravitational field and interprets the extra term as an apparent dark matter contribution. In this approach, the circular speed satisfies

$$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$0

so one may define an apparent dark component

$$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$1

Applied to the Milky Way and M31, this framework found $$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$2 and $$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$3 kpc for the Milky Way, and $$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$4 and $$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$5 kpc for M31 in the original dataset. The Milky Way fit was favored over $$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$6CDM by $$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$7, while the original M31 dataset favored $$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$8CDM with $$\Phi(r)=-\dfrac{GM}{(1+\delta)r}\left[1+\delta e^{- \left(\dfrac{r}{\Lambda} \right)} \right],$$9; with a newer M31 dataset, the discrepancy weakened to $M$0. The same paper argued that a MOND-like acceleration scale

$M$1

emerges approximately from the Yukawa correction, yielding $M$2 at $M$3 for representative parameters (D'Agostino et al., 2024).

A second galactic program modifies the interaction primarily in the dark-matter sector. In a Milky Way rotation-curve fit where the Yukawa term acts only in baryon–dark-matter gravity, the modified one-particle potential is

$M$4

with the Yukawa correction applied only to the NFW dark halo. Outside the inner $M$5 kpc, the data were found to be consistent with Newtonian gravity plus a standard NFW halo; if a Yukawa term exists, the preferred sign is negative, corresponding to an effectively repulsive baryon–dark-matter interaction with short range, and the allowed parameters follow an approximate degeneracy curve

$M$6

This implies $M$7 and shows that the data constrain a combination of coupling and range more strongly than either parameter separately (Henrichs et al., 2020).

A related 2025 Bayesian study extended Yukawa rotation-curve modeling to four scenarios—no-dark-matter Yukawa gravity, non-trivial dark-matter coupling, fully modified gravity, and Newtonian gravity—and concluded that the Milky Way can show extremely large Bayes factors in favor of short-range Yukawa models, but that these cases risk overfitting because Yukawa-modified dark matter can mimic baryonic kinematics. In M31, where photometric priors are stronger, Bayes factors instead favor Newtonian gravity. This suggests that realistic priors and independent morphology constraints are essential in any attempt to infer Yukawa parameters from rotation curves (Hassan et al., 6 Aug 2025).

5. Cosmological and strong-field formulations

Not all versions of Yukawa gravity arise as beyond-GR modifications in the usual sense. One cosmological formulation argues that, within standard general relativity on an expanding FLRW background, the effective gravitational interaction between nonrelativistic point masses is Yukawa-screened: $M$8 with screening length

$M$9

In this interpretation, the homogeneous cosmological matter background itself screens gravity. The paper’s main consistency test is the homogeneous-Universe limit, where the net peculiar gravitational acceleration should vanish. It argues that the Yukawa law passes this test exactly while a naive Newtonian cosmological force law double-counts the matter contribution. The present-day screening length was quoted as $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$00, suggesting a cosmological-scale Yukawa range rather than an astrophysical fifth force (Eingorn, 2017).

A different cosmological construction embeds a regularized Yukawa potential in a Verlinde-style entropic framework. There the effective Newton constant becomes

$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$01

the graviton mass is identified through

$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$02

and the model is used to reinterpret both dark matter and dark energy as emergent effects of Yukawa-modified baryonic gravity. That work estimated $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$03 Mpc, $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$04, $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$05 kg, and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$06, and argued that the same framework can reproduce a MOND-like acceleration scale of order $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$07 (Jusufi et al., 2023).

In strong gravity, Yukawa deformations have been applied directly to black-hole metrics. A static spherically symmetric Yukawa black hole was modeled with

$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$08

so that positive $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$09 increases the shadow size and negative $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$10 decreases it. Using EHT shadow measurements, the resulting constraints were: for Sgr A*, $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$11 for $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$12 AU with the Keck prior and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$13 with the VLTI prior; at $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$14 AU, the bounds weaken to $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$15 and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$16, respectively. For M87*, the bound $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$17 applies for $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$18 AU. No significant deviation from GR was found (Tan et al., 5 Oct 2025).

A plausible synthesis of these cosmological and strong-field strands is that “Yukawa gravity” is not a single theory but a shared functional form arising in very different contexts: screening by cosmic background matter, effective massive degrees of freedom in extended gravity, entropic-gravity phenomenology, and strong-field metric deformations.

6. Experimental, timing, and laboratory probes

Solar-System-scale laboratory and satellite experiments probe Yukawa gravity in a very different regime from Galactic-center dynamics. For a drag-free Earth satellite with $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$19 m and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$20, the predicted secular periastron advance from a Yukawa perturbation spans from several nano-arcseconds per year to hundreds of arcseconds per year over the parameter domain

$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$21

The same work showed that the Yukawa contribution to onboard clock time transfer has amplitude

$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$22

amounting only to $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$23–$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$24 s in the studied range. Its central conclusion was that orbital dynamics, especially secular periastron advance, are much more detectable than timing/redshift signals for such a satellite mission (Deng et al., 2013).

Pulsar timing around Sgr A* provides a much more sensitive relativistic probe. A forecast based on $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$25 years of observations, $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$26 TOAs, and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$27 found that a pulsar–SMBH system could improve current Yukawa tests when the interaction range varies between $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$28 and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$29 AU. In that range, the predicted sensitivity reaches $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$30–$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$31, and the graviton-mass interpretation yields $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$32, with favorable cases approaching $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$33. The analysis also emphasized the $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$34–$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$35 degeneracy when $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$36 is much larger than the orbital scale, as well as residual absorption of slowly varying signals by the pulse-number offset $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$37 (Dong et al., 2022).

S-star astrometry remains an important intermediate-scale probe. A 2022 MCMC study using the measured closeness of S-star precession to the Schwarzschild prediction found that current S-star data prefer small Yukawa strength and interaction range around $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$38 AU in joint fits, with full-sample values

$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$39

That work also introduced the scale variable

$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$40

as a classification criterion for gravitational systems in Yukawa gravity, since the inferred $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$41 scales with $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$42 when the precession benchmark and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$43 are fixed (Jovanović et al., 2022).

At meter scales, laboratory searches have also entered the Yukawa-gravity domain. A 2026 study proposed a test using the CHRONOS torsion-bar detector with a differential gravitational calibrator. The setup cancels the full Newtonian torque from two rotating source systems while retaining a residual Yukawa torque because the exponential factor breaks the exact Newtonian cancellation. The projected sensitivity reaches

$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$44

with a systematic floor reached at

$$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$45

The dominant limitation is residual Newtonian torque from imperfect cancellation, driven mainly by geometric uncertainties in the source masses and baselines (Inoue et al., 13 Apr 2026).

Together, these studies show that the experimentally accessible Yukawa range spans many orders of magnitude, from meters to kiloparsecs to gigaparsecs, but the nature of the constraint changes with scale: laboratory probes test static or modulated deviations from inverse-square gravity, orbital systems measure secular phase drift and precession, and cosmological or strong-field probes test screening laws or null geodesics.

7. Degeneracies, controversies, and scope

Several recurring issues complicate the interpretation of Yukawa gravity. The first is parameter degeneracy. Orbital analyses repeatedly find that the Yukawa amplitude and range are only partly separable: in S2 fits, $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$46 and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$47 are strongly correlated for $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$48, and in Milky Way rotation-curve fits $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$49 and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$50 are constrained mainly along a degeneracy curve rather than individually. In the $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$51 regime, much of the effect is absorbed into an effective rescaling of the force, making $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$52 weakly constrained unless the orbit or geometry probes radii comparable to the Yukawa length.

The second issue is environmental contamination. In Galactic-center dynamics, extended mass distributions, other stars, and frame systematics can mimic or obscure Yukawa signatures. In galaxy rotation curves, uncertainties in bulge and disk morphology can be traded against Yukawa parameters, producing apparent evidence for modified gravity that may instead reflect overfitting. In pulsar timing, incomplete orbital coverage and omitted spin or environment effects can degrade sensitivity.

The third issue is sign convention and internal inconsistency. One relativistic periastron paper notes a mismatch between its textual description and numerical table concerning the sign of the Yukawa correction for $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$53, warning that explicit formulae and direct numerical evaluation should be trusted over prose summaries. Other papers contain typographical corruption in equations or notation shifts between $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$54 and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$55, $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$56 and $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$57, or different graviton-mass conventions. This does not invalidate the shared physical content, but it requires careful normalization when comparing results across the literature.

Finally, the term “Yukawa gravity” itself is broader than a single theory. In some papers it means a phenomenological fifth-force potential; in others it is the weak-field limit of analytic $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$58 gravity; in still others it denotes cosmological screening within standard GR, a massive-graviton proxy, or a strong-field metric ansatz. A related but distinct usage appears in matter-sector studies of Yukawa interactions that are not modifications of gravity itself, and in source-based constructions such as Yukawa-Casimir wormholes, where the field equations remain Einstein–Cartan rather than Yukawa-modified gravity. This suggests that the defining feature of Yukawa gravity is not a unique fundamental Lagrangian, but the repeated emergence of the exponentially screened kernel $$\Phi(r)=-\frac{GM}{(1+\delta)r}\left(1+\delta e^{-r/\lambda}\right),$$59 as an effective descriptor of scale-dependent gravitational behavior.