Testing horndeski gravity with S2 star orbit

MNRAS 519, 1981–1988 (2023) https://doi.org/10.1093/mnras/stac3648 Advance Access publication 2022 December 15 Testing horndeski gravity with S2 star orbit R. Della Monica ,1 ‹ I. de Martino ,1 D. Vernieri2,3,4 and M. de Laurentis2,3 1 Universidad de Salamanca, Departamento de Fisica Fundamental, P. de la Merced S/N, 37008, Salamanca, ES di Fisica, Universitá di Napoli ‘Federico II’, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, I-80126, Napoli, Italy 3 INFN Sezione di Napoli, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, I-80126, Napoli, Italy 4 Scuola Superiore Meridionale, Largo San Marcellino 10, I-80138, Naples, Italy 2 Dipartimento Accepted 2022 December 9. Received 2022 November 17; in original form 2022 September 13 We have explored a completely new and alternative way to restrict the parameter space of Horndeski theory of gravity. Using its Newtonian limit, it is possible to test the theory at a regime where, given its complexity and the small magnitude of the expected effects, it is poorly probed. At Newtonian level, it gives rise to a generalized Yukawa-like Newtonian potential which we have tested using S2 star orbit data. Our model adds five parameters to the General Relativity model, and the analysis constrains two of them with unprecedented precision to these energy scales, while it only gives an exclusion region for the remaining parameters. We have shown the potential of weak-field tests to constrain Horndeski gravity opening, as a matter of fact, which is a new avenue that deserves to be further, and deeply, explored in the near future. Key words: Galaxy: centre – gravitation – celestial mechanics. metric gμν and one scalar field φ, is given by the following terms: 1 I N T RO D U C T I O N In the last decades, several extensions of General Relativity (GR) have been proposed in order to overcome its historical limitations, by adding extra ingredients to the gravitational action. Indeed taking into account higher-order derivatives or extra dynamical fields enriches the dynamics and allows to get new phenomenology. Hereafter, we will consider the possibility to add one extra scalar degree of freedom directly coupled to the metric tensor. In this context, we will take into account the very well-known and general class of Horndeski theory of gravity (Horndeski 1974; Kobayashi, Yamaguchi & Yokoyama 2011) which has been extensively investigated so far. Such theories have been proposed as the most general ones leading to second-order field equations for the metric tensor and the scalar field, even if they have also been extended to more general schemes without increasing the number of degrees of freedom (Horndeski 1974). Horndeski theory has been studied both in cosmological and astrophysical frameworks (for a comprehensive review, see Kobayashi (2019) and the references therein). In the latter case, several studies have been conducted in particular concerning black holes (Maselli et al. 2015; Babichev, Charmousis & Lehébel 2017) and neutron stars (Cisterna et al. 2016; Maselli et al. 2016). The gravitational action of Horndeski theory takes the following form (Kobayashi et al. 2011): S= 5   i=2 √ d 4 x −g Li [gμν , φ] + Sm [gμν , χm ], (1) where Sm is the matter action and χ m collectively denotes all matter fields. The gravitational part of the action, which depends on the ⋆ E-mail: rdellamonica@usal.es © 2022 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society L2 = G2 (φ, X) , L3 = −G3 (φ, X)φ,   L4 = G4 (φ, X)R + G4X (φ, X) (φ)2 − (∇μ ∇ν φ)2 , 1 L5 = G5 (φ, X)Gμν ∇ μ ∇ ν φ − G5X (φ, X) 6   × (φ)3 − 3(φ)(∇μ ∇ν φ)2 + 2(∇μ ∇ν φ)3 . (2) Generally speaking, the functions G2 , G3 , G4 , and G5 are free functions of the scalar field φ and its kinetic term X. Each choice of these functions determines a distinct gravity theory which is determined once a suitable mapping is imposed. Notice that it is possible to map within the general class of Horndeski theories, not just those with a scalar field appearing explicitly in the gravitational action and in the corresponding field equations, but also theories (like f(R) gravity or f (G) Gauss–Bonnet gravity) that inherit higher-order field equations but that somehow admit an equivalent scalar field representation (Kobayashi et al. 2011). Additionally, the full Newtonian and Post-Newtonian expansions of the theory have been performed in 2015 (Hohmann 2015), leading to an effective general Yukawa-like potential. Nevertheless, no test of the weak-field limit of Horndeski theory has been made so far by using such a potential, which can indeed allow to constrain the theory at astrophysical scale where it is not yet probed. Hence, we will focus on using the orbital data of the S2 star (Schödel et al. 2002; Eisenhauer et al. 2003; Ghez et al. 2003) orbiting the supermassive black hole (SMBH) in the Galactic Centre (Ghez et al. 1998; Genzel, Eisenhauer & Gillessen 2010) to bound potential deviation from GR into the Horndeski gravity framework. Although the approach is not new (Hees et al. 2017; De Martino, della Monica & De Laurentis 2021; Della Monica, de Martino & de Laurentis 2022; Della Monica & de Martino 2022), it has never been applied to such a large class of theories and, therefore, it represents a path never Downloaded from https://academic.oup.com/mnras/article/519/2/1981/6909075 by guest on 11 January 2023 ABSTRACT 1982 R. Della Monica et al. 2 GENERALIZED MODIFIED NEWTONIAN P OT E N T I A L The Newtonian limit of Horndeski gravity theory is calculated through a perturbative expansion of the field equations around a Minkowski background ημν and a constant cosmological background value  of the scalar field (Hohmann 2015): 1 gμν = ημν + hμν , φ =  + ψ , X = − ∇ μ ψ∇μ ψ. (3) 2 Besides assuming that the background is homogeneous and isotropic, we thus also assume that it is stationary, and additionally that the four free functions G2 , G3 , G4 , and G5 , can be expanded in Taylor series as Gi (φ, X) = ∞  Gi(m,n) ψ m X n , (4) m,n=0 where the coefficients Gi(m, n) are given by   1 ∂ m+n  Gi(m,n) = , G ( φ, X) i  m n m!n! ∂ φ ∂ X φ=,X=0 (5) where each term Gi(m, n) ψ m Xn is of the order O(ψ m+2n ). Notice that, in light of the multi-messenger event GW170817/GRB170817 A, the surviving theories are those for which G4, X = 0 and G5 = 0 (Creminelli & Vernizzi 2017). Additionally, it is worth to notice that G4, X and G5 do not appear at Newtonian order (we refer the reader to equation (28) in Hohmann (2015)). For a mass point source the solution at Newtonian order is (Hohmann 2015)   c1 cψ −mψ r M (2) c2 + 2 (e − 1) , (6) h00 (r) = 4π r mψ where the constants have been defined as follows: G4(1,0) G2(2,0) , c1 = −2 G4(0,0) ϒ G24(1,0) 1 c2 = + , 2G4(0,0) 2G24(0,0) ϒ MNRAS 519, 1981–1988 (2023) −2G2(2,0) G4(1,0) −1 , cψ = ϒ , (8) ϒ 2G4(0,0) where for sake of simplicity we have introduced the following definition mψ = ϒ ≡ G2(0,1) − 2G3(1,0) + 3 G24(1,0) . G4(0,0) (9) Following the seminal work of Kobayashi et al. (2011), and the more recent results of Hohmann (2015), equation (8) requires ϒ > 0, which implies G2(0, 0) < 0, reducing the allowed parameter space of the theory. Replacing equations (7) and (8) in the metric coefficient (6), one can express the time component of the metric perturbation in the following form:   M 1 1 G4(1,0) 2 −1 −mψ r (2) h00 (r) = , (10) ϒ e + 4π r 2G4(0,0) 2 G4(0,0) where G4(0, 0) plays the role of the inverse of gravitation constant G. In a general fashion, Horndeski theory of gravity has a unique limit that recovers GR. In particular, by looking at the gravitational action in equations (1) and (2), it is clear that GR (L = R/16π G) is recovered if and only if all the functions G2 , G3 , and G5 are identically zero (and hence their Taylor expansion coefficient are zero, as well), and G4 = (16π G)−1 . This is the only physical limit to GR as it involves the cancellation of the additional scalar field φ introduced to describe the gravitational interaction. On the other hand, at a metric level, the Newtonian approximation of GR for the metric of space–time around a massive point source is given by: h(2) μν (r) = diag 2GM , 0, 0, 0 . r (11) If we compare equation (10) with equation (11), we see that we can recover the expression of GR by putting G4(0, 0) = 1/16π G and performing one of the following limits: (A) G4(1, 0) → 0 which recovers equation (11) for any value of the other Taylor coefficients; (B) G2(2, 0) → ±∞ which brings mψ → ∞ and hence the cancellation of the Yukawa-like term in equation (10); (C) Any combination of the coefficients G2(0, 1) , G3(1, 0) , and G4(1, 0) which brings ϒ → 0, and hence mψ → ∞, driving the entire Yukawalike term to an indeterminate form of the kind e−∞ /0 whose limit is 0; (D) The simultaneous limit G2(0, 1) , G2(2, 0) , G4(1, 0) , G3(1, 0) , →0, which, again, recovers the expression in equation (11) for the metric coefficient. It is evident that the presence of multiple GR limits, among which only the case (D) corresponds to the physical limit of GR of the theory itself, is a spurious effect due to having performed the Newtonian approximation. Care must be taken, hence, in fixing the values of the coefficients in equation (10) as only one region of the parameter space corresponds to being close to GR, even though at a metric (and thus geodesic) level multiple regions result in a GR-like orbit. The above considerations involve notable degenerations in the parameter space, as different combinations of these parameters return the same physical limit, which cannot be broken without being able to impose a priori limits on these parameters. 3 O R B I TA L M O D E L (7) Let us start summarizing here the technique used to build the orbit of S2 in Horndeski theory of gravity. From the line element in Downloaded from https://academic.oup.com/mnras/article/519/2/1981/6909075 by guest on 11 January 2023 explored that can potentially lead to new bounds on the underlying gravitational theory. More in general, the Yukawa interaction arises in many different scenarios from unification theories (Kaplan & Wise 2000) to braneworld scenarios (Dvali, Gabadadze & Porrati 2000), modified gravity (De Martino et al. 2021; Della Monica et al. 2022), and certain models of dark matter (Carroll et al. 2009; Carroll, Mantry & Ramsey-Musolf 2010). Hees et al. (2017) used the phenomenological framework of the fifth force to constrain the coupling strength and the range or mass of the scalar interaction. Indeed, using S stars orbiting around the SMBH in the Galactic Centre, they bounded the absolute value of the strength of the fifth force to be <0.016 on a scale length of ∼150 AU. Although mapping these constraints in a specific model is not always a straightforward task, De Martino et al. (2021) found them to agree with constraints of the Yukawa interaction in f(R)-gravity. This article is organized as follows: in Section 2, we introduce the generalized modified Newtonian limit of Horndeski gravity; in Section 3, we describe the orbital model for the S2 star in this potential and provide some toy models to illustrate the effects arising in this framework; in Section 4, we give a description of the set of observational data used in order to constrain the additional parameters of the theory and our data analysis procedure; in Section 5, we report the results of our analysis; and finally in Section 6, we give our final remarks. Horndeski gravity with S2 1983 Table 1. The full set of orbital parameters for the orbit of the S2 star as estimated by Gillessen et al. (2017) used as priors of our analysis (for these parameters, we considered gaussian distributions whose FWHM are 10 times the experimental errors reported here) and the reference frame priors on x• , y• , v x, • , v y, • , and v z, • from Plewa et al. (2015). In the bottom-right part of the table, we report the set of priors used for the parameters G4(0, 0) , G4(1, 0) , G3(1, 0) , G2(0, 1) , and G2(2, 0) in our analysis. We have left G4(0, 0) vary around its GR value in an interval of amplitude 10 per cent. The other intervals have been set heuristically, starting from some guessed values and then widening or narrowing the interval as needed. The parameter G2(2, 0) has been fixed to be negative and thus we assumed an additional prior to set ϒ > 0 (Kobayashi et al. 2011; Hohmann 2015). Parameter 106 M kpc mas ◦ ◦ ◦ T tp x• yr yr mas ⊙ Value Error Parameter Unit Value Error 4.35 8.33 125.5 0.8839 134.18 65.51 226.94 16.00 2018.33 −0.2 0.012 0.0093 0.044 0.000079 0.033 0.030 0.031 0.0013 0.00017 0.2 y• v x, • v y, • v z, • mas mas yr−1 mas yr−1 km s−1 0.1 0.02 0.06 0 0.2 0.2 0.1 5 G4(0, 0) G4(1, 0) G3(1, 0) G2(0, 1) G2(2, 0) M⊙ AU s−2 M⊙ AU s−2 M⊙ AU s−2 M⊙ AU s−2 M⊙ /AU s−2 equations (3) and (10), one can compute the geodesic equations, Equations (A1–A4) that describe the motion of test particles in the Newtonian limit. These equations can be integrated once initial conditions are assigned for the four space–time coordinates and the components of the initial tangent vector to the geodesic. In doing so, one has to take into account the normalization condition for the four-velocity. Since the metric coefficients do not explicitly depend on the time coordinate t, the particular value of t(0) does not affect the integrated orbit. For the S2 star, we set for convenience t(0) = tA , where tA corresponds to the time of the last apocentre passage of the star, which occurred in ∼2010. Finally, as regards the initial condition on r, φ, and the corresponding components of the velocity, we set values that are consistent with the choice t(0). Namely, we set the radial and angular position of the star and its velocity at the apocentre. These quantities are estimated by computing the Keplerian osculating ellipse at the apocentre, as given by the Keplerian elements measured for S2 in Gravity Collaboration (2020) and reported in Table 1. Once geodesics have been integrated numerically, we have projected both spatial and velocity vectors of the mock orbit in the reference frame of a distant observer in order to be able to compare it with the observational data. Additionally, in order to properly reconstruct synthetic orbits for S2, we have to take into account on the projected quantities the classical Rømer delay and the redshift on the radial velocity coming from Special Relativity and gravitational time dilation (we refer to Grould et al. (2017), Della Monica et al. (2022), De Martino et al. (2021), Della Monica & de Martino (2022) for a detailed analysis on how these effects can be quantified on our predicted orbits). Finally, we correct our synthetic orbits for systematic effects arising from the construction of the reference frame. More details on the integration procedure we followed can be found in Appendix A. Thus, our baseline model has 18 parameters, namely, six orbital parameters (the time of pericentre passage tp , the semi-major axis of the orbit a, its eccentricity e, the inclination i, the angle of the line of nodes , and argument of the pericentre ω), two source parameters that are the mass of the central object M• and its distance R• , five system reference parameters, i.e. (x• , y• ) related to a zero-point offset of the central mass with respect to the origin of the astrometric reference frame; (v x, • , v y, • ) describing a drift over time of the central mass on the astrometric reference frame; and v z, • related to an offset in the estimation of the radial velocity of the S2 Interval [0.95, 1.05](c4 /16π G) [−100, 100] [−100000, 100000] [−100000, 100000] [−10, 0] star when correcting to the Local Standard of Rest (LSR). Finally, we have five extra parameters from the expansion of Horndeski gravity at Newtonian order: G2(2, 0) , G2(0, 1) , G3(1, 0) , G4(0, 0) , and, G4(1, 0) . 3.1 Toy models As a first step to study the motion of test particles in the Newtonian expansion of Horndeski gravity for a spherically symmetric space– time, we have selected five toy models for S2 with the same orbital parameters (the ones in Table 1) and with varying values for the Horndeski Taylor coefficients. The choice of such parameters has been made in order to assess the deformation of the orbit as a function of 1/mψ (which has dimensions of a length) and G4(1, 0) , which, together, modulate the strength and scale length of the modification to the Newtonian potential in equation (3). In particular, we have fixed the GR value of G4(0, 0) = c4 /16π GN , while we have chosen values for the parameters G4(1, 0) , G3(1, 0) , G2(0, 1) , and G2(2, 0) that would give a value of 1/mψ that is a given fraction of the typical scale length of the S2 orbit (e.g. its semi-major axis, a ∼ 1000 AU). The chosen values are reported in Table 2 and the corresponding orbits from the numerical integration are shown in Fig. 1. In particular, the model Toy 1 represent a star moving under the influence of an unperturbed Newtonian potential, due to the fact that all parameters have been set to zero. Models Toy 2 and Toy 4 have G2(2, 0) = −0.01 M⊙ /AUs2 , while for Toy 3 and Toy 5 a value of G2(2, 0) = −0.001 M⊙ /AUs2 has been chosen. On the other hand, while Toy 2 and Toy 3 have a value G4(1, 0) = 10 M⊙ AU/s2 , a value of G4(1, 0) = 100 M⊙ AU/s2 has been fixed for Toy 4 and Toy 5. This choice results in an increasingly high scale length 1/mψ for the different models, associated with a strength of the Yukawa-like potential in equation (3), modulated by G4(1, 0) , that is increasingly higher. As can be seen from the integrated orbits in Fig. 1, the main effect related to the variation of the length 1/mψ is a modification of the orbit and, in particular, its pericentre shift on the orbital plane. 4 DATA A N D DATA A N A LY S I S In order to explore the parameter space of the coefficients arising in the Newtonian approximation of Horndeski gravity for a spherically symmetric space–time, we have employed publicly available MNRAS 519, 1981–1988 (2023) Downloaded from https://academic.oup.com/mnras/article/519/2/1981/6909075 by guest on 11 January 2023 M• R• a e i ω Unit 1984 R. Della Monica et al. Table 2. The set of parameters chosen to build our toy models. The second column reports the value of the scale length of the Horndeski modification to the gravitational potential at Newtonian order. Columns 3 to 7, on the other hand, provide with the particular set of parameters resulting in the chosen value of 1/mψ . Model (1) 1/mψ (AU) (2) G4(0, 0) (M⊙ AU/s2 ) (3) G4(1, 0) (M⊙ AU/s2 ) (4) G3(1, 0) (M⊙ AU/s2 ) (5) G2(0, 1) (M⊙ AU/s2 ) (6) G2(2, 0) (M⊙ /AUs2 ) (7) Toy 1 Toy 2 Toy 3 Toy 4 Toy 5 0.0 66.0 208.6 433.4 1370.5 8.0940 8.0940 8.0940 8.0940 8.0940 0 10 10 100 100 0 0 0 0 0 0 50 50 50 50 0 −0.01 −0.001 −0.01 −0.001 kinematic data for the S2 star. In particular, the data presented in Table 5 of Gillessen et al. (2017) provide with 145 astrometric data points of the recorded sky-projected position of S2 with respect to the ‘Galactic Centre (GC) infrared reference system’ (Plewa et al. 2015) and 44 measurements of its velocity along the line of sight, corrected for the LSR, which have been retrieved from line shift through spectroscopic observations. These data cover a period spanning from ∼1992 to ∼2016. This means that, given that the period of S2 is ∼16 years, a full orbit is covered, but no data are present for the crucial pericentre passage occurred in May 2018. The data that have been recorded with GRAVITY interferometer during this passage are not publicly available, but have been used to place significant constraints on the first post-Newtonian order of GR (Gravity Collaboration 2018, 2020). An independent analysis by the Galactic Center group at the University of California, Los Angeles, was able to detect the gravitational redshift at the S2, called S-02 in their works, pericentre using data collected at the Keck observatory (Do et al. 2019). Nevertheless, the available data are sufficient for our purpose, since our aim is to perform for the first time an exploration of the parameter space for Horndeski gravity at Newtonian order. MNRAS 519, 1981–1988 (2023) We have employed the Bayesian sampling of the posterior probability distribution of these parameters treated as random variables. In particular, this was done by performing a Markov Chain Monte Carlo (MCMC) analysis, implemented in the emcee (Foreman-Mackey et al. 2013) module. In order not to force the results to a Keplerian orbit, we have not fixed the orbital elements of S2 to the values derived in previous analysis on the same dataset (Gillessen et al. 2017), but we have left them as free parameters of our fit. Thus, the sampler draws random values for all the 18 aforementioned parameters, within appropriate prior distributions, and, for each parameter set, we perform a synthetic orbit for S2, as illustrated above. The agreement between the numerically integrated orbit and the observational data is assessed with the following log-likelihood function: − 2 log L =  i + 2 i Xobs − Xi i σX, obs VZi ,obs − VZi σVi Z ,obs + i Yobs − Yi σYi ,obs 2 2 , (12) Downloaded from https://academic.oup.com/mnras/article/519/2/1981/6909075 by guest on 11 January 2023 Figure 1. The integrated orbits for the toy models in Table 2. The pink segment illustrates the varying scale length 1/mψ . The plots report how modification of the orbit is predicted, as the coefficient of the Newtonian expansion of Horndeski are varied. The most noticeable effect that we can visually see on the plots is the pericentre shift on the orbital plane of the star. The quantity 1/mψ behaves like a scale and it gets more prominent for smaller values of 1/mψ . Horndeski gravity with S2 where (Xobs , Yobs , VZ,obs ) are the observed positions and radial velocities for S2, while (X, Y, VZ ) are the ones from our numerical integration. Convergence of the algorithm is assessed with the estimation of the autocorrelation time of the Markov chains (Goodman & Weare 2010). The particular priors used for our analysis are given by a set of univariated Gaussian distributions on the orbital parameters whose central values are fixed on the best-fitting values of Gillessen et al. (2017) and whose full weight half maximum are chosen to be ten times the uncertainty on these parameters from Gillessen et al. (2017) (reported in Table 1). The parameters related to the systematic effects in the construction of the reference frame have been assigned Gaussian priors from an independent analysis in Plewa et al. (2015). Finally, for the five parameters of the Newtonian expansion of Horndeski gravity, we used uniform priors whose amplitudes have been set heuristically. All priors are listed in Table 1. 5 R E S U LT S The posterior distributions of the 18 parameters of our orbital model for the S2 star in Horndeski gravity are reported in the online supplementary materials. More specifically, we report the 68, 95, and 99.7 per cent confidence intervals for all the parameters. It is worth mentioning that all the orbital elements and the reference frame parameters for the S2 star are bounded and recover, within 3σ (at most), all fiducial values reported in Table 1. Only the semi-major axis a and the inclination angle i represent an exception, being correlated with the parameter G4(0, 0) , and are only marginally compatible with their counterpart obtained in a Newtonian framework in Gillessen et al. (2017). This is anyway rather expected since G4(0, 0) modulates the Newtonian gravitational constant G, hence affecting the posterior distribution of the proposed values of the semi-major axis and of the inclination angle. Coming to the additional Horndeski parameters, we place constraint on G4(0, 0) and G4(1, 0) , while the other parameter remains unconstrained. Nevertheless, we were able to place a constraint on a combination of the remaining parameters and determined an exclusion region. In detail, Fig. 2 is an inset of the corner plot of the full posterior reporting the 68, 95, and 99.7 per cent confidence intervals of G4(0, 0) and G4(1, 0) parameters, and their corresponding posterior distributions. As reported in the Figure, we 0.16 +32 −2 find G4(0,0) = 8.17+ and G4(1,0) = −2−26 M⊙ AU s−2 . −0.17 M⊙ AU s The parameter G4(0, 0) is fully compatible with its expected value in GR (which is depicted as a vertical red line in the top panel), and the parameter G4(1, 0) is compatible with zero at 1σ , reducing the gravitational potential to the Newtonian one: h00 = 8πGM4(0,0) r . The result was rather expected but still surprising because the orbital motion data of S2 star, which we are using to test the Newtonian limit of Horndeski theory of gravity, cannot distinguish GR from Newton theory either (Gillessen et al. 2017). Previous tests of the weak-field limit of Horndeski gravity strongly constrained G4(0, 0) , but they consider a subclass of the whole theory by setting G2(2, 0) = 0, which means mψ = 0 (Hou & Gong 2018). In this case, the Yukawa-like term disappears reducing itself to an additional constant in equation (10). On the other hand, none was able to place constraints on G4(1, 0) in the weak field limit. Additionally to those new bounds, looking at the contours of the G2(0, 1) − G3(1, 0) slice of the whole parameter space, we can see that even though we are not able to place a definitive constraint on either one of the two parameters, we are able to exclude a region of the parameter space that results in a lower limit on the combination G2(0, 1) − 2G3(1, 0) > 6200M⊙ AU s−2 as shown in Fig. 3. Finally, we were able to bound the scale length m−1 ψ . First, we randomly sampled m−1 from the posterior of Horndeski ψ gravity parameters carrying out a 1000 Monte Carlo simulation to estimate the distribution of mψ and its variance. As shown in Fig. 4, +265 we found m−1 ψ = 160−75 AU. Noticeably, this represents the best constraint up to date on m−1 ψ , since other analyses only place a lower limit m−1 ψ ≥ 9 AU using binary systems composed by a pulsar plus a white dwarf (Dyadina, Avdeev & Alexeyev 2019). Finally, let us point out that, even though other constraints exist in literature for a phenomenological Yukawa-like modification to the Newtonian potential using the S2 star at the GC (see, e.g. Hees et al. (2017)), due to the specific dependence of equation (10) on the Taylor parameter of the Horndeski functions, there is no way to consistently map constraints on the former potential into bounds on the latter. Indeed, phenomenological Yukawa-like modifications are parametrized by a scale length (λ) and a modulation factor (α); these are treated (both conceptually and statistically) as independent parameters. The equivalent scale length 1/mψ in our metric (10), on the other hand, is not independent from the coefficient 1/2ϒ(G4(1, 0) /G4(0, 0) )2 of the exponential term in the metric, due to the fact that mψ itself depends on all the four free parameters of the Newtonian limit studied here. This, along with the intrinsic non-linearity of the relation between the Horndeski parameters and those in the Yukawa-like potential, makes it conceptually impossible to map pre-existing bounds on α and λ into bounds on G4(0, 0) , G4(1, 0) , G3(1, 0) , and G2(0, 1) . 6 DISCUSSIONS AND CONCLUSIONS Horndeski theory of gravity represents the most general scalar–tensor theory (including – but not limited to – Brans–Dicke, f(R) gravity, f (G) Gauss–Bonnet gravity, covariant Galileons, and chameleons among the others), with an additional scalar degree of freedom coupled to the metric and leading to second order equations of motion. While formulated almost 50 years ago, this model has only recently been back in vogue and extensively investigated at cosmological scales (Horndeski 1974). Nevertheless, it lacks a deep investigation in the weak field limit. On one hand, this is due MNRAS 519, 1981–1988 (2023) Downloaded from https://academic.oup.com/mnras/article/519/2/1981/6909075 by guest on 11 January 2023 Figure 2. Focus on the posterior distribution of the Taylor coefficients of Horndeski gravity. The vertical red line in the top panel depicts the expected value in GR of the Taylor coefficient G4(0, 0) . 1985 1986 R. Della Monica et al. Table 3. Results of our posterior analysis on all the parameters of our orbital model for the S2 star in the Newtonian expansion of Horndeski’s gravity. In particular, we report the 68% confidence interval for all the orbital and reference frame parameters; the 68% confidence interval for G4(0, 0) (which is bounded and perfectly compatible to the Newtonian value c4 /16π GN ) and for G4(1, 0) (which is centred on 0). The parameter G2(2, 0) results unbounded from our analysis, while we are only able to put a constraint on the combination G2(0, 1) − 2G3(1, 0), which appears as an exclusion region on the corresponding contour plot. Finally, we report the 68% confidence interval for the scale length 1/mψ . Parameter Unit Best-fit Parameter Unit Best-fitting R• kpc y• mas T yr 0.13 8.53+ −0.12 0.23 15.95+ −0.15 +0.11 0.33−0.08 +0.00052 0.12323−0 .00048 +0.0027 0.8821−0.0027 0.18 135.12+ −0.17 +0.54 226.09−0.54 0.54 64.25+ −0.44 1.6 −0.4+ −2.4 1.8 0.2+ −2.1 tp − 2018 a yr mas i ( ◦) ( ◦) ω ( ◦) x• mas v y, • v z, • G4(0, 0) G4(1, 0) G2(0, 1) − 2G3(1, 0) mas yr−1 mas yr−1 km s−1 M⊙ M⊙ M⊙ AU s−2 AU s−2 AU s−2 G2(2, 0) M⊙ AU−1 s2 1/mψ AU 0.12 −0.02+ −0.08 0.09 0.00+ −0.11 3.9 7.5+ −4.0 0.16 8.17+ −0.17 32.0 −2.0+ −26.0 6200 Unbounded 265 160+ −75 Figure 3. Ninety-five per cent confidence level exclusion regions of the G2(0, 1) − G3(1, 0) slice of the parameter space from our posterior analysis. Figure 4. Posterior distribution on the inverse of mψ in AU derived from the posteriors of the four parameters in equation (8). Dashed vertical lines correspond to the median value of the posterior distribution, also reported in Table 3. to the fact that a weak field approximation was formulated only a few years ago in 2015 (Hohmann 2015). On the other hand, the intrinsic complexity of such a weak field limit, reflecting in a very high number of additional parameters with respect to GR, makes it rather complicated to extract valuable constraints from astrophysical dataset. Indeed, a previous study constraining the weak field limit of Horndeski theory was only able to bound G4(0, 0) in a subclass of models that set G2(2, 0) = 0, which means mψ = 0, using time-delay data from the Cassini mission and gravitational wave radiation from binary pulsar systems (Hou & Gong 2018). The only previous analysis investigating the whole theory uses binary systems to place a lower limit on the parameter m−1 ψ ≥ 9 AU (Dyadina et al. 2019). Therefore, a comprehensive analysis constraining the Taylor coefficients of the functions G2 , G3 , G4 , and G5 appearing in the weak field limit in equation (10) was fully missing so far. Using the orbital motion of the S2 star around the SMBH in the centre of the Milky Way, we have placed constraints on 0.16 −2 two Taylor coefficients: G4(0,0) = 8.17+ and G4(1,0) = −0.17 M⊙ AU s +32 −2−26 M⊙ AU s−2 ; and on a combination of other two coefficients G2(0, 1) − 2G3(1, 0) > 6200M⊙ AU/s2 . These results recover the value of G4(0, 0) expected in GR at 1σ , and set a new bound on the parameter G4(1, 0) that encodes the deviations from GR. Our constraint shows a full accordance with the GR value of G4(1, 0) = 0 within 1σ . Moreover, we also are capable to set a strong bound on the previously unbounded parameter m−1 ψ that tunes the scale on which the novel Yukawainteraction acts. This analysis does not only provide with unprecedented constraints on Horndeski gravity at astrophysical scales that had never been tested before, but also hints the potential of future astronomical observations to constrain the theory in both the weak field limit through the Post-Newtonian approximation that can be probed by the forthcoming GRAVITY data (Gravity Collaboration 2017), and the strong field limit that can be probed through the black hole solutions in Horndeski gravity using horizon-scale Event Horizon Telescope MNRAS 519, 1981–1988 (2023) Downloaded from https://academic.oup.com/mnras/article/519/2/1981/6909075 by guest on 11 January 2023 e v x, • Horndeski gravity with S2 observations (Doeleman et al. 2008; Broderick, Loeb & Narayan 2009). Joint observations exploring both weak and strong field limits will potentially lead to either endorse or rule out higher-order theories of gravity. AC K N OW L E D G E M E N T S DATA AVA I L A B I L I T Y Data used in this article are publicly available on the electronic version of Gillessen et al. (2017) at https://iopscience.iop.org/article /10.3847/1538-4357/aa5c41/meta#apjaa5c41t5. REFERENCES Babichev E., Charmousis C., Lehébel A., 2017, J. Cosmology Astropart. Phys., 2017, 027 Broderick A. E., Loeb A., Narayan R., 2009, ApJ, 701, 1357 Carroll S. M., Mantry S., Ramsey-Musolf M. J., Stubbs C. W., 2009, Phys. Rev. Lett., 103, 011301 Carroll S. M., Mantry S., Ramsey-Musolf M. J., 2010, Phys. Rev. D, 81, 063507 Cisterna A., Delsate T., Ducobu L., Rinaldi M., 2016, Phys. Rev. D, 93, 084046 Creminelli P., Vernizzi F., 2017, Phys. Rev. Lett., 119, 251302 De Martino I., della Monica R., De Laurentis M., 2021, Phys. Rev. D, 104, L101502 Della Monica R., de Martino I., 2022, J. Cosmology Astropart. Phys., 2022, 007 Della Monica R., de Martino I., de Laurentis M., 2022, MNRAS, 510, 4757 Do T. et al., 2019, Science, 365, 664 Doeleman S. S. et al., 2008, Nature, 455, 78 Dormand J., Prince P., 1980, Journal of Computational and Applied Mathematics, 6, 19 Dvali G., Gabadadze G., Porrati M., 2000, Physics Letters B, 484, 112 Dyadina P. I., Avdeev N. A., Alexeyev S. O., 2019, MNRAS, 483, 947 Eisenhauer F., Schödel R., Genzel R., Ott T., Tecza M., Abuter R., Eckart A., Alexander T., 2003, ApJ, 597, L121 Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306 Genzel R., Eisenhauer F., Gillessen S., 2010, Reviews of Modern Physics, 82, 3121 Ghez A. M., Klein B. L., Morris M., Becklin E. E., 1998, ApJ, 509, 678 Ghez A. M. et al., 2003, ApJ, 586, L127 Gillessen S. et al., 2017, ApJ, 837, 30 Goodman J., Weare J., 2010, Communications in Applied Mathematics and Computational Science, 5, 65 Gravity Collaboration, 2017, A&A, 602, A94 Gravity Collaboration, 2018, A&A, 615, L15 Gravity Collaboration, 2020, A&A, 636, L5 Green R. M., 1985, Spherical Astronomy, Cambridge University Press Grould M., Vincent F. H., Paumard T., Perrin G., 2017, A&A, 608, A60 Hees A. et al., 2017, Phys. Rev. Lett., 118, 211101 Hohmann M., 2015, Phys. Rev. D, 92, 064019 Horndeski G. W., 1974, Int. J. Theor. Phys., 10, 363 Hou S., Gong Y., 2018, Eur. Phys. J. C, 78, 247 Kaplan D. B., Wise M. B., 2000, Journal of High Energy Physics, 2000, 037 Kobayashi T., 2019, Rep. Prog. Phys., 82, 086901 Kobayashi T., Yamaguchi M., Yokoyama J., 2011, Prog. Theor. Phys., 126, 511 Maselli A., Silva H. O., Minamitsuji M., Berti E., 2015, Phys. Rev. D, 92, 104049 Maselli A., Silva H. O., Minamitsuji M., Berti E., 2016, Phys. Rev. D, 93, 124056 Plewa P. M. et al., 2015, MNRAS, 453, 3234 Schödel R. et al., 2002, Nature, 419, 694 Taff L. G., Szebehely V., 1986, Nature, 319, 630 S U P P O RT I N G I N F O R M AT I O N Supplementary data are available at MNRAS online. Figure S1 The figure illustrate the 1σ , 2σ , and 3σ contours of the full posterior distribution of the Newtonian limit of the Horndeski gravity. The parameters G4(0, 0) and G4(1, 0) are bound and the mean value of their distributions, plus the 68 per cent confidence interval are reported on the respective histograms. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article. A P P E N D I X A : N U M E R I C A L I N T E G R AT I O N O F T H E G E O D E S I C E Q UAT I O N S The gravitational field around a mass point source is described in Horndeski’s gravity at Newtonian order by the line element in equation (3). From it one can compute the geodesic equations that describe the motion of freely falling test particles in the Newtonian limit and for the four spherical isotropical coordinates: M t˙ṙ t¨ = − 1 r G4(0,0) ϒ + G24(1,0) e−mψ r + G24(1,0) mψ e−mψ r 8π G4(0,0) ϒr − G4(0,0) ϒM − G24(1,0) Me−mψ r t˙2 M r̈ = r φ̇ 2 sin2 θ + r θ̇ 2 + 16π G4(0,0) r   2 −mψ r 1 G4(1,0) e 1 × − − , mψ + r G4(0,0) ϒ r , (A1) (A2) 2ṙ θ̇ , (A3) r 2φ̇ ṙ 2φ̇ θ̇ cos θ − φ̈ = − , (A4) sin θ r where a dot represents a derivative with respect to an affine parameter τ. These equations can be integrated once initial conditions are assigned for the four space–time coordinates and the components of the initial tangent vector to the geodesic. In doing so, one has to take into account the normalization condition for the Gravity Collaboration, 2017, A&A, 602, A94 In particular, for a time-like geodesic, which describes the world-line of a massive particle, it is required that gμν ẋ μ ẋ ν = −1. This condition allows to assign the initial value for one of the ẋ μ components, given the other three. θ̈ = φ̇ 2 sin θ cos θ − MNRAS 519, 1981–1988 (2023) Downloaded from https://academic.oup.com/mnras/article/519/2/1981/6909075 by guest on 11 January 2023 RDM acknowledges support from Consejeria de Educación de la Junta de Castilla y León and from the Fondo Social Europeo. IDM acknowledges support from Ayuda IJCI2018-036198-I funded by MCIN/AEI/ 10.13039/501100011033 and: FSE ‘El FSE invierte en tu futuro’ o ‘Financiado por la Unión Europea ‘NextGenerationEU’/PRTR. IDM is also supported by the project PGC2018096038-B-I00 funded by the Spanish ‘Ministerio de Ciencia e Innovación’ and FEDER ‘A way of making Europe’, and by the project SA096P20 Junta de Castilla y León. DV and MDL acknowledge the support of Istituto Nazionale di Fisica Nucleare (INFN) iniziative specifiche QGSKY and TEONGRAV. DV also acknowledges the FCT project with ref. number PTDC/FIS-AST/0054/2021. 1987 1988 R. Della Monica et al. X = Bx + Gy, (A5) Y = Ax + F y, (A6) Z = Cx + Hy + D. (A7) where (x, y) are the coordinates of the star on its orbital plane, while (X, Y) and Z are the sky-projected position of the star and its distance from the observer, respectively. The constants A, B, C, F , G and H are called Thiele-Innes elements (Taff & Szebehely 1986) and are obtained from the inclination i of the orbit, the angle of the line of nodes and the argument of the pericentre ω, via A = cos cos ω − sin sin ω cos i, MNRAS 519, 1981–1988 (2023) (A8) B = sin cos ω + cos sin ω cos i, C = − sin ω sin i, (A9) (A10) F = − cos sin ω − sin cos ω cos i, (A11) G = − sin sin ω + cos cos ω cos i, (A12) H = − cos ω sin i. (A13) A similar expression leads from the orbital velocities (v x , v y ) to the sky-projected ones (VX , VZ ) and the radial velocity VZ : VX = Bvx + Gvy , (A14) VY = Avx + F vy , (A15) VZ = −(Cvx + Hvy ). (A16) Additionally, in order to properly reconstruct synthetic orbits for S2, we have to take into account on the projected quantities the classical Rømer delay and the redshift on the radial velocity coming from special relativity and gravitational time dilation (we refer to (Grould et al. 2017; Della Monica et al. 2022; De Martino et al. 2021) for a detailed analysis on how these effects can be quantified on our predicted orbits). Finally, we correct our synthetic orbits for systematic effects arising from the construction of the reference frame. These result in five additional parameters in the generation of a synthetic orbit of S2, namely, (x• , y• ) related to a zero-point offset of the central mass with respect to the origin of the astrometric reference frame; (v x, • , v y, • ) describing a drift over time of the central mass on the astrometric reference frame; v z, • related to an offset in the estimation of the radial velocity of the S2 star when correcting to LSR. This paper has been typeset from a TEX/LATEX file prepared by the author. Downloaded from https://academic.oup.com/mnras/article/519/2/1981/6909075 by guest on 11 January 2023 Furthermore, as it appears from equation (A3), if initial conditions are assigned so that θ (0) = π /2 and θ̇ = 0, a value θ̈ = 0 is obtained, identically. This implies that one can always define the reference frame so that the motion of the test particle occurs in the equatorial plane. Moreover, since the metric coefficients do not explicitly depend on the time coordinate t, the particular value of t(0) does not affect the integrated orbit. For the S2 star, we set for convenience t(0) = tA , where tA corresponds to the time of the last apocentre passage of the star, which occurred in ∼2010. Finally, as regards the initial condition on r, φ, and the corresponding components of the velocity, we set values that are consistent with the choice t(0). Namely, we set the radial and angular position of the star and its velocity at apocentre. These quantities are estimated by computing the Keplerian osculating ellipse at apocentre, as given by the Keplerian elements (Green 1985) measured for S2 from Gillessen et al. {2017) and reported in Table 1. We have integrated equations (A1–A4) numerically, by means of an adaptive step-size Runge–Gravity Collaboration, 2017, A&A, 602, A94Kutta algorithm (Dormand & Prince 1980), obtaining parametric arrays {t(τ ), r(τ ), θ(τ ), φ(τ )} and {t˙(τ ), ṙ (τ ), θ̇(τ ), φ̇(τ )} that describe the motion of S2 in a reference frame centered on the central source of the gravitational field. In order to be able to compare the integrated orbit with the observational data, a projection is required in the reference frame of a distant observer, by means of the following relations: