Advances in Mathematical Physics
Volume 2019 (2019), Article ID 6765827, 10 pages
https://doi.org/10.1155/2019/6765827

Research Article

Yukawa Potential Orbital Energy: Its Relation to Orbital Mean Motion as well to the Graviton Mediating the Interaction in Celestial Bodies

Connor Martz,¹ Sheldon Van Middelkoop,² Ioannis Gkigkitzis,³ Ioannis Haranas,⁴ and Ilias Kotsireas⁴

¹University of Waterloo, Department of Physics and Astronomy, Waterloo, ON, N2L-3G1, Canada
²University of Western Ontario, Department of Physics and Astronomy, London, ON N6A-3K7, Canada
³NOVA, Department of Mathematics, 8333 Little River Turnpike, Annandale, VA 22003, USA
⁴Wilfrid Laurier University, Department of Physics and Computer Science, Waterloo, ON, N2L-3C5, Canada

Correspondence should be addressed to Ioannis Haranas; iharanas@wlu.ca

Received 25 September 2018; Revised 25 November 2018; Accepted 4 December 2018; Published 1 January 2019

Academic Editor: Eugen Radu

Copyright © 2019 Connor Martz et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Research on gravitational theories involves several contemporary modified models that predict the existence of a non-Newtonian Yukawa-type correction to the classical gravitational potential. In this paper we consider a Yukawa potential and we calculate the time rate of change of the orbital energy as a function of the orbital mean motion for circular and elliptical orbits. In both cases we find that there is a logarithmic dependence of the orbital energy on the mean motion. Using that, we derive an expression for the mean motion as a function of the Yukawa orbital energy, as well as specific Yukawa potential parameters. Furthermore, various special cases are examined. Lastly, expressions for the Yukawa range and coupling constant are also derived. Finally, an expression for the mass of the graviton mediating the interaction is calculated using the expression its Compton wavelength (i.e., the potential range ). Numerical estimates for the mass of the graviton mediating the interaction are finally obtained at various eccentricity values and in particular at the perihelion and aphelion points of Mercury’s orbit around the sun.

1. Introduction

Any scientist will agree that Einstein’s general relativity theory (GR) is one of the most mathematically elegant theories invented in the human history. Even though the theory explains many physical phenomena, it is unable to shed light on the problem of the observed accelerating universe. To do that, GR introduces a cosmological constant lambda as well as the so-called dark energy. Various gravitational theories exist today which try to explain the observed acceleration of the universe [1]. These gravitational theories are nonsymmetric, scalar-tensor, quantum gravitational, or theories of gravity, etc. As a common denominator in the weak limit, all theories result in a Yukawa type of gravitational potential. In a paper by Chan [2], the authors put forward the idea that observational evidence for the existence of cold dark matter particles in the cores of dwarf galaxies could be explained through the interaction of a Yukawa potential. Therefore, the experimental and observational search for such deviation might result in new type of physics [3].

In a recent paper by Haranas et al. (2018), the authors examine dust particle orbits under the influence of Poynting-Robertson effect in which Newtonian gravity has been modified by a Yukawa term. Similarly, in Mukherjee and Sounda [4], the authors investigate the orbits resulting from various coupling constants (alpha) of a Yukawa correction to the Newtonian potential. Quite often in celestial mechanics, a Yukawa-type potential is proposed to modify the Newtonian gravity [5–9] (Iorio 2007), and its effect on various gravitational, astrophysical, and orbital scenarios is examined. In this contribution, we examine the time rate of the Yukawa orbital energy of circular and elliptical orbits, and from that we derive relations for the orbital energy as a function of mean motion of the orbiting body. To relate and to the orbital parameters of the secondary, expressions are derived that relate them to various orbital parameters. In particular, considering the expression for λ, we obtain a Lambert function that relates the mass of the graviton along the orbit of the secondary to the Yukawa parameters, eccentric anomaly, orbital energy, and eccentricity. This is done using the already derived expression for lambda and substituting it into the corresponding equation for the range of graviton and then solving for its mass . It is well known that if gravitation is propagated by a massive field, the velocity of the gravitational waves (gravitons) will depend upon their frequency, and the effective Newtonian potential will have a Yukawa form, i.e., , where is the graviton Compton wavelength. Today’s research for the mass of the graviton includes both theoretical and observational work. For example, in Stavridis and Will [10] the authors try to bound the graviton mass using gravitational effects and their effect in the spin precessions of massive hole binaries. Similarly, in Mureika and Mann [11], the authors use an entropic gravity approach to estimate a bound for the mass of the graviton. Finally, in Zacharov et al. [12], the authors consider Yukawa gravity interactions of S2 star orbits near the galactic plane to improve expectations for graviton mass bounds. At this point we must say that in today’s gravity research various methods have been employed in the determination of the graviton mass. Our motivation for paper emanates from the fact that this work can serve as another possible observational test in setting solar system as well as binary system bounds on graviton mass , where the bound depends on the mass of the source, which in this case is a sun like type of star.

2. The Yukawa Potential

Let us consider a two-body problem, where a secondary body of mass orbits under the influence of a primary body of mass . With the potential being central, the two-body problem can be reduced to a central-force problem, and the motion of the secondary body can be examined. The effects of gravity on the secondary in the presence of a Yukawa correction can be described by the modified potential energy per unit mass [13].

In (1), denotes the distance between the centers of the two bodies, is the Newtonian gravitational constant, , where and are the coupling constants of the new force to the bodies relative to the gravitational one, and is the range of this interaction (ibid. 2016). Next, let us now write down an expression of the orbital energy by making use of the virial theorem which states that for a bounded unperturbed system the following relation holds [14]:

where is the time average on the system’s kinetic energy and the time average of the system’s potential energy. Therefore, conservation of energy implies that the total energy is equal to (ibid. 2002)

Anticipating small perturbations, we can add the average value of the perturbing term to the average of the Newtonian potential to that of the perturbative term making also use that the semimajor axis is the average value of the radial distance along the orbit; we can further write that as

Now rewriting (1) only the Yukawa orbital energy we have is

Next differentiating (5) with respect to time we obtain that the total time rate of the energy of the Yukawa correction is

3. Circular Orbits

To examine the case of circular orbits, namely, we let the eccentricity = 0 and r = a, (6) becomes

Next, following Haranas et al. [9], we consider the unperturbed relative orbit of the secondary body, in this case a Keplerian ellipse. If is the semimajor axis, is its eccentricity and is its mean motion. On the unperturbed Keplerian ellipse, this law is expressed as [15]

Differentiating (8) w.r.t. time, the rate of change of the mean motion can then be expressed as

from which we obtain that

Therefore, substituting (10) in (7) and simplifying, we obtain

Next defining the potential energy at semimajor to be

the Newtonian part of the potential, (11) can written as follows:

Integrating, we obtain the orbital energy dependence on the mean motion to be

On the other hand, to obtain the dependence of the mean motion on the Yukawa parameters and orbital energy , we use (14) and solving for we obtain

In the case where the range of the potential equals the semimajor axis = a, (15) becomes

At this point we can see that value of the mean motion depends on the energy difference . If then and therefore the exponent remains negative and the mean motion reduces exponentially with the energy difference. On the other hand, if and , the exponent is positive and therefore the mean motion increases. Finally, if , then (15) for the mean motion becomes

decreasing for and increasing for . Moreover, solving for the range of the potential we find that

where is the Lambert function of the indicated argument. Similarly, solving for the coupling constant we obtain

In the case where we find that the range lambda of the Yukawa potential is equal to

and the coupling contract becomes

The semimajor axis is equal to the range of potential and according to (20) must satisfy the following equation:

from which we find that the gravitation satisfies the following equation:

Consequently, the energy difference satisfies the following equation:

Using (24) we find that the mean motion for must independently satisfy the equation

In particular if , (25) satisfies the equation

Furthermore, using (19) if we find that

Using, (19) it is impossible to find an expression for since the resulting equation takes the form

Moreover, if in (28) we impose the condition λ = a, we obtain that

4. Elliptical Orbits

Similarly, for elliptical orbits we have the following expression:

where is the orbital radial vector. Using (8) we find that

The resulting simplified equation is expressed in terms of the eccentric anomaly, i.e., where is the orbital eccentricity of the secondary and is the eccentric anomaly; equation (31) becomes

Next substituting (33) in (30) we find that

from which we obtain the following equation for the change of orbital energy w.r.t. the mean motion :

Integrating and rearranging (34) with respect to , we obtain that

Solving (35) for the mean motion we find that

As a check of (36), let us note that if e = 0 we can easily recover (15) for the circular orbits. In the case where , (36) simplifies in the following way:

At this point we find that value of the mean motion again depends on the energy difference . If again then implies that the mean motion reduces exponentially to the energy difference. On the other hand, if and the motion increases exponentially. Finally, if , (36) for the mean motion becomes

Next, if simultaneously and , we find that the mean motion of the orbiting body takes the following form:

Furthermore, if and and solving for the mean motion n we find that

Finally, if then we find that

Going back to (36) which is the most general equation and solving for and , respectively, we obtain that

Similarly, for we find that

5. The Yukawa Potential Range and the Particle Mediating the Interaction

At this point I would like to bring the readers’ attention to an issue that might have some of them concerned. This is the issue of closed bounded orbits under the light of Bertrand’s theorem which proves that only a Newtonian as well as a Hook type of potential results in closed bounded orbits. This is an issue that the authors of this paper will extensively deal with in a paper to come soon. As a preliminary result we say that all the applications of the Yukawa potential in our research are related to the solar system and within the boundaries of the solar system dimensions where the following condition, namely, , is true [16]. The parameter is the Compton wavelength of the exchange particle, which in the present case is a graviton. In other words, the range of the Yukawa potential is always larger than any distance away from the primary star the sun in our case that the secondary body orbits. Our preliminary analysis of nearly circular orbits has shown that the ratio of a harmonic oscillator with frequency oscillating around a point situated at to that of the orbital frequency at is given the ratio . In the case of the solar system with sun as the primary gravitating body producing the interaction, and since , we obtain to a second order approximation that and therefore the orbits are bounded and closed. As a support of our findings we refer to a paper by Mukherjee and Sounda [4]. In this paper the authors find closed and bounded orbits for various values of the coupling constant of the Yukawa potential.

One of the primary motivations for today’s interest in non-Newtonian gravity emanates from our interest to probe and understand long-range forces. This interest addresses the following question: Why a certain class of theories predict the existence of gravity - like interactions. At this point we must say that many various models in which gravity is unified with other forces have been studied in the recent years. These models are formulated along the consideration that the product μ,

determines the mass of the new field [17]. The Compton wavelength (or “range") λ associated with the field characterises the distance at which the corresponding field acts and is given by the equation below:

In Haranas et al. [13] and Moffat and Toth [16], the strength of the potential near the sun in a point source scenario is given by the following equation:

where and . Similarly, the range can be written as follows:

and kpc⁻¹. Therefore equating (42) and (45) we obtain that the mass of the mediating of the interaction graviton is

and finally using (46) the equation takes the form

In the case of circular motion, i.e., e = 0, (49) simplifies to

Similarly, in the case of circular motion and if , (50) becomes

Equations (48) to (51) give the mass of the graviton field mediating the interaction, i.e., a massive graviton field. According to (48) we see that the mass of graviton continuously varies along elliptical orbit with a maximum value at the perihelion/periastron point (binary star scenario). Its mass also depends on the ratio of as well as on ratio of the initial to final mean orbital motion and finally on the mass of the primary body producing the Yukawa field, Moffat and Toth [16].

6. The Extended Source Case

In the case of an extended mass distribution following Haranas et al. [13], we say that the equation of Yukawa potential must be modified. For an extended matter distribution and in relation to the MOG field, there are not presently any solutions. In this case, authors treat the problem phenomenologically, seeking to find an effective mass distribution that could be used in (46) and (47). What this function does simply determines an “effective mass" which determines and . This way, the gravitational influence of matter located at a point of distance on the test particle located at is calculated. Thus, the point source function can simply yield a mass proportional to the volume if the distribution is constant. Therefore, following Haranas et al. [13], we can write

where the coefficient is to be determined by observation, for example, comparison with Bullet Cluster Data. For a mass density function, ; then we will have that . Furthermore, in the case of a nonconstant density , the coupling constant and range of the potential become (ibid. 2016)

and

In the case of an extended mass distribution, calculation of α and λ will result in better estimates of the two corresponding parameters and therefore a better estimate of the effect. With reference to Haranas et al. [13] in the case of an extended source we can write the gravitational potential in the following way:

which finally integrates to

where

and , and it has the following limits:

In a paper by Branco et al. [18], by examining the component of the metric the authors identify a Newtonian potential plus a Yukawa perturbation of the form (ibid. 2014)

where is the characteristic length, is the strength of the Yukawa addition to the field, and is the form factor defined as follows [18]:

which integrates to

Equation (62) admits the following limiting cases (ibid. 2014):

For a sun like central body (ibid. 2014) the more accurate NASA model has been used for the density of the sun that obeys the following conditions and . Thus, the density function is of the form

Integrating (65) the authors obtained the following function for the form factor (http://spacemath.gsfc.nasa.gov/(2014)):

with limiting cases given by , and with the limiting cases

In the case of an extended source we only consider a sun like star. In this case according to Haranas et al. [13] and Moffat and Toth [16] m will imply and therefore . Therefore, in the case of a sta,r our analysis above does not need to involve the factor function.

7. Discussion and Numerical Results

To obtain numerical results, we first start with circular orbits. Using (13) and the relation to eliminate from the equation, we obtain that

Furthermore, using the fact that or angular momentum per unit mass and for circular orbits, i.e., e = 0, we write (69) in the following way:

In other words, we find that the rate has the units of angular momentum. Similarly, in the elliptic case, (33) results in (71) below:

where we define as “Yukawa function” the term appearing in (71) and (72) that is originally derived from the term . Assuming circular orbits and a semimajor axis equal to the semimajor axis of Mercury = 5.79 × 10⁷ km and Yukawa coupling constants α = 3.039 × 10⁻⁸ Moffat and Toth [16] and = 3.57 × 10⁻¹⁰ and range λ = 4.937 × 10¹⁵ m as in [13] in the case where we obtain that

Similarly, if we obtain that

Finally, if and using as it is given in Moffat and Toth [16], we obtain the following relation:

Similarly, in Table 1 looking at elliptical orbits, we consider Mercury’s actual orbit with semimajor axis a = 57,909,050 km, and eccentricity e = 0.205 and thus we have obtained the values shown in Table 1 for the rate as a function of eccentric anomaly.

Table 1: Table of the rate as a function of eccentric anomaly for an elliptical orbit in a Yukawa potential.

Eccentric Anomaly
°
[kgm²s⁻¹]

40
-4.738×10¹³h

60
-4.450×10¹³h

80
-4.141×10¹³h

100
-3.856×10¹³h

120
-3.622×10¹³h

140
-3.451×10¹³h

160
-3.348×10¹³h

180
-3.314×10¹³h

200
-3.348×10¹³h

220
-3.451×10¹³h

240
-3.622×10¹³h

260
-3.856×10¹³h

280
-4.141×10¹³h

300
4.450×10¹³h

320
-4.738×10¹³h

340
-4.946×10¹³h

360
-5.023×10¹³h

In Figure 1 the Yukawa function is plotted for elliptical orbits of various eccentricities, i.e., blue e = 0.01, red e = 0.205, olive green e = 0.6, and green e = 0.9, and for eccentric anomaly values during one full revolution. Similarly, in Figure 2 we plot the Yukawa function for planet Mercury as a function of eccentric anomaly and during one full revolution. Next, in Table 2 we tabulate the parameters and as related to the graviton mass at perihelion as is “felt" from an orbiting around a sun like star body. For this numerical calculation for planet Mercury, α = 3.04 × 10⁻⁸, λ = 4.937 × 10¹⁵ m, and e = 0.205. Similarly, in Table 3 we tabulate the same parameters as in Table 2 for Mercury using α = 3.57 × 10⁻¹⁰ and λ = 4.937 × 10¹⁵ m. In Table 4 we tabulate the ratio of the mediating of the interaction graviton mass at perihelion and aphelion points as a function of various eccentricities using = 138.326. In Figure 3 the mass ratio of graviton particle mediating the Yukawa interaction as is felt from Mercury at perihelion and aphelion position is plotted, as a function of eccentric anomaly for = 5.79 × 10¹⁵ m and , respectively. Here and indicate the mass of the massive graviton mediating the interaction at the perihelion and aphelion points on the orbit of the secondary. In Table 5, we tabulate the ratio of the graviton mass ratios as a function of orbital eccentricities and for = 13.83408952. Using (48) we find that at perihelion the mass of the graviton is given by

Similarly, at aphelion we find that

Moreover, since the graviton mass must be a positive number, (48) implies that the following equation must be satisfied:

This can be true if the following conditions are met for all the corresponding parameters separately. First from (77) solving for the coupling constant , we find that (83) is satisfied if

where is the exponential function. For circular orbits, (78) takes the form

Similarly, for the mean motion from Eq. (77) we find that satisfied if:

which for circular orbits becomes

Next, in relation to the final orbital energy, (77) is satisfied if

which for circular orbits takes the form

Moreover, solving for the value of satisfying (77) we obtain

which for circular orbits becomes

Next using (75) and (76) we can obtain a relation for the ratio of the graviton mediating the interaction at perihelion and aphelion correspondingly, to be

where the is given by

Next calculating we calculate the orbital energy difference , for planet Mercury using α = 3.04 × 10⁻⁸ and λ = 4.937 × 10¹⁵ m and for various values of the eccentric anomaly. We tabulate the results in Tables 2 and 3.

Table 2: Tabulation of parameters and relating to the mediating graviton mass at perihelion for a sun like star as is “felt" by the secondary i.e. Mercury for α = 3.04×10⁻⁸ and λ = 4.937 ×10¹⁵ m, e= 0.205.

Coupling constant

Eccentric Anomaly
°

J

0.0
-7.708×10²⁴

3.04×10⁻⁸
45
-7.685×10²⁴
13.83408952

90
-7.631×10²⁴

180
-7.555×10²⁴

270
-7.593×10²⁴

360
-7.708×10²⁴

Table 3: Tabulation of parameters relating to the ratio of the mediating graviton mass at perihelion for a sun like star as is “felt” from the orbiting secondary i.e. Mercury for = 3.57×10⁻¹⁰ and = 4.937 × 10¹⁵ m and for various values of the eccentric anomaly.

Coupling constant

Eccentric Anomaly
°

0.0
-9.052×10²²

3.57×10⁻¹⁰
45
-9.025×10²²
138.326

90
-8.961×10²²

180
-8.872×10²²

270
-8.916×10²²

360
-9.052×10²²

Table 4: Ratio of the mediating the interaction graviton as a function of various eccentricities as is “felt” by the orbiting secondary body = 138.326 between perihelion and aphelion positions.

Eccentricity

0.0
1.000000

0.001
1.001660

0.1
1.182141

0.2
1.400340

0.3
1.671480

0.4
2.018540

0.5
2.482622

0.6
3.142260

0.7
4.170730

0.8
6.049680

0.9
10.90870

Table 5: Ratio of the mediating the interaction particle as a function of various eccentricities as is “felt” by the orbiting secondary body for = 13.83408952.

Eccentricity

  g
  g

0.0
1.00000
1.00000×10⁻⁶²
1.000×10⁻⁶²

0.001
1.00155
1.00155×10⁻⁶²
9.984×10⁻⁶³

0.1
1.16851
1.16851×10⁻⁶²
8.558×10⁻⁶³

0.2
1.36960
1.36959×10⁻⁶²
7.301×10⁻⁶³

0.3
1.61576
1.61574×10⁻⁶²
6.190×10⁻⁶³

0.4
1.92698
1.92695×10⁻⁶²
5.190×10⁻⁶³

0.5
2.33765
2.33761×10⁻⁶²
4.277×10⁻⁶³

0.6
2.91303
2.91297×10⁻⁶²
3.433×10⁻⁶³

0.7
3.79632
3.79623×10⁻⁶²
2.634×10⁻⁶³

0.8
5.38523
5.38508×10⁻⁶²
1.856×10⁻⁶³

0.9
9.46986
9.46964×10⁻⁶²
1.056×10⁻⁶³

Figure 1: Plot of the “Yukawa function” for elliptical orbits of various eccentricities and eccentric anomaly during one full blue e = 0.01, red e = 0.205, olive green e = 0.6, green e = 0.9.

Figure 2: Plot of the Yukawa function for planet Mercury as a function of eccentric anomaly during a full revolution.

Figure 3: Mass ratio of the mediating the Yukawa interaction particle as it is felt by Mercury at perihelion and aphelion is plotted as a function of eccentric anomaly for = 5.79 × 10¹⁵ m and .

Taking the numerical values of into account as in (87) we can write the mass of the graviton at perihelion and aphelion to be

and

where is the Lambert function of the indicated argument. In Table 4 we calculate and plot the variation of the mass of the mediating particle for the fixed values of the coupling constant alpha and range lambda of the potential as a function of eccentricity e.

With reference to Ohanian and Ruffini [19] looking at (91) and referring to Yukawa gravitational potential according to relativistic quantum theory, we say that the mass of the graviton is inversely proportional to the range . Relying on the observational limit of > 10²⁴ cm we find that g. Next, assuming the graviton mass to be we tabulate the corresponding mass values of graviton at perihelion and aphelion as “felt" by the orbiting Mercury.

8. Conclusions

In this paper we have derived the rate of change of the energy of an orbiting body with respect to mean motion in the influence of Yukawa potential of coupling constant and range . We have found that the rate of change of the orbital energy of circular and elliptical orbits w.r.t. the mean motion is a logarithmic function of the mean anomaly . Furthermore, using this expression we have derived expressions for the mean anomaly as a function of the Yukawa parameters and the orbital energy difference which simplifies considerably in the case where that range of the Yukawa potential is equal to the semimajor axis, i.e., λ = a. Furthermore we have found that in the case where the energy difference is negative the orbital mean motion increases exponentially. If is positive, then mean motion decreases exponentially. In the special case where , the mean motion also decreases exponentially for > 0 and increases for . Similarly, for circular orbits we have derived expressions for the Yukawa parameters and . We have found that is given in terms of a Lambert function of argument . Similarly, for circular orbits we find that . Analogous relations have been derived for elliptical orbits. Next expressions for the mass of graviton have been derived in terms of the orbital elements of the orbiting secondary for circular and elliptical orbits, from which special conditions have been extracted for a positive graviton mass. We have found that there is no graviton mass difference effect for circular orbits, as to what the mass of the graviton between perihelion and aphelion is. Therefore, the strength of the gravitational interaction along the orbit propagates with the help of a constant mass graviton. As the eccentricity increases, the mass of the massive graviton appears to increase, being larger at perihelion where at aphelion it is always less by at least two orders of magnitude from its assumed mass, i.e., . Our motivation for this paper emanates from the hope that this work can add an extra possibility in the estimation of the graviton mass. This will be the addition of one more observational test in setting solar system as well as binary system bounds for the graviton mass in terms of orbital energy mean motion relations. This is in the case where the bound depends on the mass of the source, which in this case is a sun like type of star. Thus near the massive body we would expect the graviton mass to have its higher value where its lowest value will be at the aphelion point, respectively. In other words, the strength of the interaction of gravity will reduce and this is something we already know. The only difference is that now we can have an explanation for the strength of gravitational interaction in terms of a massive graviton mass change as the new theories predict.

Appendix

The Lambert function is defined as the inverse function of the mapping and thus solving the equation . This solution is given in the form of the Lambert function, which according to the relation satisfies the following equation: [20]. The function is also known as Omega or Product Log function.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

F. Piazza and C. Marinoni, “Models for Gravitational Interraction between Dark Matter and Baryons,” vol. 91, page 14, 2003.

H. M. Chan, “Observational evidence of the Yukawa Potential Interacting Dark Matter,” The Astrophysical Journal Letters, 2013.

E. G. Adelberger, B. R. Heckel, and A. E. Nelson, “Tests of the Gravitational Inverse-Square Law,” Annual Review of Nuclear and Particle Science, vol. 53, pp. 77–121, 2003.

R. Mukherjee and S. Sounda, “Single particle closed orbits in Yukawa potential,” Indian Journal of Physics, vol. 92, no. 2, pp. 197–203, 2018.

L. Iorio, “Constraints on a Yukawa gravitational potential from laser data of LAGEOS satellites,” Physics Letters A, vol. 298, no. 5-6, pp. 315–318, 2002.

J. R. Brownstein and J. W. Moffat, “The Bullet Cluster IE0657-558 Evidence shows Modified Gravity in the absence of Dark Matter,” vol. 23, page 3427, 2006.

J. R. Brownstein and J. W. Moffat, “Gravitational solution to the Pioneer 10/11 anomaly,” Classical and Quantum Gravity, vol. 23, no. 10, pp. 3427–3436, 2006.

I. Haranas and O. Ragos, “Yukawa-type effects in satellite dynamics,” Astrophysics and Space Science, vol. 331, no. 1, pp. 115–119, 2010.

I. Haranas, O. Ragos, and V. Mioc, “Yukawa-type potential effects in the anomalistic period of celestial bodies,” Astrophysics and Space Science, vol. 332, no. 1, pp. 107–113, 2010.

A. Stavridis and C. Will, “Bounding the mass of the graviton with gravitational waves: Effect of spin precessions in massive black hole binaries,” Physical Review D: Particles and Fields, vol. 80, no. 4, 2009.

J. R. Mureika and R. B. Mann, “Does entropic gravity bound the masses of the photon and graviton?” Modern Physics Letters A, vol. 26, no. 3, pp. 171–181, 2011.

A. Zakharov, P. Jovanović, D. Borka, and V. B. Jovanović, “Constraining the range of Yukawa gravity interaction from S2 star orbits III: improvement expectations for graviton mass bounds,” Journal of Cosmology and Astroparticle Physics, vol. 2018, no. 04, 2018.

I. Haranas, I. Kotsireas, G. Gomez, M. J. Fullana, and I. Gkigkitzis, “Yukawa effects on the mean motion of an orbiting body,” Astrophysics and Space Science. An International Journal of Astronomy, Astrophysics and Space Science, vol. 10, 2016.

H. Goldstein, C. P. Poole Jr., and J. L. Safko, Classical Mechanics, Addison Wesley, 2002.

C. D. Murray and S. F. Dermott, Solar System Dynamics, Cambridge University Press, Cambridge, UK, 1999.

J. W. Moffat and V. T. Toth, “Fundamental parameter-free solutions in modified gravity,” Classical and Quantum Gravity, vol. 26, no. 8, 2009.

E. Fischbach and C. L. Talmadge, The Search for Non-Newtonian Gravity, Springer, Berlin, Germany, 1999.

N. C. Branco, J. Paramos, and R. March, “Perturbation of the metric around a spherical body from a nonminimal coupling between matter and curvature,” Physics Letters. B. Particle Physics, Nuclear Physics and Cosmology, vol. 735, pp. 25–32, 2014.

H. C. Ohanian and R. Ruffini, Gravitation and Spacetime, Cambridge University Press, Cambridge, UK, 2012.

D. Veberič, “Lambert W function for applications in physics,” Computer Physics Communications, vol. 183, no. 12, pp. 2622–2628, 2012.