© The Authors 2026

1 Introduction

High-precision measurements of the stellar dynamics in the Galactic centre have allowed for direct tests of modified theories of gravity (Borka et al. 2013; Della Monica et al. 2021; Dialektopoulos et al. 2019; Rahman & Sen 2019; Zakharov 2018) and dark mass (DM) (DM includes dark matter and any other faint source (Heißel, et al. 2022; Lechien et al. 2024; Yuan et al. 2022)). Namely, the contribution of Schwarzschild and mass pre-cessions on the orbit of S2 around Sgr A* has been studied (Heißel, et al. 2022), and a method for constraining the DM was set up assuming two different density profiles around the super-massive black hole. A more general density profile has been used as well to reconstruct DM (Lechien et al. 2024), which requires significantly more observations to achieve the same level of accuracy in setting bounds on the mass, however.

Constraints on the parametrised post-Newtonian (PPN) parameters β_PPN and γ_PPN were also imposed using the orbits of S-stars (Gainutdinov 2020), and constraints were forecast for future S2 and S62 observations (de Mora Losada et al. 2025).

The PPN formalism is independent of a particular modified theory because it provides a general framework that parametrises deviations from GR. When strong constraints are to be placed on a particular theory of gravity, however, PPN is not optimal. As a more general framework, it considers more parameters than in the case of a particular theory. For instance, based on the orbital precession of the S2 star in particular, scalar-tensor-vector gravity (STVG; see below) has been subject to tests (Della Monica et al. 2021) under a fully relativistic framework. Even though the latter is formally complete, it is not necessarily the most intuitive and handy approach for investigating effects of a modified theory. When far enough from the black hole, as is the case for the currently known S-stars, the post-Newtonian framework offers a greater flexibility, such as an easier way to add perturbative accelerations, and the ability to investigate the impact of these perturbations, as well as of other effects, in isolation.

In an attempt to construct a modified theory of gravity that can explain the missing mass in the galaxies without dark matter, John Moffat constructed the STVG, also known as modified gravity (MOG) (Moffat 2006). In this theory, the gravitational constant G, a vector field coupling constant ω as well as its mass µ_ϕ, are dynamical scalar fields that can vary in space and time. It should be noted that unlike modified Newtonian dynamics (MOND) (Milgrom 1983), MOG theory is relativistic by construction. It can therefore be tested in a wide range of astrophysical and cosmological environments. Namely, MOG passes the wide binary test (Moffat 2024) and accounts for galactic dynamics without postulating dark matter (Moffat & Rahvar 2013; Brownstein & Moffat 2006; Moffat & Rahvar 2014; Davari & Rahvar 2020; Moffat & Toth 2023; Green & Moffat 2019; Moffat & Toth 2015). It can also explain the Bullet cluster (Israel & Moffat 2018; Brownstein & Moffat 2007), reproduces lensing solely with baryon mass (Moffat & Toth 2021), and fits or potentially explains cosmological data (Moffat 2026a,b; Davari & Rahvar 2021; Moffat 2020). Furthermore, two novel observational signatures have recently been investigated in MOG black hole shadows perturbed by passing gravitational waves (Lobos & Rodulfo 2026). We focus on a particle orbiting a MOG black hole (BH). To do this, the space-time geometry around compact objects is described in the next paragraph.

2 Background: Modified gravity, and the MOG black hole

Following Moffat (2021), we review the theoretical basis of MOG in Sect. 2.1, which is a truncated version of the initial theory proposed by Moffat (2006). In Sect. 2.2 we briefly review the equations governing test particle motion, the fifth force, and the classical limit of the theory. We then focus on the case of a MOG black hole in Sect. 2.3, which is studied in more detail throughout the rest of this work. In the current paper, the metric signature convention (+, −, −, −) is used.

2.1 Modified gravity

In General Relativity, the gravitational field is carried by the metric tensor g_µν. In addition to this symmetric field component, MOG adds two further aspects to the gravitational interaction: a massive vector field ϕ_µ, and a scalar field χ = 1/G, where G is the coupling strength of gravity. With c = 1 units, the action is expressed as $$S=S_G+S_{\varphi }+S_M,$$(1)

where $$S_G=\frac{1}{16\pi }{{\displaystyle \int }}^{\text{}}{d}^{4}x\sqrt{-g}\left(\chi R+\frac{{\omega }_M}{\chi }{\nabla }^{\mu }\chi {\nabla }_{\mu }\chi +2\Lambda \right),$$(2a) $$S_{\varphi }={{\displaystyle \int }}^{\text{}}{d}^{4}x\sqrt{-g}\left(-\frac{1}{4}{B}^{\mu \nu }{B}_{\mu \nu }+\frac{1}{2}{\mu }_{\varphi }^2{\varphi }^{\mu }{\varphi }_{\mu }\right),$$(2b)

with g = det g_µν, R the Ricci scalar, Λ the cosmological constant, S_M the matter action, B_µν ≡ ∂_µϕ_ν − ∂_νϕ_µ an anti-symmetric Faraday-like tensor, and ω_M a constant. Moreover, G is expanded by G = G_N(1 + α), and µ_ϕ is the effective running mass of the spin 1 graviton vector field. G_N denotes Newton’s gravitational constant.

To obtain the field equations, the action with respect to g_µν, ϕ_µ and χ has to be varied, which leads to $$\begin{array}{cc}{G}_{\mu \nu }& =-\frac{{\omega }_M}{{\chi }^2}({\nabla }_{\mu }\chi {\nabla }_{\nu }\chi -\frac{1}{2}{g}_{\mu \nu }{\nabla }^{\alpha }\chi {\nabla }_{\alpha \chi })\\ & -\frac{1}{\chi }\left({\nabla }_{\mu }\chi {\nabla }_{\nu }\chi -g_{\mu \nu }\square \chi \right)+\frac{8\pi }{\chi }{T}_{\mu \nu },\end{array}$$(3a) $${\nabla }_{\nu }{B}^{\mu \nu }+{\mu }_{\varphi }^2{\varphi }^{\mu }=J^{\mu },$$(3b) $$\square \chi =\frac{8\pi }{\left(2{\omega }_M+3\right)}T,$$(3c)

respectively, where G_µν is the usual Einstein tensor, □ = ∇^µ∇_µ, J_µ is the current matter source of the vector field, and finally, the energy-momentum tensor is $$T_{\mu \nu }=T_{\mu \nu }^M+T_{\mu \nu }^{\varphi }+g_{\mu \nu }\frac{\chi \Lambda }{8\pi }.$$(4)

The $T_{\mu \nu }^{\varphi }$ term is defined as $$T_{\mu \nu }^{\varphi }=-\left(B_{\mu }{}^{\alpha }{B}_{\alpha \nu }-\frac{1}{4}{g}_{\mu \nu }{B}^{\alpha \beta }{B}_{\alpha \beta }+{\mu }_{\varphi }^2{\varphi }_{\mu }{\varphi }_{\nu }-\frac{1}{2}{g}_{\mu \nu }{\varphi }^{\alpha }{\varphi }_{\alpha }\right),$$(5)

and T is the conventional trace T = g^µνT_µν.

2.2 Test-particle motion, the fifth force, and the classical limit

The action of a test-particle with mass m is $$S_{TP}=-m{c}^2{{\displaystyle \int }}^{\text{}}d\tau -\frac{q}{c}{{\displaystyle \int }}^{\text{}}{\varphi }_{\mu }\frac{d{x}^{\mu }}{d\chi }d\tau ,$$(6)

where the second term is added due to the fifth force $(ℱℱ)$ (Moffat 2006). This additional interaction is only relevant for massive particles, where $q=\sqrt{\alpha G_{N}m}$ is the particle $ℱℱ$ charge, and ϕ_µ is the vector field. The original scalar field is the inverse gravitational coupling strength $\chi =G^{-1}={(}^1{G}_N^{-1}$, which for a MOG black hole is constant (Moffat 2021). We can therefore treat α as a dimensionless constant. In a general space-time characterised by Christoffel symbols ${\Gamma }_{\alpha \beta }^{\mu }$, the equation of motion of a test particle with four-velocity U^µ is obtained from the action (6) (Moffat 2015) $${\dot{U}}^{\mu }+{\Gamma }_{\alpha \beta }^{\mu }{U}^{\alpha }{U}^{\beta }=\{\begin{array}{ll}0\hfill & \text{massless\hspace{0.17em}test\hspace{0.17em}particle}\hfill \\ \frac{q}{m}{B}_{\nu }^{\mu }{U}^{\nu }\hfill & \text{massive\hspace{0.17em}test\hspace{0.17em}particle}\hfill \end{array}.$$(7)

Due to the $ℱℱ$ term on the right-hand side, the motion is non-geodesic in the massive case.

As derived by Moffat (2006), the acceleration of a massive particle in the weak-field regime due to a central point mass M is given by $$a_{\text{MOG}}(r)=-\frac{G_{N}M}{r^2}\left[1+\alpha -\alpha e^{-{\mu }_{\varphi }r}(1+{\mu }_{\varphi }r)\right].$$(8)

At large distances, when µ_ϕr ≫ 1, the term proportional to the exponential function vanishes, and what survives is the Newtonian gravitational acceleration, but enhanced by a factor (1 + α). Via this feature, the theory is able to describe the observed flatness of galaxy rotation curves without postulating additional mass (Moffat & Rahvar 2013; Brownstein & Moffat 2006; Davari & Rahvar 2020; Moffat & Toth 2023; Green & Moffat 2019; Moffat & Toth 2015). When µ_ϕr ≪ 1, we can linearise the exponential term, resulting in a cancellation of α, and recovering the Newtonian classical limit. This means that for small length scales (but still in the weak field), the repulsive vector field exactly cancels out the enhancement of the gravitational constant.

2.3 MOG black holes

When we consider a BH in this theory, additional simplifications can be introduced. Namely, Moffat (2021) demonstrated that starting from the field equations (3) for the matter-free ϕ_µ field-vacuum case with zero cosmological constant, we can further set the running mass to zero (i.e. µ_ϕ = 0), which causes χ to be just a constant. The latter indicates that α from now on is a parameter that quantifies the fractional increase of the MOG gravitational constant G compared to the Newtonian one G_N. Moffat & Toth (2009) also justified why the same simplifications can be admitted starting from the full general version of the theory.

The static spherically symmetric solution can then be written as (Moffat 2015) $$d{s}^2=\Psi (r)c^{2}d{t}^2-{\Psi }^{-1}(r)d{r}^2-r^{2}d{\Omega }^2,$$(9)

where dΩ² = dθ² + sin² θdϕ² and $$\Psi (r)=1-\frac{2GM}{c^{2}r}+\frac{G{Q}^2}{c^4{r}^2}.$$(10)

M is the BH mass, and $Q=\sqrt{\alpha G_{N}M}$ denotes its $ℱℱ$ charge. Setting α = 0 recovers the Schwarzschild solution in GR. In Sect. 3, we comment on the relation between the line element in Eq. (9) and the Reissner-Nordström solution for a charged BH, as well as for the motion of test particles in both theories and space-times.

Moreover, the vector field for a MOG black hole with mass M in Schwarzschild-like coordinates is given by (Moffat 2021; Della Monica et al. 2021) $${\varphi }_{\mu }=\left(\frac{\sqrt{\alpha G_{N}M}}{r},0,0,0\right),$$(11)

where the choice of the sign is explained in Appendix A.

3 Relation to Reissner–Nordström space-time

In GR, the space-time of a BH with mass $\tilde{M}$ and electric charge $\tilde{Q}$ is given by the Reissner–Nordström (RN) solution (Reissner 1916; Nordström 1918). The line element in standard spherical coordinates and in Gaussian units (so that 1/(4πε₀) = 1) has the same form as Eq. (9), with $$\Psi (r)\mapsto {\Psi }_{\text{RN}}(r)=1-\frac{2G_N\tilde{M}}{c^{2}r}+\frac{G_N{\tilde{Q}}^2}{c^4{r}^2}.$$(12)

This relation with a MOG BH was pointed out (Moffat & Toth 2009; Moffat 2015) under the straightforward relabelling $$G\mapsto G_N\text{and}Q\mapsto \tilde{Q}.$$(13)

Here, we present a different substitution, which clearly shows that not only are there striking similarities, but the two theories produce indistinguishable geometries and equations of motions under a certain fine-tuning of their parameters. In particular, writing Eq. (10) as $${\Psi }_{\text{MOG}}(r)=1-\frac{2G_N[(1+\alpha )M]}{c^{2}r}+\frac{G_N{\left[\sqrt{1+\alpha }Q\right]}^2}{c^4{r}^2},$$(14)

and identifying $$\tilde{M}=(1+\alpha )M\text{and}\tilde{Q}=\sqrt{1+\alpha }Q,$$(15)

leads to $${\Psi }_{\text{MOG}}(r;M,Q)={\Psi }_{\text{RN}}(r;\tilde{M},\tilde{Q}).$$(16)

In other words, under the fine-tuning of Eq. (15), the line elements of the two theories are indistinguishable, and hence, so are the left hand sides (geodesic part) of their equations of motion (expressed by Eq. (7)). With the above straightforward dictionary of Eq. (13), this equivalence is not achieved because the value of α affects G and Q both.

Furthermore, the right-hand sides of the equations of motion (Eq. (7)) for massive test particles become indistinguishable for the two theories when the fine-tuning of Eq. (15) is also applied to the test particle, $$\tilde{m}=(1+\alpha )m\text{and}\tilde{q}=1+\alpha q.$$(17)

In the RN case, this right-hand side term for a massive $(\tilde{m})$ and electrically charged $(\tilde{q})$ test particle is $${\text{RHS}}_{\text{RN}}=\frac{\tilde{q}}{\tilde{m}}{g}^{\mu \rho }({\partial }_{\rho }{A}_{\nu }-{\partial }_{\nu }{A}_{\rho })U^{\nu },\text{with}{A}_0=\frac{\tilde{Q}}{r}$$(18)

being the only non-zero component of the four potential. In the MOG case, the right-hand side for a test particle with mass m and $ℱℱ$ charge q is $${\text{RHS}}_{\text{MOG}}=\frac{q}{m}{g}^{\mu \rho }({\partial }_{\rho }{\varphi }_{\nu }-{\partial }_{\nu }{\varphi }_{\rho })U^{\nu },\text{with}{\varphi }_0=\frac{Q}{r}$$(19)

being the only non-zero component of the four potential. Imposing the fine-tuning of Eqs. (15) and (18) then yields the equivalence of the two right-hand sides.

We summarise some of the observations and implications of the above discussion below.

• Under the fine-tuning of Eq. (15), geodesics around a MOG BH with mass and fifth-force charge (M, Q) are indistinguishable from the geodesics around a GR BH $(\tilde{M},\tilde{Q})$. Hence, massless test particles follow the same dynamics around a BH in both universes.

• Massive but electrically neutral test particles also follow the geodesics of the Reissner–Nordström space–time, while massive test particles around a MOG BH deviate from these geodesics due to the proportionality of the fifth-force charge q to the mass. In other words, there is no MOG analogue to the electrically neutral massive test particle of GR.

• The additional fine-tuning of Eq. (18) also establishes an equivalence between the (non-geodesic) dynamics of massive particles around a MOG BH (M, Q) and a charged GR BH $(\tilde{M},\tilde{Q})$.

It should be emphasised, however, that the MOG theory around a BH corresponds to a reduced version of the full theory, as discussed in Sect. 2.3. Consequently, our fine-tuning procedure may not remain directly applicable in other systems, particularly when α is not constant.

After establishing the above correspondence between the two theories, we demonstrate its utility through a brief discussion of the horizons of a MOG BH. The Reissner–Nordström BH horizons are given by $$r_{\pm }=\frac{G_N\tilde{M}}{c^2}\pm \frac{\sqrt{G_N^2{\tilde{M}}^2-G_N{\tilde{Q}}^2}}{c^2},$$(20)

and since any $ℱℱ$ charge is given by $Q=\sqrt{\alpha G_{N}M}$, starting from (15), for a positive α, we have $$\frac{G_N\tilde{M}}{\sqrt{G_N\tilde{Q}}}=\sqrt{\frac{1+\alpha }{\alpha }}>1,$$(21)

meaning that a MOG BH always has two horizons (Moffat 2015), and it can only approach the extremal case when α → ∞.

In conclusion, a MOG BH is a very specific case of the RN solution. The latter is free to take any mass $\tilde{M}$ and electric charge $\tilde{Q}$ (the same holds for every other massive object), meaning that on the mass-charge plane, every point is allowed. In contrast, a MOG BH (and again, every other massive object) will live on the straight line $\tilde{Q}=\sqrt{\alpha G_N/(1+\alpha )}\tilde{M}$ on this plane. Thus, by observing trajectories around BHs, we would not be able to infer whether our universe follows MOG or Einstein-Maxwell under the fine-tuning (15), but now for every mass and charge. An alternative interpretation might be that a dark sector is associated with dark photons and charges, leading to a modified RN space-time (Morris 2023). A compact object then has a total mass M_tot = M_D + M_S and a total charge Q_tot = Q_D + Q_S, where the D and S indices refer to the dark and standard model charge, respectively. In this context, we can set Q_S = 0 and identify M_S = M, M_D = αM and $Q_D=\sqrt{1+\alpha }Q$, such that $M_{\text{tot}}=\tilde{M}$ and $Q_{\text{tot}}=\tilde{Q}$, leading to indistinguishable dynamics once again.

4 Post-Newtonian dynamics of a massive test particles orbiting a MOG black hole

4.1 Derivation of the equation of motion

From now on, we recover the powers of c, such that the massive test-particle action of Eq. (22) is $$S_{TP}=-m{c}^2{{\displaystyle \int }}^{\text{}}d\tau -\frac{q}{c}{{\displaystyle \int }}^{\text{}}{\varphi }_{\mu }\frac{d{x}^{\mu }}{d\tau }d\tau .$$(22)

The relative normalisation of powers of c for the two components, as well as the relative signs, are consistent with a charged particle in an electromagnetic field (Jackson 1998). For the test particle orbiting a MOG black hole with mass M (see Sect. 2.3), the action S_TP in Schwarzschild-like coordinates can be written as $$S_{TP}={{\displaystyle \int }}^{\text{}}dt{L}_{TP}=-{{\displaystyle \int }}^{\text{}}dtm{c}^2\left[\frac{d\tau }{dt}+\alpha \frac{G_{N}M}{r{c}^2}\right],$$(23)

where L_TP denotes the test-particle Lagrangian. Our goal is to derive the equations of motion up to first post-Newtonian (1PN) order from this action. Following the method of Gainutdinov (2020), we achieved this by first transforming to isotropic coordinates, then developing L_TP to the desired powers of c⁻², and finally, transforming to Cartesian-like coordinates. The full steps can be found in Appendix B. The Lagrangian, obtained by combining Eqs. (B.7) and (B.8), is $$\begin{array}{ll}\frac{L_{TP}}{m}=\hfill & \frac{{\dot{x}}^2}{2}\left[1+\frac{{\dot{x}}^2}{4c^2}-3(\alpha ℳ+1)\frac{{\phi }_N}{c^2}\right]\hfill \\ \hfill & -(\alpha ℳ+1){\phi }_N\left[1+(2\alpha ℳ+1)\frac{{\phi }_N}{2c^2}\right]\hfill \\ \hfill & +Q\left\{\alpha {\phi }_N+\alpha (1+\alpha )\frac{{\phi }_N^2}{c^2}\right\}+\mathcal{O}(c^{-4}),\hfill \end{array}$$(24)

where x denotes the spatial position vector in Cartesian-like coordinates, and φ_N = −G_N M/|x| denotes the Newtonian central mass potential. The two parameters ℳ, 𝒬 were inserted by hand as book-keeping of the origin of the various terms. In the full theory, they have a numerical value equal to 1. However, in Sect. 5, we occasionally also set them to zero in order to investigate each contribution in isolation. Terms not proportional to either bookkeeping parameter are due to the Schwarzschild metric, and are thus also present in the 1PN expansion of GR. These are indeed the only surviving terms for α = 0, as expected. Terms proportional to ℳ are due to the deviation of the MOG BH metric from the Schwarzschild metric when α > 0. Terms proportional to 𝒬 are due to the fifth force.

Finally, we plugged the Lagrangian of Eq. (24) into the Euler–Lagrange equations to obtain the MOG equation of motion up to 1PN order (see also Gainutdinov 2020 for a related but less explicit result), $$\begin{array}{llll}\ddot{x}\hfill & =\hfill & -\nabla {\phi }_N-\alpha {\cancel{\left(ℳ-𝒬\right)}}^0\nabla {\phi }_N\hfill & \left(\text{Kepler}\right)\hfill \\ \hfill & \hfill & -\nabla {\phi }_N\left(4\frac{{\phi }_N}{c^2}+\frac{{\ddot{x}}^2}{c^2}\right)+4\left(\nabla {\phi }_N\cdot \frac{\dot{x}}{c}\right)\frac{\dot{x}}{c}\hfill & \left(\text{GR\hspace{0.17em}1\hspace{0.17em}PN}\right)\hfill \\ \hfill & \hfill & +\alpha \left(4ℳ-𝒬\right)\left(\nabla {\phi }_N\cdot \frac{\dot{x}}{c}\right)\frac{\dot{x}}{c}\hfill & \left(\text{MOG\hspace{0.17em}1PN}\right)\hfill \\ \hfill & \hfill & -\alpha \nabla {\phi }_N\left[\left(2ℳ+𝒬\right)\frac{{\dot{x}}^2}{2c^2}\right]\hfill & \left(\text{MOG\hspace{0.17em}1PN}\right)\hfill \\ \hfill & \hfill & -\alpha \nabla {\phi }_N\left[\left\{5\left(\alpha +1\right){\cancel{\left(ℳ-𝒬\right)}}^0+4ℳ\right\}\frac{{\phi }_N}{c^2}\right].\hfill & \left(\text{MOG\hspace{0.17em}1PN}\right)\hfill \end{array}$$(25)

We ordered the terms by powers in 1/c⁻² and marked the physical interpretation of each line on the right. Moreover, recalling that the parameters ℳ and 𝒬 are numerically equal to 1, we stress the fact that (ℳ − 𝒬) is in reality 0, but we included it to recall the book-keeping order of the expression.

This derivation can be directly generalised to any modified theory of gravity, for instance, Yukawa gravity (Tan & Lu 2024). We conclude this section with a few key observations.

• There is no MOG correction to the classical (Kepler) limit because the respective individual contributions of metric and fifth force cancel each other out. This is consistent with Eq. (8) in the limit µ_ϕr ≫ 1.

• Although the fifth force also counters the metric effect on the 1PN level to some extent, there are surviving terms. Therefore, although MOG is constructed to account for the missing mass problem via its weak-field limit Eq. (8), its deviation from GR is also apparent in the strong-field regime.

• For α = 0, we recover the GR equations of motion up to 1PN order (Gainutdinov 2020), as it should be.

• We recall that Eq. (25) constitutes an expression in isotropic coordinates (modulo the transformation to the Cartesian-like spatial basis).

4.2 Formulation in terms of osculating Kepler orbits

Motivated by the final point of the previous section, we prepared Eq. (25) for an investigation into the effect of MOG on an orbit that probes the PN regime of a central mass M. To do this, we viewed our equations of motion as a perturbed Kepler problem, $$\ddot{x}=-\frac{G_{N}M}{{\rho }^2}n+a_{\text{p}},$$(26)

where the perturbative acceleration a_p abbreviates all terms on the right-hand side of Eq. (25) but the first one, ρ = |x| and n = x/ρ. We still kept track of the terms proportional to $(ℳ-\mathcal{Q})$ for later utility. We projected a_p onto the Gaussian frame (Poisson & Will 2014; Merritt 2013), $$a_{\text{p}}=ℛn+S\lambda +\mathcal{W}{e}_{\text{z}},$$(27)

where e_z is the unit vector aligned with the orbital angular momentum, and λ completes the right-handed orthonormal basis (Fig. 1). The components can be further expressed as $$ℛ={ℛ}_{\text{GR\hspace{0.17em}1PN}}+{ℛ}_{\text{MOG\hspace{0.17em}1PN}},$$(28a) $$\mathcal{S}={\mathcal{S}}_{\text{GR\hspace{0.17em}1PN}}+{\mathcal{S}}_{\text{MOG\hspace{0.17em}1PN}},$$(28b) $$\mathcal{W}={\mathcal{W}}_{\text{GR\hspace{0.17em}1PN}}+{\mathcal{W}}_{\text{MOG\hspace{0.17em}1PN}},$$(28c)

where the GR 1PN components read (Merritt 2013) $${ℛ}_{\text{GR\hspace{0.17em}1PN}}=\frac{G_{N}M}{c^2{\rho }^2}\left[\frac{4G_{N}M}{\rho }-{\dot{x}}^2+4{\left(n\cdot \dot{x}\right)}^2\right],$$(29a) $${\mathcal{S}}_{\text{GR\hspace{0.17em}1PN}}=\frac{4G_{N}M}{c^2{\rho }^2}\left(n\cdot \dot{x}\right)\left(\lambda \cdot \dot{x}\right),$$(29b) $${\mathcal{W}}_{\text{GR\hspace{0.17em}1PN}}=0,$$(29c)

with $$\rho =\frac{p}{1+e\cos f}=\frac{a(1-e^2)}{1+e\cos f},$$(30a) $${\dot{x}}^2=\frac{G_{N}M}{p}(1+e^2+2e\cos f),$$(30b) $$\left(n\cdot \dot{x}\right)=\sqrt{\frac{G_{N}M}{p}}e\sin f,$$(30c) $$\left(\lambda \cdot \dot{x}\right)=\sqrt{\frac{G_{N}M}{p}}(1+e\cos f).$$(30d)

For the MOG components, we derived $$\begin{array}{lll}{\Re }_{\text{MOG}}\hfill & =\hfill & -\alpha \frac{G_Nℳ}{{\rho }^2}{\cancel{\left(ℳ-𝒬\right)}}^0\hfill \\ \hfill & \hfill & +\alpha \frac{G_Nℳ}{{\rho }^2}\{\left(4ℳ-𝒬\right)\frac{{\left(n\cdot \dot{x}\right)}^2}{c^2}-\left(2ℳ+𝒬\right)\frac{{\dot{x}}^2}{2c^2}\hfill \\ \hfill & \hfill & +\frac{G_Nℳ}{c^2\rho }\left[5\left(\alpha +1\right){\cancel{\left(ℳ-𝒬\right)}}^0+4ℳ\right]\},\hfill \end{array}$$(31a) $${\mathcal{S}}_{\text{MOG}}=\alpha (4ℳ-Q)\frac{G_{N}M}{c^2{\rho }^2}\left(n\cdot \dot{x}\right)\left(\lambda \cdot \dot{x}\right),$$(31b) $${\mathcal{W}}_{\text{MOG\hspace{0.17em}1PN}}=0.$$(31c)

Following Poisson & Will (2014), we used the formalism of osculating orbits, in which the perturbed Kepler problem (given by Eq. (26)) is reformulated as a set of evolution equations for the six orbital elements (p, e, ι, Ω, ω, and f) correspondent to semi-latus rectum, eccentricity, inclination, argument of the ascending node, argument of pericentre, and true anomaly, respectively (see Fig. 1). The full equations are explicitly reported in Eq. (3.64) of Poisson & Will (2014), while they are here written for the special case 𝒲 = 0, $$\dot{p}=\sqrt{\frac{p^3}{G_{N}M}}\frac{2}{1+e\cos f}\mathcal{S},$$(32a) $$\dot{e}=\sqrt{\frac{p}{G_{N}M}}\left[\sin fℛ+\frac{2\cos f+e(1+{\cos }^{2}f)}{1+e\cos f}\mathcal{S}\right],$$(32b) $$\dot{\omega }=\frac{1}{e}\sqrt{\frac{p}{G_{N}M}}\left[-\cos fℛ+\frac{2+e\cos f}{1+e\cos f}\sin f\mathcal{S}\right],$$(32c)

while $\dot{\iota }=0$ and $\dot{\Omega }=0$. For our problem, the remaining Gaussian frame components ℛ and 𝒮 are given by Eqs. (28), (29) and (31). Finally, we obtained closed-form expressions by using the relations (30).

At any point in time, the transformation between the orbital elements, position, and velocity is given by Eqs. (C.1). Some qualitative features of the MOG orbital mechanics to 1PN order can be seen directly from Eqs. (32):

• Because of the vanishing out-of-plane component of the per-turbative acceleration 𝒲, the elements ι and Ω are constants of motion, such that MOG does not alter the orbital plane.

• Because of the non-vanishing in-plane components of the perturbative acceleration ℛ_MOG, 𝒮_MOG (Eqs. (31)), MOG contributes to the evolution of p, e, and ω. The latter implies in particular that it contributes to the rotation of the orbit within its plane.

Fig. 1 Definition of the Euler angles (Ω, ι, ω) giving the orbital orientation and of the true anomaly f giving the current position within the orbit. The static fundamental frame (X, Y, Z) (blue) is adapted to the observables, while the co-rotated Gaussian frame (n,λ, ez) (orange) is adapted to the symmetry of the (perturbed) two-body problem (reproduced from Heißel, et al. 2022).

4.3 Secular evolution

Next, we investigated whether the evolutions of p, e, and ω exhibited secular components, that is, whether there are non-vanishing changes in these elements per revolution, for instance, ∆ω = ω(f + 2π) − ω(f), and hence, accumulating effects over several orbits. To do this, we assumed the trajectories to be close to Keplerian, which is valid on the distance scales in which we are interested. Then, we considered the first-order approximation of the evolution Eq. (32), $$\frac{dp}{df}\simeq \frac{p^3}{G_{N}M}\frac{2}{{\left(1+e\cos f\right)}^3}\mathcal{S}$$(33a) $$\frac{de}{df}\simeq \frac{p^2}{G_{N}M}\left[\frac{\sin f}{{\left(1+e\cos f\right)}^2}ℛ+\frac{2\cos f+e(1+{\cos }^{2}f)}{{\left(1+e\cos f\right)}^3}\mathcal{S}\right]$$(33b) $$\frac{d\omega }{df}\simeq \frac{p^2}{G_{N}M}\frac{1}{e}\left[-\frac{\cos f}{{\left(1+e\cos f\right)}^2}ℛ+\frac{2+e\cos f}{{\left(1+e\cos f\right)}^3}\sin f\mathcal{S}\right],$$(33c)

in which the orbital elements on the left-hand side are viewed as functions of f , all orbital elements on the right-hand side except f are assumed constant, and $\dot{f}$˙ is approximated by its non-perturbed Keplerian expression (Poisson & Will 2014, p. 160).

Plugging in ℛ and 𝒮 from Eqs. (28)–(31), we then integrated Eq. (33) analytically over one period. More precisely, the secular change of an orbital element µ^a over a complete orbit is given by $$\Delta {\mu }^a={\displaystyle {\int }_0^P\frac{d{\mu }^a}{dt}dt}={\displaystyle {\int }_0^{2\pi }\frac{d{\mu }^a}{df}df,}$$(34)

where P is the orbital period. Thus, we found ∆p = 0, ∆e = 0, and $$\Delta \omega =\frac{6\pi G_{N}M}{p{c}^2}\left(1+\frac{5}{6}\alpha \right).$$(35)

Moreover, for an orbital element ${\mu }_i^a$, denoting a combination of effects i (as $i\in \{\text{GR},\mathcal{S}\mathcal{T},ℱℱ,\dots \}$, see Table 1), we defined $\delta {\mu }_i^a$ as $$\delta {\mu }_i^a(f)\equiv \Delta {\mu }_i^a(f)-\Delta {\mu }_{\text{GR}}^a(f).$$(36)

This is useful in our discussion of Fig. 5 below.

This shows that

• there is no secular evolution of p and e, such that these elements return to their original values after each revolution,

• the only element other than f that exhibits a secular evolution is ω, and hence, there is a precession of the orbit within its plane. For α = 0, we recover the GR expression, famous ever since its first derivation by Einstein (1915). For α > 0, MOG amplifies this precession by the above factor.

• We note, however, that MOG and GR can both very well yield the same ∆ω since the gain in precession by an α > 0 can be counterbalanced by a lower central mass. Through correlations like this, that knowledge of the secular behavior alone is insufficient to distinguish between different theories of gravity and dark matter. To do this, an understanding of the distinct non-secular signatures of each model is key (see Sect. 5 and Heißel, et al. (2022).

• Our result (Eq. (35)) confirms the result of Della Monica et al. (2021, 2023), who performed an alternative derivation. In Appendix. A, we also settle a disagreement between these authors and Turimov (2022).

4.4 Observational effects: Transverse Doppler shift and gravitational redshift

When the interest is not only in the dynamics of test particles around a MOG black hole as such, but furthermore, in how this dynamics would be observed from a distance, MOG effects on the light by which we observe have to be considered as well. We restricted our considerations here to two such effects, namely the transverse Doppler shift and the gravitational redshift, since these are currently the only relativistic effects of relevance to the observed light of star S2 in the Galactic centre. Moreover, we consider this latter system as the example to investigate MOG features in the non-secular evolution in Sect. 5.

First, we note that as a special relativistic effect, the transverse Doppler effect is unaltered by MOG. Furthermore, since photons are massless, there is no fifth-force contribution to the gravitational redshift. The MOG BH metric (see Eqs. (9)–(10)), however, does yield modifications to the well-known GR red-shift. Labelling the event at which a photon is emitted at the star by P and the event at which this photon is detected on Earth by Q, we conclude according to Misner et al. (1973) that for the gravitational redshift z, we have $$1+z=\frac{{\nu }_P}{{\nu }_Q}=\sqrt{\frac{g_{00}{|}_Q}{g_{00}{|}_P}}=\sqrt{\frac{\Psi (|x_Q|)}{\Psi (|x_P|)}}=\sqrt{\frac{\Psi (R)}{\Psi (\rho )}},$$(37)

where ν_P, ν_Q are the emitted and observed frequencies, respectively, R is the distance from Earth to the Galactic centre, and Ψ is given by Eq. (10) transformed into isotropic coordinates (see Eq. (B.5)). On the other hand, the radial (line-of-sight) velocity, namely RV, is inferred from redshift via $$1+z=\sqrt{\frac{1+\Delta RV/c}{1-\Delta RV/c}}\approx 1+\frac{\Delta RV}{c}$$

for |∆RV| ≪ c , by which we obtain $$\frac{\Delta RV}{c}\approx \sqrt{\frac{1-2(\alpha +1)\frac{G_{N}M}{R_0{c}^2}}{1-2(\alpha +1)\frac{G_{N}M}{\rho c^2}}}-1$$(38)

due to the gravitational redshift, where we kept Ψ up to 1PN order. Therefore, in order to obtain the radial velocity that the star had at emission, we need to subtract the ∆RV (given by Eq. (38)) from the RV_Q measured on Earth, obtaining RV_P = RV_Q − ∆RV. For α = 0, Eq. (38) reduces to the familiar result of GR.

5 Effects of MOG on stellar orbits in the Galactic centre

5.1 Model setup: Star S2 in orbit around Sagittarius A*

For the system in the example of which we conducted our numerical investigations, we chose star S2, which is in orbit around the massive black hole Sgr A* in the Galactic centre (Ghez et al. 2008; Gillessen et al. 2009; GRAVITY Collaboration 2022). With a pericentre distance of about 1400 Schwarzschild radii, it probes the post-Newtonian regime of Sgr A*. Coupled with a high eccentricity of about 0.88, its orbit serves as an excellent laboratory for testing relativity (GRAVITY Collaboration 2018, 2020; Do et al. 2019; Saida et al. 2019), dark matter (Heißel, et al. 2022; Chan et al. 2022; Della Monica & de Martino 2023; GRAVITY Collaboration 2022), and alternative theories of gravity (Borka et al. 2013; De Martino et al. 2021; Della Monica et al. 2022a,b; Abd El Dayem et al. 2025).

The hypothetical scenario in which Sgr A* is assumed to be a MOG black hole was investigated by (Della Monica et al. 2021), who obtained the upper bound α < 0.662 with a confidence of 99.7%. The measurement was based on fitting the fully relativistic equations of motion (Eq. (7)) to publicly available data of astrometric and spectroscopic S2 observations. We instead investigated which distinct signature a MOG black hole inscribes on the S2 orbit as compared to a Schwarzschild black hole, and we compared the astrometric and spectroscopic observables, as well as the evolution of the orbital elements. We hence based our model on the post-Newtonian osculating-orbit framework discussed in Sect. 4.2. This choice also facilitates comparisons to dark matter scenarios (see Sect. 5.4) and prepares MOG for potential future Solar System tests, in which Newtonian dominate relativistic perturbations, such that a fully relativistic model is practically impossible to set up.

In the following, model orbits are based on the integration of the osculating Eqs. (32) with ℛ, 𝒮 given by Eqs. (28)–(31). For the initial point (t₀), we chose the (at the time of writing) last apocentre passage of the star in 2010. For the other initial orbital elements, the black hole mass and the distance to the Galactic centre, we chose the same values as in GRAVITY Collaboration (2022, Table B.1). In summary, our model parameters common to all orbits were $$p_0=0.00109214\text{\hspace{0.17em}pc},e_0=0.88441,{\iota }_0={134.70}^{\circ },$$(39a) $${\omega }_0={66.25}^{\circ },{\Omega }_0={228.19}^{\circ },f_0=\pi ,$$(39b) $$M=4.297\times {10}^6{M}_{\odot },R=8277\text{\hspace{0.17em}pc},$$(39c)

where the orbital elements correspond to the initial time t₀. The sole purpose of the distance to the Galactic centre R is to translate between physical distances and angular distances on sky. In one investigation, which we point out explicitly, we chose different values for the black hole mass. The parameters that we varied to investigate the theory were $$\alpha \ge 0,ℳ\in \{0,1\}\text{\hspace{0.17em}and\hspace{0.17em}}Q\in \{0,1\},$$(40)

where the latter two were set to 1, unless stated otherwise, and were occasionally set to 0 individually in order to investigate the fifth force from the MOG space-time contribution in isolation.

Since we did not fit observations, but focused on features of MOG, we did not model systematic effects such as the motion of the Solar System or the offset and drift of instrument reference frames (as in e.g. GRAVITY Collaboration (2020) and Reid & Brunthaler (2020)). Moreover, in comparison plots, where we subtract a quantity of a GR orbit from the corresponding quantity of a MOG orbit, these contributions would cancel out in any case because they affect the two orbits identically. To investigate the effects on the observables (RA, DEC, and RV), which are the right ascension, declination, and radial velocity, we transformed the orbital elements in position and velocity using Eqs. (C.1), such that (RA, Dec, RV) = (r^Y, r^X, −v^Z) (Fig. 1).

5.2 Comparing MOG orbits to GR orbits

In Sect. 4.3, we determined that for the same central mass M, MOG enhances the secular pericentre advance of GR (Schwarzschild precession) by a factor 1 + 5α/6 (Eq. (35)). Figure 2 demonstrates how this alters the appearance of the observed orbit on sky for the example value α = 5, yielding an almost sixfold Schwarzschild precession. Since publicly available observational data readily constrain α < α_*:= 0.662 (Della Monica et al. 2021), this is an extreme choice for the sole demonstrative purpose in which the effect is visible by naked eye.

To place into perspective what is actually observable, we note that the current 1σ measurement uncertainties for S2 observations are about 50 µas in astrometry using the GRAVITY instrument at the Very Large Telescope Interferometer (GRAVITY Collaboration 2017) and 10 km/s in radial velocity using the Spectrograph for Integral Field Observations in the Near Infrared (Eisenhauer et al. 2005) and the Enhanced Resolution Imager and Spectrograph (Davies et al. 2023) at the Very Large Telescope. Therefore, in order to probe the detection threshold for the deviation of a MOG orbit from a GR orbit for realistic α values, we present Fig. 3. We subtracted the observables of the GR orbit from the observables of the MOG orbit on a logarithmic scale as functions of time for three values of α, the highest of which is the bound α_* found by Della Monica et al. (2021). By our choice of initial parameters (Eq. (39)), the orbits coincide at the 2010 apocentre passage and then deviate with a certain signature within each 16-year period and an overall accumulation of the deviation per revolution according to Eq. (35).

The precise form of the non-secular deviation has to be interpreted with some caution. In Fig. 3, which shows the observables in the time domain, the sharp features around pericentre are in part due to the fact that the orbital period changes with α. This ought to be automatically corrected in the true anomaly domain. A certain f value might correspond to different times t for orbits in different gravity theories. Furthermore, due to the high eccentricity we considered (see Eq. (39)), the part of the orbit around the pericentre is traversed in a rather short window in the time domain, which is stretched out slightly further in the true anomaly domain (see the grey area in Fig. 6).

Overall, we observe that for α equal to the current bound (right column of Fig. 3), the MOG deviation from GR (solid blue lines) surpasses the 1σ instrument threshold (dashed horizontal line) after a single revolution (≈16 years), and the 3σ limit (horizontal black line) shortly after. For lower α values (left and middle columns of Fig. 3), the effect has to accumulate for longer in order to exceed the thresholds. However, an improvement of the (Della Monica et al. 2021) limit by a factor 1/10, namely α = 6.62 · 10⁻², represented in the middle column of Fig. 3, seems achievable within realistic timescales.

Furthermore, the RA deviation has a positive offset with respect to the DEC deviation because of the spatial orientation of the orbit. Finally, we also note that the RV deviation needs to accumulate longest in order to breach the thresholds. However, we recall that the deviation plotted here is solely due to the effects of MOG on the orbital motion itself. In addition to this, spectral observations are subject to the gravitational redshift, which is altered by MOG in comparison to GR as well (Sect. 4.4). We discuss the MOG redshift contribution to RV in isolation in Sect. 6.

We finally focused on the representation of the orbits in terms of orbital elements (p, e, ω, a) and show their change over one revolution in Fig. 4. Firstly, we note that the secular evolution agrees with what we derived in Sect. 4.3, such that ∆p = ∆e = 0, and ω changes from f = 0 to 2π according to Eq. (35). Secondly, the coincidence of all orbits at f = π is again due to our arbitrary choice of initial parameters (Eq. (39)). The absolute deviation between the MOG and GR orbits thus must not be interpreted as an indicator of where (in terms of true anomaly) MOG impacts the orbit. Instead, the slope of these curves marks where deviations from Keplerian motion are driven by both GR and MOG.

Fig. 2 Oribt of star S2 (yellow star) around the BH (black dot) in the centre of our Galaxy for GR (left panel) and an extreme MOG case (right panel).

Fig. 3 Evolution of the observables’ (RA, Dec, and RV in the first, middle, and bottom row, respectively) absolute difference between GR and i) full MOG (solid blue lines), ii) fifth force (solid red lines), and iii) space-time (dashed yellow lines) contribution for different values of α. We fixed α as α = 6.62 · 10−3 in the left column, α = 6.62 · 10−2 in the middle one and α = 6.62 10−1 in the right column. An explanation of the labels can be found in Table 1. The horizontal dashed lines correspond to the 1σ (black) and 3σ ·(grey) accuracy of the instruments. For the spatial observables (RA and Dec), these are 50 µas and 150 µas, respectively, and for the radial velocity (RV), they were fixed to 10 km/s and 30 km/s.

5.3 Fifth force and metric effects in isolation

As was already apparent from the fully relativistic equations of motion (Eqs. (7)), MOG accelerates massive test particles via two components: through its distinct metric from GR, and due to the fifth force, by which the motion is non-geodesic. In Sect. 4 we kept track of these separate contributions throughout our derivations of the PN equations of motion (Eq. (25)) and of the perturbative acceleration components ℛ, 𝒮, and 𝒲 (Eq. (28)) via the book-keeping parameters ℳ, 𝒬. We used these then to investigate the individual contributions to the dynamics in isolation by setting them to zero individually. We stress that this is merely a diagnostic tool and does not reflect a physical scenario in the MOG theory, since any massive object carries a fifth-force charge. Therefore, the only physical scenario is the full MOG, or equivalently, 𝒮𝒯 & $ℱℱ((ℳ,\mathcal{Q})=(1,1))$, and GR ((ℳ, 𝒬) = (0, 0)). However, due to the relation to Reissner–Nordström, shown in Sect. 3, the 𝒮𝒯-only scenario around a MOG BH can be considered equivalent to a massive but electrically neutral particle following a geodesic around an electrically charged BH in GR. In order to select an α value for this investigation, we recall that for (ℳ, 𝒬) = (1, 0) or (0, 1), the weak-field acceleration of MOG does not cancel, such that the Keplerian two-body acceleration is modified by a factor 1 ± α (first lines in Eqs. (25) and (31)). Being of 0PN order, this term then completely dominates the perturbative accelerations, such that we had to choose a very small α to yield curves on reasonable scales. We chose α = 10⁻³. For the same reason, we waited for the isolated effects to appear equal and opposite in our plot, up to small 1PN differences not visible by naked eye on these scales for such a small α. Consequently, the difference between the curves for GR (dashed black) and the full MOG (solid blue) theory are also expected to be insignificant and overlap by naked eye. All of the above can indeed be observed in Fig. 5.

The semi-latus rectum p (first column of the figure) does not strongly deviate from GR even for single contributions. The maximum effect is about 10⁻³ AU on the pericentre. In general, the space-time geometry mainly tends to slightly reduce p, but drastically increases the eccentricity away from the apocentre, and thus, the semi-major axis as well (see the second and last columns, respectively). In contrast, the repulsive force $(ℱℱ)$ has the opposite effect and cancels out 𝒮𝒯 almost totally.

Regarding the argument of pericentre (third column), the 𝒮𝒯 produces a characteristic feature of a negative slope around the apocentre, which is compensated for by the opposite feature produced by $ℱℱ$. Near pericentre, the effects compete again, leading in the overall compliance with GR. Eq. (35) indicates that ∆ω_MOG ≃ ∆ω_GR for this value of α, which explains why we do not see a discrepancy between the solid blue and dashed black curve by eye. When we consider the findings from this figure, it becomes clear that the origin of the deviation between MOG and GR in Fig. 4 is the overcompensation of $ℱℱ$ for the 𝒮𝒯 effect around the pericentre for high α values, while around the apocentre (grey shaded region), the contributions cancel out almost totally.

Fig. 4 Orbital element evolution with true anomaly. Here, ∆p stands for the offset of p from its GR value at f = 0, and analogously for the other elements. The dashed black line represents GR, and the solid coloured lines represent a MOG case with a particular α value in the range (0.001, 1). For each orbital element µa, the definition of ∆µa is given in Eq. (34).

5.4 MOG versus DM

Considering the inspiration behind MOG, the interesting question arises whether the effect of a dark component on the stellar orbits is distinguishable from that of MOG. Heißel, et al. (2022) showed that a DM profile demonstrates a characteristic feature on the evolution of the argument of pericentre, mainly in the range f ∈ (5π/6, 7π/6), and that the overall precession ∆ω_DM also decreases. We might therefore relax M as a free variable and impose ∆ω_DM = ∆ω_GR for a particular DM profile. Doing so for the MOG case as well, by relaxing α and M and imposing ∆ω_MOG = ∆ω_GR, would lead to a direct comparison between the DM and MOG signature within an orbit.

We refer to four cases, presented in Fig. 6: (i) full MOG (𝒮𝒯 and $ℱℱ$), (ii) 𝒮𝒯 contribution alone, (iii) GR plus an extended mass (DM,) and (iv) GR only. All the relativistic effects were taken into account up to 1PN. In the GR case, we selected M = 4.261 × 10⁶ M_⊙, leading to the evolution of ω(f) represented with a dashed black line (Fig. 6). Regarding MOG, we first took the full theory (𝒮𝒯 & $ℱℱ$) for α = α_* = 0.662 (Della Monica et al. 2021) and set M = 2.205 M_⊙, such that the overall precession matched (blue line).

Regarding the extended mass, we considered the Plummer profile (Plummer 1911), $$\rho (r)={\rho }_0{\left(1+\frac{r^2}{r_0^2}\right)}^{-5/2},$$(41)

where ρ₀ is the density parameter, and and r₀ = 0.3″ is the length-scale parameter (Mouawad et al. 2005). Following the most recent upper bounds on the extended mass (GRAVITY Collaboration 2022), we chose the density parameter such that the total mass up to a radial distance of 230mas was 3500 M_⊙, corresponding roughly to 1σ. To match the orbital precession, we selected M = 5.43 × 10⁶ M_⊙, leading to the red line in Fig. 6. Finally, we isolated the 𝒮𝒯 effect (yellow line), and fixed α = 0.758 × 10⁻³ and M = 4.252 × 10⁶ M_⊙, such that the slope at f = π and the overall precession matched. All the other parameters were kept the same for all four cases.

The full MOG theory, even for this high α = 0.662 (which indicates an increase of 66.2% in the gravitational strength), is almost identical to GR for almost the whole orbit. On the other hand, the DM case generates a characteristic feature around the apocentre, which was extensively discussed in Heißel, et al. (2022). It is significant that the 𝒮𝒯 effect achieves the same feature, although it deviates strongly far from the apocentre. In MOG, it is impossible for a massive test particle to follow a geodesic because of its $ℱℱ$ charge q_TP ~ m_TP. However, when we take the relation between MOG BH and the Reissner–Nordström space-time (Sect. 3) into account, it is evident that an electrically neutral test particle around a charged black hole in GR follows geodesics on a space–time that produces exactly the same feature. At very far distances compared to the BH scales, the $ℱℱ$ weakens exponentially because the mass of the vector field is non-zero, and the 𝒮𝒯 effect therefore accounts for the missing mass. Therefore, these feature similarities should be expected.

Fig. 5 Effects on the evolution of orbital elements (top panel) and the deviation from GR (bottom panel). Each case is integrated independently. The cases are (i) full MOG (solid blue line), (ii) fifth-force contribution alone (solid red line), iii) space-time contribution alone (solid yellow line), and (iv) GR (dashed black line). Cases (ii) and (iii) are not physical in MOG, but are used for diagnostic purposes, as explained in Sect. 5.3. An explanation of the labels can be found in Table 1. For each orbital element µa, the definitions of ∆ µa and δµa are given in Eqs. (34) and (36), respectively.

Fig. 6 Evolution of the argument of pericentre with true anomaly. The full MOG theory is represented with a solid blue line, its 𝒮𝒯 effect with the solid yellow line, the DM case with the solid red line, and the GR case with the dashed black line. An explanation of the labels can be found in Table 1.

Fig. 7 Evolution of RV deviation for three values of α. The opaque lines (blue, red, and yellow, with spikes at the bottom) represent the RV absolute difference between (i) MOG EoM with MOG ℛ𝒮 and GR EoM with GR ℛ𝒮 (blue), (ii) MOG EoM without ℛ𝒮 and GR EoM without ℛ𝒮 (red), and (iii) GR EoM with MOG ℛ𝒮 and GR EoM with GR ℛ𝒮 (yellow). An explanation of the labels can be found in Table 2. Instead, the transparent lines, namely, the progressive red and blue lines and the constant yellow line that connect the peaks of the respective opaque lines, represent the evolution of the corresponding RV absolute difference peaks. The horizontal dashed black and grey lines represent the 1σ and 3σ precision of the instrument, respectively.

6 Impact of redshift MOG versus GR

In all the investigations until now, we only used the equation of motion (EoM) up to 1PN (Eq. (25)). However, for real data, the gravitational redshift (ℛ𝒮) has to be taken into account as well. The spatial coordinates in Fig. 3 are not affected by this phenomenon, unlike the radial velocity. Fig. 7 demonstrates three different comparisons to help us understand whether the gravitational redshift dominates the deviation of the local velocity due to the EoM. The comparisons are summarized in Table 2, where each of the two orbits (i.e. considering different configurations for EoM and ℛ𝒮) are integrated separately and their resulting RV is then subtracted.

Fig. 7 shows that although the deviation between the MOG and GR EoM (i.e. 𝒮𝒯 and $ℱℱ$) is progressive (red line), meaning that it increases with each period, the deviation due to the ℛ𝒮 is constant with time. Initially, the ℛ𝒮 has a stronger effect by one order of magnitude (yellow line), but the progressive deviation due to the EoM quickly kicks in and becomes dominant after some periods. Hence, for a sufficiently high α value and for a limited observation time, we have more opportunities to observe a discrepancy due to the ℛ𝒮 (e.g. α = 0.662), while for lower α values and for a long enough observation time, the EoM eventually produces a stronger deviation (e.g. α = 0.0662). For very low α values, we can calculate how much observation time is required in order for the total deviation (blue line) to reach the instrument precision (dashed horizontal lines).

7 Summary and discussion

We have performed a detailed investigation of the first post-Newtonian effects of MOG on stellar orbits around BHs and compared the results with those predicted by GR. We found similarities between MOG and GR with electromagnetism, and we quantified MOG signatures on observables and orbital elements.

Namely, we confirmed the orbital precession (∆ω) formula for MOG presented by Della Monica et al. (2021, 2023) using an alternative derivation (Eq. (35)). We further settled the disagreement that arose with Turimov (2022), in Appendix. A.

We proved for the first time, to the best of our knowledge, an exact equivalence between dynamics around MOG and Reissner–Nordström BHs by proposing a fine-tuning (Eqs. (15) and (18)).

On the side of observable signatures, our analysis showed that for low values of α, deviations from GR are minimal and difficult to detect within current precision. However, for higher values, the differences become significant, particularly in the orbital precession and in the evolution of orbital elements such as the semi-latus rectum and eccentricity. These effects manifest prominently in the pericentre half of the orbit and become more pronounced over multiple orbital periods.

By isolating the contributions of the MOG space-time geometry (𝒮𝒯) and fifth force $(ℱℱ)$, we found that the two components can partially cancel each other out, leading to minimal deviations in the overall orbit. However, for high α values, the cancellation is incomplete, resulting in noticeable differences from GR predictions. This behaviour highlights the non-linear nature of MOG and the complex interplay between its components.

Additionally, we compared MOG to a DM scenario with an extended mass distribution. Interestingly, we found that certain MOG effects, particularly those related to space–time geometry, produce features similar to those seen by DM, such as the characteristic behaviour of the argument of pericentre around the apocentre. This suggests a potential degeneracy between MOG and DM, which might complicate efforts to distinguish between the two using current observational data.

Finally, we explored the effect of gravitational redshift (ℛ𝒮) on radial velocity measurements. Our results indicate that while the redshift effect dominates initially, deviations due to the MOG EoM become significant over longer observation periods. This implies that precise and extended monitoring of stellar orbits is necessary to fully capture the MOG effect.

In conclusion, while MOG introduces subtle but potentially detectable deviations from GR, especially for higher values of α, it remains challenging to distinguish between MOG and DM effects. Future observations, particularly those involving long-term monitoring of stellar orbits and high-precision measurements of radial velocity, will be crucial in determining the validity of MOG as an alternative to dark matter and in refining our understanding of gravity in the strong-field regime.

Acknowledgements

We thank N. Aimar, T. Paumard, G. Perrin and F. Vincent for fruitful discussions and access to the Galactic Centre orbit model code library of Paris Observatory/LESIA. This scientific paper was supported by the Onassis Foundation – Scholarship ID: F ZU 042/2 2025-2026. I. Liodis would also like to thank the Lilian Voudouri Foundation for their master’s scholarship, the Onas-sis Foundation for their master’s and doctoral scholarship, as well as Deutsches Elektronen-Synchrotron (DESY).

References

Appendix A Fifth force and electromagnetism

We delve into a detailed analysis of the choice of sign in Eq. (11), due to the confusion arose in the comment papers Turimov (2022) and Della Monica et al. (2023).

Lopez Armengol & Romero (2017) studied the STVG-Kerr space time and invoked the resemblance of the vector potential ϕ_µ around a black hole, where µ_ϕ = 0 with the electromagnetic four potential A_µ of the Einstein-Maxwell theory (Misner et al. 1973). For a black hole with zero spin, the only non-zero covariant component of the vector potential is always negative, that is, $${\varphi }_0=-\frac{Q_{ff}}{r},$$(A.1)

where Q_ff is the fifth force charge of the black hole and r the Schwarzschild-like radial coordinate. Note that we now use natural units, since we only care about the direction of the extra force. It is important that Lopez Armengol & Romero (2017) and Misner et al. (1973) used the signature (-,+,+,+), from now on called 𝒮₋. This is of course in accordance with electromagnetism in special relativity Srednicki (2007), since the four potential around a point source of charge Q_em is defined as $$A^{\mu }=(\Phi ,A),$$(A.2)

where Φ = Q_em/r, leading to $$A_0=-\frac{Q_{em}}{r}\text{for\hspace{0.17em}}{\mathcal{S}}_{-},$$(A.3)

due to the metric with 𝒮_- signature.

If one chooses the (+,-,-,-) signature (namely 𝒮₊), the definition of the contravariant components of the four potential remains the same as in Eq. (A.2) (Peskin & Schroeder 1995), but now the covariant time component is given by $$A_0=\frac{Q_{em}}{r}\text{for\hspace{0.17em}}{\mathcal{S}}_{+},$$(A.4)

and the spatial components change sign. Note that in curved space-time, the covariant time components of the four potential have the same form with Eq. (A.3) or (A.4), depending on the signature (see Misner et al. (1973)).

Going back to GR, it can be easily shown that in both cases, that is, using either 𝒮₋ or 𝒮₊, the physics is of course the same, obeying to the following equation of motion (see Misner et al. (1973); Mehdizadeh & Jalili (2016)): $${\stackrel{\cdot }{U}}^{\mu }+{\Gamma }_{\alpha \beta }^{\mu }{U}^{\alpha }{U}^{\beta }=\frac{q_{em}}{m}{F}_{\nu }^{\mu }{U}^{\nu },$$(A.5)

where U^µ is the four velocity of a particle with mass m and charge q_em, moving in a space-time with Christoffel symbols Γ^µαβ and in an electromagnetic field F_µν = ∂_µA_ν ∂_νA_µ. Finally, the dot denotes derivative with respect to an affine parameter. For a point charge source placed at the origin, the only non zero components of the Faraday tensor are $$F_{10}=-F_{01}={\partial }_r{A}_0=\{\begin{array}{cc}+Q_{em}{r}^{-2}& \text{for\hspace{0.17em}}{\mathcal{S}}_{-}\\ -Q_{em}{r}^{-2}& \text{for\hspace{0.17em}}{\mathcal{S}}_{+}\end{array}$$(A.6)

while the radial component of the inverse metric can be expressed as $$g^{11}=\{\begin{array}{cc}+{\Psi }_{\text{RN}}(r)& \text{for\hspace{0.17em}}{\mathcal{S}}_{-}\\ -{\Psi }_{\text{RN}}(r)& \text{for\hspace{0.17em}}{\mathcal{S}}_{+}\end{array}$$(A.7)

where Ψ_RN follows the Reissner–Nordström metric (Reissner 1916; Nordström 1918). Therefore, the right hand side of the radial component of Eq. (A.5) is written as $$\frac{q_{em}}{m}{g}^{11}{F}_{10}{U}^0=+\frac{q_{em}{Q}_{em}}{m}\frac{\Psi (r)}{r^2}{U}^0,$$(A.8)

considering 𝒮₋ and 𝒮₊. The only thing that defines the sign of this expression outside the horizon is therefore the product q_emQ_em, which is positive for like charges leading to repulsive force, and negative for unlike charges leading to attracting force.

The $ℱℱ$ near a MOG BH has an equivalent picture if one identifies $$\begin{array}{ccc}{F}_{\mu \nu }& \to & B_{\mu \nu },\\ A_{\mu }& \to & {\varphi }_{\mu },\\ Q_{em}& \to & Q_{ff},\\ q_{em}& \to & q_{ff},\\ {\Psi }_{\text{RN}}(r)& \to & {\Psi }_{\text{MOG}}(r).\end{array}$$(A.9)

However, now both charges are always positive (Moffat 2015), and therefore, the $ℱℱ$ is always repulsive, as desired.

We highlight that the aforementioned correspondence leads to $${\varphi }_0=+\frac{Q_{ff}}{r},$$(A.10)

for the signature that we used (𝒮₊), and that a different choice of sign would require either a non-physical absorption of the minus inside a charge, or to postulate a relative minus sign in the action (22).

Retrospectively, it is evident that the claim of Turimov (2022) about the inconsistency of the signature used in Della Monica et al. (2021), is correct. However, Turimov (2022) gave an expression for the orbital precession using 𝒮₊ and a negative sign in the vector field, leading to an attracting force as demonstrated in this section, and as correctly pointed out in Della Monica et al. (2023). Finally, the latter study seems to have postulated a different sign in the $ℱℱ$, in order to make it repulsive, which of course results in the correct equation of motion.

Appendix B Derivation of the 1PN Lagrangian

The first term of Eq. (23) can be obtained using the method followed by Gainutdinov (2020), where the line element is given by (9) and (10). To comply with the notation in Gainutdinov (2020), Ψ(r) is written as $$\Psi (r)=1-2(1+\alpha )\frac{\mu }{r}+\alpha (1+\alpha )\frac{{\mu }^2}{r^2},$$(B.1)

with µ = G_N M/c². One should now choose a proper coordinate transformation to go to isotropic coordinates with radial element ρ. Let the arbitrary transformation $$r=\rho f(\rho ),\text{so\hspace{0.17em}that}{\Psi }^{-1}(\rho ){\left(\frac{dr}{d\rho }\right)}^2=f^2(\rho ).$$(B.2)

For positive α, the differential equation obtained has the following solution $$f(\rho )=1+\frac{(1+\alpha )\mu }{\rho }+\frac{(1+\alpha ){\mu }^2}{4{\rho }^2},$$(B.3)

and the line element is therefore written as $$\begin{array}{l}d{s}^2=\Psi (\rho )c^{2}d{t}^2-f^2(\rho )\left(d{\rho }^2+{\rho }^{2}d{\Omega }^2\right)\Rightarrow \hfill \\ d{s}^2=\Psi (\left|x\right|)c^{2}d{t}^2-f^2(\left|x\right|)d{x}^2,\hfill \end{array}$$(B.4)

where x denotes the spatial position vector in Cartesian-like coordinates. Now, writing φ_N = −GNM/ρ = −c²μ∕∣x∣ and expanding Ψ(∣x∣) and f²(∣x∣) with an accuracy of 𝒪(c⁻⁶) and 𝒪(c⁻⁴), respectively, yields $$\Psi (|x|)=1+2(\alpha +1)\frac{{\phi }_N(|x|)}{c^2}+(\alpha +1)(3\alpha +2)\frac{{\phi }_N^2(|x|)}{c^4}+\mathcal{O}(c^{-6}),$$(B.5a) $$\Psi (|x|)=1+2(\alpha +1)\frac{{\phi }_N(|x|)}{c^2}+(\alpha +1)(3\alpha +2)\frac{{\phi }_N^2(|x|)}{c^4}+\mathcal{O}(c^{-6}),$$(B.5b)

For simplicity, from now on φ_N(|x|) will be referred to as just φ_N, implicitly having a dependence on |x|. The line element is now written as $$\begin{array}{ll}\frac{1}{c}\frac{ds}{dt}=\hfill & 1-\frac{{\dot{x}}^2}{2c^2}\left[1+\frac{{\dot{x}}^2}{4c^2}-3(\alpha +1)\frac{{\phi }_N}{c^2}\right]\hfill \\ \hfill & +(\alpha +1)\frac{{\phi }_N}{c^2}\left[1+(2\alpha +1)\frac{{\phi }_N}{2c^2}\right]+\mathcal{O}(c^{-6})\hfill \end{array}$$(B.6)

where the square root has been taken to an accuracy of 𝒪(c⁻⁶). Now, one can identify the geometrical term of the action (23) as −c ds/dt, leading to $$\begin{array}{ll}-m{c}^2\frac{d\tau }{dt}=\hfill & -m{c}^2+m\frac{{\dot{x}}^2}{2}\left[1+\frac{{\dot{x}}^2}{4c^2}-3(\alpha +1)\frac{{\phi }_N}{c^2}\right]\hfill \\ \hfill & -m(\alpha +1){\phi }_N\left[1+(2\alpha +1)\frac{{\phi }_N}{2c^2}\right]+\mathcal{O}(c^{-4}).\hfill \end{array}$$(B.7)

Finally, in the isotropic coordinates, the $ℱℱ$ term of the action (23) is written as $$-m{c}^2\alpha \frac{G_{N}M}{r(\rho )c^2}=m\alpha {\phi }_N+m\alpha (1+\alpha )\frac{{\phi }_N^2}{c^2}+\mathcal{O}(c^{-4}).$$(B.8)

The full Lagrangian, including the $ℱℱ$ (Moffat & Toth 2009), can be now read off from the action in the Cartesian-like coordinates, with the result given by Eq. (24).

Appendix C Auxiliary

The transformation rules between the osculating orbital elements and position and velocity in Cartesian coordinates, at any instant of time, are given by: $$r=\left(r^X,r^Y,r^Z\right)\text{and}v=\left(v^X,v^Y,v^Z\right)$$(C.1)

Where $$\begin{array}{ll}{r}^X=\hfill & r\left(\cos \Omega \cos (\omega +f)-\cos \iota \sin \Omega \sin (\omega +f)\right)\hfill \\ r^Y=\hfill & r\left(\sin \Omega \cos (\omega +f)+\cos \iota \cos \Omega \sin (\omega +f)\right)\hfill \\ r^Z=\hfill & r\sin \iota \sin (\omega +f)\hfill \end{array}$$

and $$\begin{array}{ll}{v}^X\hfill & =-C_p(\cos \Omega (\sin (\omega +f)+e\sin \omega )\hfill \\ \hfill & =C_p\left(+\cos \iota \sin \Omega (\cos (\omega +f)+e\cos \omega )\right)\hfill \\ v^Y\hfill & =-C_p(\sin \Omega (\sin (\omega +f)+e\sin \omega )\hfill \\ \hfill & =(-\cos \iota \cos \Omega (\cos (\omega +f)+e\cos \omega )\hfill \\ v^Z\hfill & =C_p\sin \iota \left(\cos (\omega +f)+e\cos \omega \right),\hfill \end{array}$$

with r = p/(1 + e cos f ) and $C_p=\sqrt{Gm/p}$. For their derivation see, for instance, (Poisson & Will 2014, Sect. 3.2.5).

All Tables

All Figures

	Fig. 1 Definition of the Euler angles (Ω, ι, ω) giving the orbital orientation and of the true anomaly f giving the current position within the orbit. The static fundamental frame (X, Y, Z) (blue) is adapted to the observables, while the co-rotated Gaussian frame (n,λ, ez) (orange) is adapted to the symmetry of the (perturbed) two-body problem (reproduced from Heißel, et al. 2022).
In the text

	Fig. 2 Oribt of star S2 (yellow star) around the BH (black dot) in the centre of our Galaxy for GR (left panel) and an extreme MOG case (right panel).
In the text

Fig. 3 Evolution of the observables’ (RA, Dec, and RV in the first, middle, and bottom row, respectively) absolute difference between GR and i) full MOG (solid blue lines), ii) fifth force (solid red lines), and iii) space-time (dashed yellow lines) contribution for different values of α. We fixed α as α = 6.62 · 10−3 in the left column, α = 6.62 · 10−2 in the middle one and α = 6.62 10−1 in the right column. An explanation of the labels can be found in Table 1. The horizontal dashed lines correspond to the 1σ (black) and 3σ ·(grey) accuracy of the instruments. For the spatial observables (RA and Dec), these are 50 µas and 150 µas, respectively, and for the radial velocity (RV), they were fixed to 10 km/s and 30 km/s.

In the text

	Fig. 4 Orbital element evolution with true anomaly. Here, ∆p stands for the offset of p from its GR value at f = 0, and analogously for the other elements. The dashed black line represents GR, and the solid coloured lines represent a MOG case with a particular α value in the range (0.001, 1). For each orbital element µa, the definition of ∆µa is given in Eq. (34).
In the text

Fig. 5 Effects on the evolution of orbital elements (top panel) and the deviation from GR (bottom panel). Each case is integrated independently. The cases are (i) full MOG (solid blue line), (ii) fifth-force contribution alone (solid red line), iii) space-time contribution alone (solid yellow line), and (iv) GR (dashed black line). Cases (ii) and (iii) are not physical in MOG, but are used for diagnostic purposes, as explained in Sect. 5.3. An explanation of the labels can be found in Table 1. For each orbital element µa, the definitions of ∆ µa and δµa are given in Eqs. (34) and (36), respectively.

In the text

	Fig. 6 Evolution of the argument of pericentre with true anomaly. The full MOG theory is represented with a solid blue line, its 𝒮𝒯 effect with the solid yellow line, the DM case with the solid red line, and the GR case with the dashed black line. An explanation of the labels can be found in Table 1.
In the text

Fig. 7 Evolution of RV deviation for three values of α. The opaque lines (blue, red, and yellow, with spikes at the bottom) represent the RV absolute difference between (i) MOG EoM with MOG ℛ𝒮 and GR EoM with GR ℛ𝒮 (blue), (ii) MOG EoM without ℛ𝒮 and GR EoM without ℛ𝒮 (red), and (iii) GR EoM with MOG ℛ𝒮 and GR EoM with GR ℛ𝒮 (yellow). An explanation of the labels can be found in Table 2. Instead, the transparent lines, namely, the progressive red and blue lines and the constant yellow line that connect the peaks of the respective opaque lines, represent the evolution of the corresponding RV absolute difference peaks. The horizontal dashed black and grey lines represent the 1σ and 3σ precision of the instrument, respectively.

In the text