© ESO 2022

1. Introduction

The observed kinematical properties of the baryonic matter in galaxies suggest that galaxies are embedded within extended halos of dark matter whose distribution and contribution to the total mass of the galaxy largely vary with the size and the morphological type of the galaxy (Del Popolo 2014). This result consistently fits within the standard cosmological model, where the properties of the Universe on large scales, including the temperature anisotropies of the cosmic microwave background (CMB) (Planck Collaboration V 2020), the large-scale distribution of galaxies (Dodelson et al. 2002), and the abundance of light elements (Cyburt et al. 2016), imply that the matter content of the Universe is dominated by cold dark matter.

The standard cosmological model, however, suffers from some tensions both on large scales and on the scale of galaxies. On large scales, the discrepancy between the value of the Hubble constant, H₀, derived from observations in the local Universe and the value derived from the CMB (Verde et al. 2019), the unlikely features in the CMB temperature anisotropies (Schwarz et al. 2016), and the deficiency of ⁷Li (Mathews et al. 2020) have recently challenged the standard model. On the scale of galaxies, the tensions have been long-lasting: in disk galaxies, fine-tuned interactions between baryonic and dark matter are necessary to describe the radial acceleration relation (RAR) or the baryonic Tully-Fisher relation (McGaugh 2020); the missing satellite problem, the plane of satellite galaxies around the Milky Way, and the cusp-core problem in dwarf galaxies pose additional challenges (see e.g., Kroupa 2012; Del Popolo & Le Delliou 2017; de Martino et al. 2020, for reviews). Moreover, the search for the elementary particles constituting dark matter remains inconclusive (Tanabashi et al. 2018).

On the scale of galaxies, a number of models beyond the standard model have been suggested. Dark matter particles might not be cold and collisionless, but rather self-interacting (Tulin & Yu 2018), fuzzy (Hui et al. 2017), axion-like (Marsh 2016), or a superfluid (Ferreira et al. 2019). Alternatively, the theory of gravity might break down on galaxy and larger cosmic scales and thus dark matter might be unnecessary. Modified gravity theories proposed in the literature are numerous: for example, f(R) gravity, where the Ricci scalar R in the Einstein-Hilbert action is replaced by an arbitrary function of R (Nojiri & Odintsov 2011) or conformal gravity, where R is replaced by the contraction of the fourth-rank Weyl tensor (Mannheim 2019). However, both modified dark matter models and modified gravity theories have encountered a number of challenges. For example, neither f(R) nor conformal gravity appear to be able to yield the correct abundance of light elements (Azevedo & Avelino 2018; Elizondo & Yepes 1994). Conformal gravity also appears to suffer from severe fine-tuning in the description of the kinematics of disk galaxies (Campigotto et al. 2019) and requires an intergalactic medium in galaxy clusters that is too hot (Horne 2006; Diaferio & Ostorero 2009). MOdified Newtonian Dynamics (MOND) has been proved to be the most successful phenomenological model on the scale of galaxies (Sanders & McGaugh 2002; McGaugh 2020): attempts to build a covariant formulation of MOND have been challenging (Famaey & McGaugh 2012; Milgrom 2015; Hernandez et al. 2019), but recent developments appear to provide directions for viable solutions (Skordis & Złośnik 2019, 2021).

Matsakos & Diaferio (2016) suggested refracted gravity (RG), a theory of gravity where the Poisson equation is modified by the introduction of the gravitational permittivity, a monotonic increasing function of the local mass density. Refracted gravity appears to properly describe the observed phenomenology on the scale of galaxies and clusters of galaxies, with the role of dark matter mimicked by the gravitational permittivity (Matsakos & Diaferio 2016; Cesare et al. 2020). In addition, RG provides a straightfoward covariant extension within the family of the scalar-tensor theories (Sanna et al. 2021). In the covariant extension, the role of the gravitational permittivity is played by the scalar field. This same scalar field is also responsible for the accelerated expansion of the Universe. Thus, RG has the attractive feature of unifying the two dark sectors of the Universe into a single scalar field (Carneiro 2018; Bertacca et al. 2010). In addition, in its ultra-weak field limit, the covariant extension of RG seems to suggest the emergence of an acceleration scale below which Newtonian gravity breaks down. This acceleration scale, originally proposed by Milgrom (1983) in MOND, is strongly supported by the observed dynamics on the scale of galaxies (see e.g., de Martino et al. 2020; Merritt 2020, for reviews). However, the emergence of an acceleration scale in the covariant extension of RG still requires a solid confirmation.

In RG, the presence of regions where the local mass density is sufficiently low distinctively affects the gravitational field: in flat systems, such as disk galaxies, the gravitational permittivity bends the lines of the gravitational field towards the system, thus determining a boost of the gravitational field in the plane perpendicular to the minor axis of the system. In spherically symmetric systems, the field lines remain radial, but the gravitational field is enhanced by the inverse of the gravitational permittivity. In Cesare et al. (2020), we investigated the dynamics of 30 disk galaxies from the DiskMass Survey (DMS, Bershady et al. 2010). These galaxies are almost face-on, with an inclination angle with respect to the line-of-sight in the range 6 − 46 degrees. This inclination enables the measurements of both the rotation curve and the profile of the stellar velocity dispersion perpendicular to the plane of the disk, although the measures of the rotation velocities are generally less accurate than for edge-on galaxies. Cesare et al. (2020) proved that the bending of the gravitational field lines occurring in RG is sufficient to describe the kinematics of these 30 DMS galaxies without resorting to dark matter. This description depends on three parameters appearing in the gravitational permittivity. A single set of parameters describes the entire sample of 30 galaxies, supporting the expectation that these parameters should be universal. These encouraging results on flat systems motivated us to test whether RG is also able to describe the dynamics of spherical systems with the same set of parameters. In other words, we wish to probe whether the boost of the gravitational field can be set by the gravitational permittivity alone, independently of the redirection of the field lines.

Here, we illustrate the result of this test on three E0 galaxies, which are approximately spherical. We considered NGC 1407, NGC 4486, alias M87, and NGC 5846 from the SLUGGS survey¹, a spectrophotometric survey of more than 4000 extragalactic globular clusters (GCs) around 27 early-type galaxies (Pota et al. 2013; Brodie et al. 2014; Forbes et al. 2017). To our knowledge, SLUGGS is the only survey with publicly available photometric and kinematic information at large radii of early-type galaxies. The three galaxies we consider here are the only E0 galaxies in the SLUGGS survey. In elliptical galaxies, the baryonic matter within R_e, the effective radius at half projected luminosity, suffices to describe the galaxy dynamics with Newtonian gravity (e.g., Dekel et al. 2005; Mamon & Łokas 2005; Proctor et al. 2009; Cappellari et al. 2013). Thus, to probe RG, we need kinematic information in the outer regions, where Newtonian gravity breaks down unless dark matter is assumed to exist. X-ray emitting gas, planetary nebulae, and globular clusters (GCs) have been adopted as kinematic tracers in these regions, out to ∼10R_e, where the stellar luminosity fades out (Danziger 1997; Mathews & Brighenti 2003; Werner et al. 2019; Pota et al. 2013; Pulsoni et al. 2018). In addition, GCs usually appear separated into two distinct populations, according to their colour. The red sub-population of the GC system is generally more spatially concentrated and has a smaller velocity dispersion than the blue sub-population (Geisler et al. 1996; Ashman & Zepf 1998; Bassino et al. 2006; Brodie & Strader 2006; Faifer et al. 2011; Strader et al. 2011; Forbes et al. 2012; Lee et al. 2008; Pota et al. 2013). By adopting the SLUGGS survey and thus the GCs as tracers of the velocity field in the outer regions of ellipticals, we actually have two tracers rather than one and thus we can more strongly constrain the model. Therefore, the three E0 galaxies of our sample provide a suitable dataset for our test of RG.

Studying the outer regions of ellipticals might be, however, particularly insidious. Ellipticals tend to live in dense environments: NGC 1407 and NGC 5846 are within groups and NGC 4486 is the central galaxy of the Virgo cluster. Therefore, their outer regions are subject to the gravitational field of nearby galaxies and of the hosting system as a whole. Nevertheless, in our model below we neglect the effects of the environment and we model each galaxy as an isolated system. Despite this severe assumption, our RG model can still describe reasonably well the overall dynamics of these elliptical galaxies.

Section 2 summarises the main features of RG. Sections 3 and 4 describe the photometric and spectroscopic data that we used in our analysis. In Sect. 5, we illustrate our RG dynamical model, and in Sect. 6 we report our results. We conclude in Sect. 7. We adopt a Hubble constant H₀ = 73 km s⁻¹ Mpc⁻¹ (Riess et al. 2016) throughout the paper.

2. Refracted gravity

Refracted gravity is a classical theory of modified gravity where dark matter is mimicked by a gravitational permittivity, ϵ(ρ), a monotonic increasing function of the local mass density ρ. The gravitational permittivity appears in the modified Poisson equation

$$\begin{array}{c}\hfill \mathrm{\nabla }·[ϵ(\rho )\mathrm{\nabla }\varphi ]=4\pi G\rho ,\end{array}$$(1)

where ϕ is the gravitational potential and G the gravitational constant; the permittivity ϵ(ρ) has the asymptotic values

$$\begin{array}{c}\hfill ϵ(\rho )=\{\begin{array}{cc}1,\hfill & \rho \gg {\rho }_{\mathrm{c}}\hfill \\ {ϵ}_0,\hfill & \rho \ll {\rho }_{\mathrm{c}},\hfill \end{array}\end{array}$$(2)

where 0 < ϵ₀ ≤ 1 is the gravitational permittivity in vacuum and ρ_c is a critical density.

Adopting a permittivity depending on the mass density introduces a spatial scale over which we need to average the mass distribution. Here, we are interested in the dynamics of galaxies and this scale has to be sufficiently smaller than the size of a galaxy. Clearly, the identification of this scale remains an open issue that needs to be addressed when RG is developed beyond its current stage of a phenomenological model that is only applied to galactic and extragalactic dynamics (Matsakos & Diaferio 2016).

According to Eqs. (1) and (2), when the density is much larger than ρ_c, we recover the Newtonian Poisson equation ∇²ϕ = 4πGρ. On the contrary, in low-density regions, the gravitational permittivity enhances the gravitational field. Refracted gravity predicts a different behaviour for the gravitational field in spherical and flattened systems. For spherical systems, the integration of Eq. (1) yields

$$\begin{array}{c}\hfill \frac{\mathrm{d}\varphi }{\mathrm{d}r}=\frac{GM(r)}{ϵ(\rho )r^2},\end{array}$$(3)

where M(r) is the total mass of the system within radius r. The RG gravitational field has the same direction and r-dependence as in the Newtonian case, but it is enhanced by the inverse of the permittivity.

Equation (3) shows that the gravitational field is ∝r⁻² in the vacuum, where both M(r) and ϵ(ρ) are constant. This behaviour of the RG gravitational field does not necessarily contradict the observed velocity dispersion profiles of elliptical galaxies (Jiménez et al. 2013; Durazo et al. 2017, 2018) and of GCs (Hernandez & Jiménez 2012; Hernandez et al. 2013a; Durazo et al. 2017; Hernandez & Lara-D I 2020), or the external ρ ∝ r⁻³ mass density profiles of spherical stellar systems (Hernandez et al. 2013b) that suggest a gravitational field falling off as r⁻¹. No astrophysical system is completely isolated and exactly satisfies the conditions M(r) = const and ϵ(ρ) = const in its outer regions. Therefore, in principle, the permittivity ϵ(ρ) can determine an r⁻¹ behaviour of the gravitational field in each of these systems depending on the actual density field of its environment.

For flat systems, the lines of the gravitational field are focussed towards the mid-plane of the source (Matsakos & Diaferio 2016; Cesare et al. 2020): the flatter the system, the stronger the redirection effect. This feature explains the rotation curves of disk galaxies and might be consistent with the claimed correlation between the dark matter content of elliptical galaxies and their ellipticity (Deur 2014, 2020).

In this work, we solved Eq. (3) for three approximately spherical galaxies to determine whether the gravitational permittivity can describe the dynamics of spherical systems independently of the redirection of the field lines. In Eq. (3), the mass density ρ and, consequently, the cumulative mass profile M(r) of each galaxy are associated with the baryonic matter alone. For the gravitational permittivity, according to Matsakos & Diaferio (2016) and Cesare et al. (2020), we adopted the smooth step function

$$\begin{array}{c}\hfill ϵ(\rho )={ϵ}_0+(1-{ϵ}_0)\frac{1}{2}\{tanh\left[\text{ln}{\left(\frac{\rho }{{\rho }_{\mathrm{c}}}\right)}^Q\right]+1\}.\end{array}$$(4)

The parameter Q regulates the steepness of the transition between the Newtonian and the RG regimes: the larger Q, the steeper the transition. With this formulation, RG has three free parameters that we expect to be universal: ϵ₀, ρ_c, and Q.

3. Photometric data

In this section, we describe the observables that entered our mass model of each galaxy: the surface brightness of the stars, the surface number densities of the GCs populations, the mass density of the hot X-ray emitting gas, and the mass of the central supermassive black hole (SMBH). Table 1 lists the quantities characterising the three galaxies and the colour cut that we adopted in Sect. 4.2 to separate the GC samples into the red and blue populations.

The galaxies NGC 1407, NGC 4486, and NGC 5846 are at redshift z < 0.007. Thus, we could neglect any k-correction in the photometric measures and in the distance modulus m − M = 5log₁₀[D(pc)] − 5, with m and M the apparent and absolute magnitudes of the source. Similarly, to convert the radial coordinate R projected on the sky from angular units, in arcsec, to physical units, in kpc, we adopted the relation, valid in a non-expanding Euclidean geometry,

$$\begin{array}{c}\hfill R(\mathrm{kpc})=4.84814\times {10}^{-6}D(\mathrm{kpc})R(\mathrm{arcsec}),\end{array}$$(5)

with D the distance to the galaxy. For elliptical galaxies, the radial coordinate R is the circularised radius

$$\begin{array}{c}\hfill R=\sqrt{qx{{}^{\prime }}^2+\frac{y{{}^{\prime }}^2}{q}},\end{array}$$(6)

where q is the minor-to-major axis ratio of the galaxy on the sky, and the x′ and y′ coordinates are orientated along the galaxy major and minor axes, respectively (e.g., Cappellari 2008; Cappellari et al. 2013; Napolitano et al. 2009; Pota et al. 2013, 2015).

3.1. Stellar surface brightness profiles

For the surface brightness of the stars, we adopted the models derived by Pota et al. (2015) for NGC 1407 and by Scott et al. (2013) for NGC 4486 and NGC 5846.

3.1.1. NGC 1407

The stellar surface brightness profile of NGC 1407 was measured in the B-band with the Hubble Space telescope/ACS (Spolaor et al. 2008; Rusli et al. 2013) and in the g-band with the Subaru/Suprime-Cam (Pota et al. 2013). The B- and g-bands largely overlap; therefore, Pota et al. (2015) derived a unique surface brightness profile in the B-band by a proper transformation of the g-band data. NGC 1407 has ellipticity ε = 1 − q = 0.05. Pota et al. (2015) modelled the surface brightness profile with the Sérsic profile

$$\begin{array}{c}\hfill I_{\ast }(R)=I_{\mathrm{e}}exp\{-b_{n_{\mathrm{s}}}[{\left(\frac{R}{R_{\mathrm{e}}}\right)}^{\frac{1}{n_{\mathrm{s}}}}-1]\},\end{array}$$(7)

where

$$\begin{array}{c}\hfill b_{n_{\mathrm{s}}}\approx 2n_{\mathrm{s}}-\frac{1}{3}+\frac{4}{405n_{\mathrm{s}}}+\frac{46}{25515n_{\mathrm{s}}^2}\end{array}$$(8)

is such that the effective radius R_e contains half of the total luminosity of the stellar distribution (Ciotti & Bertin 1999), I_e is the surface brightness at radius R_e, and n_s is the Sérsic index. After deconvolving for the seeing, Pota et al. (2015) derived R_e = (100 ± 3) arcsec, I_e = 3.5 × 10⁵ L_⊙, B arcsec⁻², and n_s = 4.67 ± 0.15.

3.1.2. NGC 4486 and NGC 5846

Scott et al. (2013) measured the surface brightness of NGC 4486 and NGC 5846 in the ugriz bands from the Sloan Digital Sky Survey DR7 (Abazajian et al. 2009) and the Wide Field Camera on the 2.5-m Isaac Newton Telescope at the Roque de los Muchachos Observatory. Scott et al. (2013) modelled the two-dimensional map of the surface brightness in the r-band; this band reduces the dust contamination and provides images with the optimal signal-to-noise value (Cappellari et al. 2013). They adopted the axisymmetric Multi-Gaussian Expansion (MGE) approach (Emsellem et al. 1994), which yields the surface brightness map (Cappellari 2008; Cappellari et al. 2013)

$$\begin{array}{c}\hfill I_{\ast }(x{}^{\prime },y{}^{\prime })={\displaystyle \sum _{k=1}^N\frac{L_k}{2\pi {\sigma }_k^{2}q{{}^{\prime }}_k}exp[-\frac{1}{2{\sigma }_k^2}(x{{}^{\prime }}^2+\frac{y{{}^{\prime }}^2}{q_k^{\prime 2}})],}\end{array}$$(9)

where N is the number of Gaussian components in the MGE fit, and L_k, σ_k, and $0\le q_k^{\prime }\le 1$ are the luminosity, the standard deviation along the major axis, and the observed minor-to-major axis ratio of each Gaussian, respectively; x′ and y′ are the coordinates on the plane of the sky, where the x′-axis is orientated along the major axis of the galaxy. Scott et al. (2013) found $q_k^{\prime }=1$ within 5% for almost all the Gaussians of the model of each galaxy. Indeed, NGC 4486 has ellipticity ε = 1 − q = 0.14, and NGC 5846 has ε = 0.08. Therefore, we modelled these galaxies as spherical systems, set $q_k^{\prime }=1$ for all the Gaussians, and adopted the parameters L_k and σ_k determined by Scott et al. (2013) for Eq. (9) in the surface brightness profile (Cappellari 2008)

$$\begin{array}{c}\hfill I_{\ast }(R)={\displaystyle \sum _{k=1}^N\frac{L_k}{2\pi {\sigma }_k^2}exp(-\frac{R^2}{2{\sigma }_k^2}).}\end{array}$$(10)

Scott et al. (2013) found N = 9 and N = 7 Gaussian components for NGC 4486 and NGC 5846, respectively. Tables A.1 and A.2, list L_k and σ_k after deconvolving for the seeing. Finally, the total stellar luminosity of the galaxy is (Cappellari et al. 2013)

$$\begin{array}{c}\hfill L_{\ast ,\mathrm{tot}}={\displaystyle \sum _{k=1}^N{L}_k.}\end{array}$$(11)

3.1.3. Stellar 3D luminosity density profiles

For NGC 1407, the three-dimensional (3D) luminosity density of the stars is the Abel integral, valid in spherical symmetry,

$$\begin{array}{c}\hfill {\nu }_{\ast }(r)=-\frac{1}{\pi }{\displaystyle {\int }_r^{+\infty }\frac{\mathrm{d}{I}_{\ast }}{\mathrm{d}R}\frac{\mathrm{d}R}{\sqrt{R^2-r^2}},}\end{array}$$(12)

where I_*(R) is given by Eq. (7) and r is the radial coordinate in three dimensions. For NGC 4486 and NGC 5846, the deprojected MGE 3D luminosity density is (Cappellari 2008)

$$\begin{array}{c}\hfill {\nu }_{\ast }(r)={\displaystyle \sum _{k=1}^N\frac{L_k}{{(\sqrt{2\pi }{\sigma }_k)}^3}exp(-\frac{r^2}{2{\sigma }_k^2})·}\end{array}$$(13)

The 3D luminosity density of the stars, ν_*(r), yields the cumulative stellar luminosity profile

$$\begin{array}{c}\hfill L_{\ast }(r)=4\pi {\displaystyle {\int }_0^r{\nu }_{\ast }(r^{\prime })r^{\prime 2}\mathrm{d}{r}^{\prime }.}\end{array}$$(14)

Setting r = +∞ in Eq. (14) yields the total stellar luminosity L_*, tot = 8.53 × 10¹⁰ L_⊙ for NGC 1407 (Pota et al. 2015). Adopting Eq. (11) yields L_*, tot = 7.40 × 10¹⁰ L_⊙ and 4.62 × 10¹⁰ L_⊙ for NGC 4486 and NGC 5846, respectively (Scott et al. 2013).

3.2. Number density profiles of GCs

The number density profiles of GCs were measured by Pota et al. (2015), Strader et al. (2011), and Pota et al. (2013), for NGC 1407, NGC 4486, and NGC 5846, respectively. The profiles were estimated separately for blue and red GCs, according to the colour thresholds listed in Table 1. Below, we remind the criteria for the GC selection because we applied the same criteria in Sect. 4.2.1 to derive the kinematic profiles of the GCs in NGC 4486 and NGC 5846.

To ensure the magnitude completeness of the GCs samples in all the three galaxies, Pota et al. (2013, 2015) and Strader et al. (2011) removed the sources fainter than M_i = −8.0 in the i-band. This magnitude limit corresponds to the peak of the GC luminosity function, which has a roughly Gaussian shape (e.g., Harris et al. 1991; Whitmore et al. 1995; Kundu & Whitmore 1998). Moreover, the samples might be contaminated by ultracompact dwarf galaxies (UCDs), which are a kinematically and spatially distinct population from the GCs. Therefore, Pota et al. (2013, 2015) and Strader et al. (2011) removed the sources brighter than M_i = −11.6. This luminosity is 1 mag brighter than ω Cen, the brightest GC in the Milky Way, and brighter objects are likely to be UCDs. Imposing the magnitude range −11.6 < M_i < −8.0 removes 72 objects from the NGC 1407 sample, and 53 objects from the NGC 4486 sample, whereas the NGC 5846 sample remains unaffected.

In addition to this photometric selection, Pota et al. (2013, 2015) and Strader et al. (2011) applied a kinematic selection on the two separated populations of blue and red GCs. Despite Galactic stars and GCs are two populations whose velocity distributions are usually well distinct, the population of stars might overlap on the low-velocity tail of the GC distribution so that a Galactic star could erroneously be identified as a GC. Thus, Pota et al. (2013, 2015) and Strader et al. (2011) removed the sources whose velocity is more than 3σ discrepant from the mean velocity of the 20 closest neighbours in the GC sample, where σ is the velocity dispersion of these 20 neighbours (Merrett et al. 2003). This kinematic criterion removed six, five, and two objects from the NGC 1407, NGC 4486, and NGC 5846 samples, respectively. The sizes of the final samples of blue and red GCs for each galaxy are listed in Table 1.

Pota et al. (2013, 2015) and Strader et al. (2011) binned each GC sample into circularised annuli, according to Eq. (6), around the galaxy centre. The background-subtracted surface number density profiles, along with their Poissonian uncertainties, are shown in Fig. 1. These uncertainties do not include possible biases due to selection effects: GC counts can miss objects either in the central region of the galaxy, because of the high surface brightness of the host galaxy, or in the galaxy outskirts, because of the star-shaped footprint of the SLUGGS survey.

Fig. 1. Surface number density profiles of blue (blue symbols and blue lines) and red (red symbols and red lines) GCs. Blue and red dots with error bars show the binned profiles; the solid lines are the models.

3.2.1. NGC 5846

We modelled the binned profiles of the surface number density of the blue and red GCs in NGC 5846 with the Sérsic profile

$$\begin{array}{c}\hfill N_{\mathrm{GC}}(R)=N_{\mathrm{e}}exp\{-b_{n_{\mathrm{s}}}[{\left(\frac{R}{R_{\mathrm{e}}}\right)}^{\frac{1}{n_{\mathrm{s}}}}-1]\},\end{array}$$(15)

where the parameters R_e, N_e, and n_s have similar meaning of those in Eq. (7).

We estimated the free parameters of the surface number density profiles of the blue and red GCs with a Monte Carlo Markov chain (MCMC) method with a Metropolis–Hastings acceptance criterion. Details of this algorithm are in Sect. 5.2. We ran the MCMC for 2 × 10⁶ steps with 10⁵ burn-in elements to achieve the convergence of the chains, according to the Geweke (1992) diagnostics. We adopted flat priors on the free parameters of the model in the ranges listed in Table 2. The upper limit of R_e is the radius of the most external data point of the surface number density; the upper limit of N_e is slightly larger than the surface number density at the smallest R. Therefore, the upper limits of the priors for R_e and N_e are different for the two GC populations. The distribution of the GCs, particularly the red population, generally follows the distribution of the stars of the host galaxy (e.g., Brodie & Strader 2006). The Sérsic index of the stellar component in giant elliptical galaxies is ≳4 (Kormendy et al. 2009). We thus set to 10 the upper limit of the prior on the Sérsic index for the two GCs sub-populations.

Our posterior distributions are single-peaked. Thus, we adopted their medians as our parameter estimates and the range between the 15.9 and 84.1 percentiles, which includes 68% of the posterior cumulative distribution, as our parameter uncertainties. These parameters are listed in Table 3. They provide the curves shown in the right panel of Fig. 1.

3.2.2. NGC 1407 and NGC 4486

For NGC 1407 and NGC 4486, we adopted the models available in the literature. Pota et al. (2015) modelled their data of NGC 1407 with Eq. (15). For NGC 4486, Strader et al. (2011) adopted the slightly different parametrisation

$$\begin{array}{c}\hfill N_{\mathrm{GC}}(R)=N_{0}exp[-{\left(\frac{R}{R_{\mathrm{s}}}\right)}^{\frac{1}{m}}].\end{array}$$(16)

Pota et al. (2015) and Strader et al. (2011) derived the parameters listed in Table 3 that yield the curves shown in Fig. 1.

In Table 3, we also list the χ² of these models, for comparing the goodness of the fits of the six GC samples. The 3D number density of each GC sub-sample is

$$\begin{array}{c}\hfill {\nu }_{\mathrm{GC}}(r)=-\frac{1}{\pi }{\displaystyle {\int }_r^{+\infty }\frac{\mathrm{d}{N}_{\mathrm{GC}}}{\mathrm{d}R}\frac{\mathrm{d}R}{\sqrt{R^2-r^2}}·}\end{array}$$(17)

3.3. Mass density profiles of the X-ray gas

For the 3D mass density profile of the X-ray emitting gas, we adopted results from the literature for NGC 1407 (Zhang et al. 2007) and NGC 4486 (Fabricant et al. 1980). For NGC 5846, we derived our model of the density profile from the estimate of Paggi et al. (2017) of a discrete profile.

3.3.1. NGC 1407

Zhang et al. (2007) assumed spherical symmetry and modelled the 3D number density profile of the gas, n_g(r), with the two-β function

$$\begin{array}{c}\hfill n_{\mathrm{g}}^2(r)=n_{\mathrm{g},1}^2{[1+{\left(\frac{r}{R_{\mathrm{c},1}}\right)}^2]}^{-3{\beta }_1}+n_{\mathrm{g},2}^2{[1+{\left(\frac{r}{R_{\mathrm{c},2}}\right)}^2]}^{-3{\beta }_2}.\end{array}$$(18)

Zhang et al. (2007) derived the parameters n_g, 1 = 0.1 cm⁻³, n_g, 2 = 3.65 × 10⁻³ cm⁻³, R_c, 1 = (8.42 ± 0.70) arcsec, R_c, 2 = (58.3 ± 0.7) arcsec, β₁ = 0.70 ± 0.01, and β₂ = 0.45 ± 0.01 by fitting the X-ray surface brightness profile measured with Chandra ACIS in the 0.7–7 keV band and with ROSAT PSPC in the 0.2–2 keV band. Zhang et al. (2007) reported n_g, 1 and n_g, 2 without uncertainty.

The 3D mass density profile of the gas entering our model is

$$\begin{array}{c}\hfill {\rho }_{\mathrm{g}}(r)=\mu m_{\mathrm{H}}{n}_{\mathrm{g}}(r),\end{array}$$(19)

where m_H = 1.66054 × 10⁻²⁷ kg is the atomic unit mass and μ = 0.6 is the mean molecular weight, which assumes a completely ionised gas and a chemical composition with zero metallicity. The normalisations n_g, 1 and n_g, 2 yield ρ_g, 1 = 9.96 × 10⁻²⁶ g cm⁻³ and ρ_g, 2 = 3.64 × 10⁻²⁷ g cm⁻³.

3.3.2. NGC 4486

Fabricant et al. (1980) modelled the 3D mass density profile of the gas with the function

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_0}{{(1+b{}^{\prime }{r}^2+c{}^{\prime }{r}^4+d{}^{\prime }{r}^6)}^{n{}^{\prime }}},\end{array}$$(20)

obtained by deprojecting the model for the X-ray surface brightness measured in the 0.7–3.0 keV energy band with the Einstein Observatory. Fabricant et al. (1980) provided the model parameters without uncertainty: b′ = 9.724 × 10⁻¹ arcmin⁻², c′ = 3.810 × 10⁻³ arcmin⁻⁴, d′ = 2.753 × 10⁻⁸ arcmin⁻⁶, n′ = 0.59, and ρ₀ = 1.0 × 10⁻²⁵ g cm⁻³.

3.3.3. NGC 5846

We derived the continuous 3D gas mass density profile from the discrete mass density profile estimated by Paggi et al. (2017). Paggi et al. (2017) used the X-ray surface brightness measurements obtained with the Chandra ACIS and XMM-Newton MOS instruments. We modelled this discrete profile with the two-β function of Eq. (18), where the two normalisation constants are now ρ_g, 1 and ρ_g, 2. We estimated the parameters of the gas mass density with the MCMC method described in Sect. 5.2 adopting the Metropolis-Hasting acceptance criterion. We ran the chains for 2 × 10⁶ steps with 10⁵ burn-in elements. To test the chain convergence, we adopted the Geweke (1992) diagnostics. We adopted uniform priors in the intervals listed in Table 4. The upper limits of R_c, 1 and R_c, 2 are the radii of the most external data point of the measured gas mass density profile. The upper limits of ρ_g, 1 and ρ_g, 2 are slightly larger than the data point with the smallest R. Our posterior distributions are single-peaked, and thus we adopted the medians of the posterior distributions as our parameter estimates and the range between the 15.9 and 84.1 percentiles as our parameter uncertainties. The parameters of the model, along with their uncertainties, are listed in Table 4. Figure 2 shows the model superimposed on the data.

Fig. 2. Three-dimensional mass density profile of the hot gas of NGC 5846. The solid circles with error bars show the discrete profile derived by Paggi et al. (2017); the curve shows our model.

3.4. SMBH mass

We included the contribution of the central SMBH in our mass model. We modelled each SMBH as a point mass, adopting the SMBH masses provided by Rusli et al. (2013): $M_{\bullet }=4.5_{-0.4}^{+0.9}\times {10}^9{M}_{\odot }$, $6.2_{-0.5}^{+0.4}\times {10}^9{M}_{\odot }$, and $1.1_{-0.1}^{+0.1}\times {10}^9{M}_{\odot }$ for NGC 1407, NGC 4486, and NGC 5846, respectively.

4. Spectroscopic data

We used the photometric information described in the previous section to build our model of the mass distribution for each galaxy. We obtained our dynamical model by combining this information with the spectroscopic information on the velocity dispersion profiles of the stars (Sect. 4.1) and the GCs (Sect. 4.2).

4.1. Stellar velocity dispersion profiles

For NGC 1407, we used the velocity dispersion profile provided by Pota et al. (2015). For NGC 4486 and NGC 5846, we derived the velocity dispersion profiles from the two-dimensional data of the ATLAS^3D survey (Cappellari et al. 2011).

4.1.1. NGC 1407

Pota et al. (2015) derived the root-mean-square velocity dispersion profile from two different data sets, depending on the radial range. For radii in the range [0, 30] arcsec, they used the long-slit data, along the major axis of the galaxy, from the European Southern Observatory Faint Object Spectrograph and Camera (v.2) (EFOSC2) (Proctor et al. 2009). For radii in the range [40, 110] arcsec, they used the multislit Keck/DEIMOS data (Foster et al. 2016), again along the galaxy major axis. In the overlapping range [30, 40] arcsec, they kept both data sets.

The root-mean-square velocity dispersion profile of the stars at the circularised projected radius R is

$$\begin{array}{c}\hfill V_{\mathrm{rms}}(R)={[V_{\mathrm{rot}}^2(R)+{\sigma }^2(R)]}^{1/2},\end{array}$$(21)

where σ(R) is the stellar velocity dispersion and V_rot(R) = V(R)−V_sys is the major axis stellar rotation velocity in the galaxy frame, with V(R) the observed rotation velocity, and V_sys the galaxy systemic velocity. The analysis of Proctor et al. (2009) confirms that NGC 1407 is a slow rotator, namely V_rot ≪ σ; therefore, Eq. (21) yields V_rms(R) ≈ σ(R).

In our dynamical model described in Sect. 5, we assumed an orbital anisotropy parameter β independent of the 3D radius r, as Pota et al. (2015) did. However, by fitting Schwarzschild orbit-superposition models to their kinematic data, Thomas et al. (2014) inferred a variable β(r) within R = 2 arcsec of NGC 1407. Therefore, Pota et al. (2015) only considered the V_rms profile beyond R = 2 arcsec = 0.272 kpc. The top-left panel of Fig. 3 shows the V_rms(R) profile of the stars obtained by Pota et al. (2015).

Fig. 3. Root-mean-square velocity dispersion profiles of the stars (upper panels) and of the GCs (lower panels). Each column refers to the galaxy indicated in the upper panels. Upper middle and right panels: grey dots and black dots with error bars show the unbinned and binned profiles, respectively. Lower panels: red and blue dots with error bars refer to the blue and red GC populations, respectively. We note the different radial and velocity ranges of the panels.

4.1.2. NGC 4486 and NGC 5846

The data publicly available on the ATLAS^3D website² provide the two-dimensional map of V and σ of the galaxy on the plane of the sky, where V and σ are the mean line-of-sight velocity of the stars and their velocity dispersion. As V and σ, we adopted the values estimated from the simple Gaussian parametrisation, rather than from the Gauss-Hermite parametrisation, of the velocity distribution of the stars derived from the spectral line profile (Cappellari 2008; Cappellari et al. 2013). In NGC 4486 and NGC 5846, 99.6% and 97.9% of the data points, respectively, have V/σ ≤ 0.2. Therefore these galaxies, like NGC 1407, are slow rotators, and we have V_rms ≈ σ.

From the map, we iteratively removed the values of V_rms that deviate more than 3 standard deviations from the mean V_rms of the entire map (Cappellari et al. 2013). For NGC 4486 and NGC 5846, the procedure converges after seven and six iterations and removes 155 and 201 data points, respectively. This removal is necessary because some data points could be spurious, because of the presence of Milky Way stars or problematic bins at the edge of the field of view (Cappellari et al. 2013). To fold and average the data with respect to the galaxy centre and thus transform the two-dimensional map V_rms into a one-dimensional profile, V_rms(R), we rotated the galaxy map counter-clockwise from the celestial coordinates to the system of coordinates (x′,y′) of Eq. (6), such that x′ is aligned with the photometric major axis of the galaxy (Krajnović et al. 2011).

This procedure yields V_rms(R) profiles with 2444 and 1752 data points for NGC 4486 and NGC 5846, respectively. These profiles are shown in the upper middle and right panels of Fig. 3. These numbers are ∼2 orders of magnitude larger than the number of data points in the kinematic profiles of the GC populations of the two galaxies (see Sect. 4.2.1). As confirmed by preliminary tests, with these unbalanced data sets the dynamical models are driven by the V_rms(R) of the stars and are basically insensitive to the GC kinematic profiles. To overcome this problem, we binned the stellar kinematic data. Each bin contains a constant number of data points: N = 111 for NGC 4486 and N = 103 for NGC 5846. NGC 4486 and NGC 5846 have now 22 and 17 binned data, respectively. These numbers of binned data are now comparable to the numbers of bins adopted for the blue and red GCs profiles (see Sect. 4.2.1) and guarantee that the stars and GCs data have a comparable weight in constraining our dynamical model. The values associated to the 22 and 17 binned stellar data are the medians, both in the V_rms and R directions, of the values of the N data points in the bin. The uncertainties are the semi-interval between the 15.9 and the 84.1 percentiles of the distributions, which are nearly symmetric, of all the data points in the bin. The binned data are shown by the black dots with error bars in the upper middle and right panels of Fig. 3.

In general, massive early-type galaxies, such as these two galaxies, show negative orbital anisotropy parameters in their central regions, indicating tangential orbits; on the contrary, in the outer regions, the orbits appear radial (Thomas et al. 2014; Rantala et al. 2019). Numerical simulations show that this difference in the stellar orbits might originate if massive elliptical galaxies form from the merging of two progenitors with mass ratios larger than 1/3 (Rantala et al. 2019): at the centre of the new-born galaxy, the gravitational torques caused by the merging of the SMBHs of the two parent galaxies cause a reversal of the orbit behaviour, from radial to tangential. Therefore, because we assumed below an orbital anisotropy parameter β independent of radius r, we excluded the innermost 2 arcsec in the kinematic profiles of NGC 4486 and NGC 5846, consistently with the prescription adopted for NGC 1407 by Pota et al. (2015) (Sect. 4.1.1). This projected radius of 2 arcsec corresponds to 0.167 kpc and 0.235 kpc for NGC 4486 and NGC 5846, respectively, and thus is comparable to the radius of 0.272 kpc of NGC 1407.

4.2. Velocity dispersion profiles of the GCs

For NGC 4486 and NGC 5846, the velocity dispersion profiles of the two populations of GCs are not available in the literature. Therefore, we derived these profiles from the GC catalogues of Strader et al. (2011) and Pota et al. (2013) for NGC 4486 and NGC 5846, respectively. For the velocity dispersions of the GCs in NGC 1407, we adopted the profiles estimated by Pota et al. (2015).

4.2.1. NGC 4486 and NGC 5846

To derive the V_rms(R) profiles of the GCs in these two galaxies, we separated the GC sample into the blue and red sub-populations, according to the colour thresholds listed in Table 1. To each sub-sample, we applied the photometric and kinematic selection criteria described in Sect. 3.2 to recover the same sub-samples used to derive the number density profiles.

We binned each GC population into circularised annuli centred on the galaxy centre. Each annulus contains the same number of GCs. Specifically, each bin contains N = 30 and N = 25 GCs for the blue and red populations in NGC 4486 and N = 26 GCs for both populations in NGC 5846. These numbers are a trade-off between bins that are too poor and profiles that are too coarse. For NGC 4486, we obtained 16 and 8 bins for the blue and the red GCs, respectively. For NGC 5846, we obtained six and seven bins for the blue and the red GCs, respectively. In NGC 5846, some GCs are common to contiguous bins, so that we have a sufficient number of bins despite the small total number of GCs. The radius R of each bin is the median of the radial coordinates of the GCs in the bin.

The V_rms of the GCs in each bin is (Pota et al. 2013, 2015)

$$\begin{array}{c}\hfill V_{\mathrm{rms}}^2=\frac{1}{N}{\displaystyle \sum _{i=1}^N[{(V_{\mathrm{rad},i}-V_{\mathrm{sys}})}^2-{(\mathrm{\Delta }{V}_{\mathrm{rad},i})}^2],}\end{array}$$(22)

where V_rad, i is the radial velocity of the i-th GC in the radial bin and ΔV_rad, i is its uncertainty. We took V_rad and ΔV_rad from Strader et al. (2011), for NGC 4486, and Pota et al. (2013), for NGC 5846. To estimate the errors on each value of V_rms, we used the procedure reported in Danese et al. (1980), as performed by Pota et al. (2015). For convenience, we report this procedure in Appendix B.

The kinematic profiles of the GCs are shown in Fig. 3. These profiles are ∼48 and ∼13 times more extended than the kinematic profiles of the stars for NGC 4486 and NGC 5846, respectively. In NGC 4486, the V_rms(R) of the GC populations show a bell-like shape that peaks at large radii, specifically at ∼500 arcsec = 42 kpc. This feature suggests that the GC kinematics is influenced by the gravitational potential well of the Virgo cluster whose central giant elliptical galaxy is exactly NGC 4486 (Strader et al. 2011). This shape is common to the stellar velocity dispersion profiles of other elliptical galaxies located at the centre of galaxy clusters, such as the brightest cluster galaxy of Abell 383 (Geller et al. 2014).

4.2.2. NGC 1407

Pota et al. (2015) derived the GC velocity dispersion profiles with the same procedure that we adopted here for NGC 4486 and NGC 5846. Pota et al. (2015) obtained their kinematic data from the Keck/DEep Imaging Multi-Object Spectrograph (DEIMOS) instrument. The V_rms(R) profiles of the two GC samples are based on 153 blue GCs and 148 red GCs (see Table 1), and they are shown in Fig. 3. They are ∼5 times more extended than the star kinematic profile.

5. Dynamical model

We now describe our dynamical model (Sect. 5.1) and the MCMC approach we adopted to estimate our model parameters (Sect. 5.2).

5.1. Basic equations

We modelled the kinematics of the three E0 galaxies by assuming spherical symmetry and no net rotation, as suggested by the observed slow rotation of the galaxies. We modelled the velocity dispersion profile of each dynamical tracer t = { * ,R, B}, namely stars, red GCs, and blue GCs, as

$$\begin{array}{c}\hfill V_{\mathrm{rms},\mathrm{t}}^2(R)=\frac{2}{I_{\mathrm{t}}(R)}{\displaystyle {\int }_R^{+\infty }K({\beta }_{\mathrm{t}},\frac{r}{R}){\nu }_{\mathrm{t}}(r)\frac{\mathrm{d}\varphi }{\mathrm{d}r}r\mathrm{d}r,}\end{array}$$(23)

which is the solution to the spherical Jeans equations (Jeans 1915; Cappellari 2008; Mamon & Łokas 2005; Pota et al. 2015). In Eq. (23), R is the circularised radius projected onto the sky according to Eq. (6), whereas r is the 3D radius; I_t(R) is the surface brightness of the stars or the surface number density of GCs; ν_t(r) is the 3D luminosity density of the stars or the 3D number density of GCs; ϕ(r) is the gravitational potential; and ${\beta }_{\text{t}}=1-{\sigma }_{\theta }^2/{\sigma }_r^2$ is the orbital anisotropy parameter, with σ_θ and σ_r the velocity dispersions in the tangential and the radial directions, respectively. We assumed β_t to be independent of r, as discussed at the end of Sects. 4.1.1 and 4.1.2. K is the kernel

$$\begin{array}{cc}\hfill K({\beta }_{\mathrm{t}},\frac{r}{R})& =\frac{1}{2}{\left(\frac{r}{R}\right)}^{2{\beta }_{\mathrm{t}}-1}[(\frac{3}{2}-{\beta }_{\mathrm{t}})\sqrt{\pi }\frac{\mathrm{\Gamma }({\beta }_{\mathrm{t}}-\frac{1}{2})}{\mathrm{\Gamma }({\beta }_{\mathrm{t}})}\hfill \\ \hfill & +{\beta }_{\mathrm{t}}{B}_{\frac{R^2}{r^2}}({\beta }_{\mathrm{t}}+\frac{1}{2},\frac{1}{2})-B_{\frac{R^2}{r^2}}({\beta }_{\mathrm{t}}-\frac{1}{2},\frac{1}{2})],\hfill \end{array}$$(24)

(see Eq. (A16) in Mamon & Łokas 2005) with Γ(z) the Euler Γ function and B_x(a, b) the incomplete beta function. By inserting the RG gravitational field, Eq. (3), into Eq. (23), we obtain

$$\begin{array}{c}\hfill V_{\mathrm{rms},\mathrm{t}}^2(R)=\frac{2G}{I_{\mathrm{t}}(R)}{\displaystyle {\int }_R^{+\infty }K({\beta }_{\mathrm{t}},\frac{r}{R}){\nu }_{\mathrm{t}}(r)\frac{M(r)}{ϵ(\rho )}\frac{\mathrm{d}r}{r}.}\end{array}$$(25)

In our model, M(r) is the baryonic mass alone, namely

$$\begin{array}{c}\hfill M(r)=M_{\ast }(r)+M_{\mathrm{g}}(r)+M_{\bullet }(r),\end{array}$$(26)

where M_*(r), M_g(r), and M_•(r) are the cumulative mass profiles of the stars, X-ray emitting gas, and SMBH, respectively. We verified that the contribution of the GCs to the total galaxy mass profile is smaller than 0.5% at all radii (upper panels of Fig. 4) and that the maximal contribution of the GCs to the total mass density profile is ∼1% (lower panels of Fig. 4). We therefore ignored the contribution of the GCs to the mass profile of Eq. (26). We also ignored the GC contribution to the density profile

$$\begin{array}{c}\hfill \rho (r)={\rho }_{\ast }(r)+{\rho }_{\mathrm{g}}(r)+{\rho }_{\bullet }(r),\end{array}$$(27)

Fig. 4. Cumulative mass (upper panels) and mass density (lower panels) profiles of the baryonic components of each galaxy: stars (black solid line), blue GCs (blue solid line), red GCs (red solid line), and hot X-ray emitting gas (magenta solid line); the green solid lines in the upper panels show the SMBH contributions. The purple solid lines show the sum of all the baryonic contributions. The stellar profiles assume Υ = 7.6 M⊙/L⊙ for NGC 1407, and Υ = 3.6 M⊙/L⊙ for NGC 4486 and NGC 5846. The grey shaded area around the stellar profiles shows the mass variation by adopting mass-to-light ratios in the range [4.0, 11.2] M⊙/L⊙, for NGC 1407, and [1.7, 5.5] M⊙/L⊙, for NGC 4486 and NGC 5846, in the B- and r-band, respectively. The profiles of blue and red GCs assume a GC mass MGC = 105 M⊙; the blue and red shaded areas show the mass variation by adopting MGC in the range [104, 106] M⊙. The magenta shaded areas show the uncertainties on the adopted gas mass density profiles; the gas profiles of NGC 4486 do not have a shaded area because the uncertainties on the profiles are unavailable. The green solid lines and shaded areas in the upper panels are the masses of the SMBH and their uncertainties. The purple shaded areas show the possible range of the total baryonic profiles. Each column refers to the galaxy indicated in the upper panels.

which appears as the argument of the permittivity ϵ(ρ) in the velocity dispersion profile of Eq. (25).

We estimated the mass profile of the stars as

$$\begin{array}{c}\hfill M_{\ast }(r)=\mathrm{{\rm Y}}{L}_{\ast }(r),\end{array}$$(28)

where L_*(r) is the cumulative luminosity profile, Eq. (14), and Υ is the stellar mass-to-light ratio that we assumed to be independent of r. In our dynamical model, Υ is a free parameter whose prior is described in Sect. 5.2.

The mass profile of the gas, M_g(r), derives from the integration of its density profile. This density profile is indicated in Eq. (18), for NGC 1407 and NGC 5846, and in Eq. (20) for NGC 4486. For each galaxy, the mass of the SMBH, M_•, is the value reported in Sect. 3.4.

Figure 4 shows the mass and the mass density profiles of each component. The stars provide the largest contribution to the total mass of each galaxy over most radial range, except in the very centre, where the SMBH dominates, and in the outskirts of NGC 4486 and NGC 5846, where the gas contribution overcomes the star contribution.

5.2. The MCMC approach

For each galaxy, the dynamical model has seven free parameters: four parameters, Υ, ϵ₀, Q, and ρ_c, contribute to the gravitational potential well of the galaxy, and three parameters, β_*, β_B, and β_R, are specific for each tracer. We note that the contributions of the X-ray emitting gas and the SMBH to the gravitational potential have no free parameters. For convenience, hereafter, we use the parameters ℬ_* = −log₁₀(1 − β_*), ℬ_B = −log₁₀(1 − β_B), ℬ_R = −log₁₀(1 − β_R), and 𝒫_c = log₁₀(ρ_c [g cm⁻³]). Tangential and radial orbits correspond to ℬ_t < 0 and ℬ_t > 0, respectively, where t = ( * ,B, R).

We explored this seven-dimensional parameter space with a MCMC algorithm, where we adopted a Metropolis-Hastings acceptance criterion: the random variate x_i + 1 = (Υ, ϵ₀, Q, 𝒫_c, ℬ_*, ℬ_B, ℬ_R), at step i + 1 of the chain, is obtained from the probability density G(x|x_i), which depends on the random variate x_i at the previous step. For G(x|x_i), we adopted a multi-variate Gaussian density distribution with mean value x_i. As the standard deviations of this Gaussian distribution, we chose 5% of the ranges of the uniform priors that we adopted for the parameters of the model.

For each single tracer t = ( * ,B, R), we adopted the likelihood function

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{t}}(\mathit{x})={\displaystyle \left[\prod _{i=1}^{N_{\mathrm{t}}}\frac{1}{\sqrt{2\pi \mathrm{\Delta }{V}_{\mathrm{rms},\mathrm{data},\mathrm{t}}^2(R_i)}}\right]exp[-\frac{{\chi }_{\mathrm{t}}^2(\mathit{x})}{2}],}\end{array}$$(29)

where

$$\begin{array}{c}\hfill {\chi }_{\mathrm{t}}^2(\mathit{x})={\displaystyle \sum _{i=1}^{N_{\mathrm{t}}}\frac{{[V_{\mathrm{rms},\mathrm{mod},\mathrm{t}}(R_i,\mathit{x})-V_{\mathrm{rms},\mathrm{data},\mathrm{t}}(R_i)]}^2}{\mathrm{\Delta }{V}_{\mathrm{rms},\mathrm{data},\mathrm{t}}^2(R_i)},}\end{array}$$(30)

N_t is the number of data points of the tracer, V_{rms, data, t}(R_i) and ΔV_{rms, data, t}(R_i) are the measured root-mean-square velocity dispersions at the projected distances R_i and their uncertainties, and V_{rms, mod, t}(R_i) is the expected root-mean-square velocity dispersion derived from the integral in Eq. (25) at the same radius. We estimated this integral by adopting two independent linear grids in the R and the r coordinates. Both grids cover the range [10⁻⁴, 4800] arcsec, with a step of 4 arcsec. This grid size guarantees that our model includes more than 99.4% of the total stellar luminosity of each galaxy.

To consider the three tracers at the same time, we adopted the joint likelihood function

$$\begin{array}{c}\hfill \mathcal{L}(\mathit{x})={\displaystyle \prod _{\mathrm{t}=1}^3{\mathcal{L}}_{\mathrm{t}}(\mathit{x})=\frac{1}{\mathcal{N}}exp[-\frac{{\chi }_{\mathrm{tot}}^2(\mathit{x})}{2}],}\end{array}$$(31)

with 𝒩 the resulting normalisation factor, and

$$\begin{array}{c}\hfill {\chi }_{\mathrm{tot}}^2(\mathit{x})={\chi }_{\ast }^2(\mathit{x})+{\chi }_{\mathrm{B}}^2(\mathit{x})+{\chi }_{\mathrm{R}}^2(\mathit{x}).\end{array}$$(32)

The Metropolis–Hastings ratio is

$$\begin{array}{c}\hfill A=\frac{p(\mathit{x})\times \mathcal{L}(\mathit{x})}{p({\mathit{x}}_i)\times \mathcal{L}({\mathit{x}}_i)}\frac{G(\mathit{x}|{\mathit{x}}_i)}{G({\mathit{x}}_i|\mathit{x})},\end{array}$$(33)

where p(x) is the product of the priors of the components of x. In the case of uniform priors, as we adopted in our analysis, p(x) = const. If A ≥ 1, we accept the proposed combination of free parameters and thus we set x_i + 1 = x; otherwise, we either set x_i + 1 = x, with probability A, or x_i + 1 = x_i, with probability 1 − A. For each galaxy, we ran two Markov chains with two different random seeds. Each chain has 8 × 10⁴ steps, where the first 6 × 10³ steps are discarded as the burn-in chain. We assessed the convergence of the chains with the Gelman & Rubin (1992) diagnostics, and the convergence of each chain was further confirmed according to the Geweke (1992) diagnostics. In the following, we consider the posterior distributions of the parameters from the two chains combined.

To model each galaxy with our MCMC analysis, we need to adopt a prior for each of our seven free parameters. For NGC 1407, the photometric information available in the B-band suggests the uniform prior [4.0, 11.2] M_⊙/L_⊙ for Υ. For NGC 4486 and NGC 5846, with photometric information in the r-band, we adopted the uniform prior [1.7, 5.5] M_⊙/L_⊙. We chose both prior intervals according to the ranges expected from the SPS models of Humphrey et al. (2006) and Zhang et al. (2007), for the B-band, and of Bell et al. (2003) and Zibetti et al. (2009), for the r-band. These ranges depend on the initial mass function: for the B-band, the lower and upper limits are set by the Kroupa (2001) and Salpeter (1955) initial mass functions, respectively, and for the r-band, the lower and upper limits are set by the Bottema (1993) and Salpeter (1955) initial mass functions, respectively. For all the three orbital anisotropy parameters, ℬ_*, ℬ_B, and ℬ_R, we adopted the uniform prior [ − 1.5, 1.0], thus including almost all possible orbits, from very tangential to very radial. For NGC 1407, our priors coincide with those adopted by Pota et al. (2015), who modelled the kinematics of this galaxy in Newtonian gravity with a generalised Navarro-Frenk-White dark matter halo (Hernquist 1990; Zhao 1996).

For the three permittivity parameters, we adopted uniform priors in the ranges (0, 1], [0.01, 2], and [ − 27, −23], for ϵ₀, Q, and 𝒫_c, respectively. These priors are those adopted in Cesare et al. (2020), except for ϵ₀: our range here is wider, compared to the range [0.1, 1] of Cesare et al. (2020). Our range for 𝒫_c includes the extreme values adopted by Matsakos & Diaferio (2016). This range includes the typical galaxy densities where Newtonian gravity does not suffice to describe the kinematic data. In passing, the range [10⁻²⁷, 10⁻²³] g cm⁻³ includes the lowest density of the interstellar medium (∼10⁻²⁵ g cm⁻³, Wood & Linsky 1997; Thornton et al. 1998), but it is well above the average density of the intergalactic medium (∼10⁻³¹ g cm⁻³, Meiksin 2009). The smallest pressure reached in ultra-high vacuum chambers on Earth is ∼10⁻¹⁰ mbar, corresponding to a density ∼10⁻¹⁷ − 10⁻¹⁶ g cm⁻³ for standard air with molar mass 29 g mol⁻¹ (Cooper et al. 2019; Liu et al. 2021). The priors on the free parameters of the model are summarised in Table 5.

6. Results

We now illustrate how our model describes the velocity dispersion profiles of the three galaxies (Sect. 6.1) and discuss our estimates of the mass-to-light ratios and the anisotropy parameters (Sect. 6.2) and of the parameters of the RG permittivity (Sect. 6.3).

6.1. Velocity dispersion profiles

Figures 5–7 show the posterior distributions of the parameters of the model obtained from the two joint Markov chains for each galaxy. Figure 8 compares the kinematic data with the models of the V_rms profiles whose parameters are the medians of the posterior distributions (solid lines). These parameters are listed in Table 5. The shape with a single peak of most of the posterior distributions justifies the use of the medians as estimators of the free parameters.

Fig. 5. Posterior distributions of the seven parameters derived from the two joint Markov chains for NGC 1407. The parameters are the stellar mass-to-light ratio, Υ, the velocity anisotropy parameters of the three tracers, ℬ*, ℬB, ℬR, and the three permittivity parameters, ϵ0, Q, and 𝒫c. The red squares with error bars show the medians and their uncertainites, set by the 15.9 and 84.1 percentiles of the posterior distributions. The yellow, cyan, and blue contours limit the 1σ, 2σ, and 3σ regions, respectively. The medians and their uncertainites are also reported above each column and in the one-dimensional posterior distributions in the top panels of each column as red solid and black dashed lines, respectively.

Fig. 6. Same as Fig. 5, but for NGC 4486.

Fig. 7. Same as Fig. 5, but for NGC 5846.

Figure 8 shows that, globally, the RG models provide an adequate description of the kinematic data of the three tracers for all the three galaxies. The reduced χ²’s listed in Table 5 support this conclusion. Nevertheless, in NGC 5846, the model V_rms of the blue GCs is systematically smaller than the data; similarly, in NGC 4486, the model is unable to describe the peak, at R ∼ 40 kpc, of the V_rms profile of the blue GCs. This poor modelling might originate from the assumption that these two galaxies are isolated rather than embedded in a larger system, as we discuss below. These discrepancies might also be reduced if we remove our assumption that the velocity anisotropy parameters β are independent of radius r or that the galaxies have null net rotation: albeit weak, the net galaxy rotation is not completely absent in the real galaxies.

Fig. 8. Models of the root-mean-square velocity dispersion profiles emerging from our MCMC analysis in RG (solid lines). Upper panels: the green, blue, and red solid lines show the model profiles for the stars, blue GCs, and red GCs, respectively. Lower panels: zoom-in of the stellar profiles. The solid circles with error bars are the measured profiles from Fig. 3. The green, blue, and red dashed lines show the kinematic profiles of the stars, blue GCs, and red GCs computed with Eq. (23) where the values of the mass-to-light ratio Υ and the velocity anisotropy parameters β are those of the RG model but the gravitational field is Newtonian with no dark matter. The discrepancy with the RG curves highlights the gravitational boosting generated by the RG gravitational permittivity.

6.2. Mass-to-light ratios and anisotropy parameters

We set the ranges of the uniform prior for the mass-to-light ratios according to the SPS models. For NGC 1407, the posterior distribution of the mass-to-light ratio peaks around the centre of this range. For NGC 4486 and NGC 5846, the posterior distribution of the mass-to-light ratio peaks close to the upper limit of the range of the prior, but the probability that the actual mass-to-light ratio lies in the rest of the range is substantial. Thus, we conclude that RG can indeed describe the data with stellar mass-to-light ratios consistent with the SPS models.

For NGC 1407, the values of Υ, ℬ_*, ℬ_B, and ℬ_R that we obtained are at −0.38σ, +1.00σ, −0.21σ, and +0.097σ, respectively, from the parameters found by Pota et al. (2015) in their analysis with Newtonian gravity and a dark matter halo. The RG models require radial orbits for the stars and tangential orbits for the blue GCs in all the three galaxies. The red GCs have tangential orbits in NGC 4486 and NGC 5846. In NGC 1407, the red GCs have radial orbits, in agreement with the result found by Pota et al. (2015) in Newtonian gravity.

6.3. Parameters of the RG permittivity

The orange squares with error bars in Fig. 9 show the permittivity parameters estimated for the three E0 galaxies: the parameters are consistent with each other within 2σ, suggesting their universality. Their mean values are $\{{ϵ}_0,Q,{\mathcal{P}}_{\mathrm{c}}\}=\{0.{089}_{-0.035}^{+0.038},0.{47}_{-0.21}^{+0.29},-24.{25}_{-0.20}^{+0.28}\}$. Figure 9 also shows the permittivity parameters derived from 30 disk galaxies of the DMS sample (Cesare et al. 2020): the purple squares with error bars show the means of the parameters estimated for each individual DMS galaxy, whereas the green dots at the centre of the shaded areas show the global values derived with an approximate procedure from the entire sample of 30 galaxies.

Fig. 9. Permittivity parameters estimated from the three E0 galaxies in our sample (orange squares with error bars) compared with the permittivity parameters estimated by Cesare et al. (2020) from the DMS disk galaxies. The purple squares with error bars show the means of the permittivity parameters found for the individual DMS galaxies. The light blue shaded areas show the posterior distributions of the three permittivity parameters found with an approximate procedure from the entire DMS sample, with the green dots indicating their median values and the yellow, red, and black contours indicating the 1σ, 2σ, and 3σ levels, respectively.

Our estimates of 𝒫_c and Q are within 1σ from the mean values derived from the individual DMS galaxies, whereas the mean DMS permittivity of the vacuum ϵ₀ is within 3σ from our three estimates. This marginal discrepancy might originate from our simplistic modelling of the E0 galaxies. Indeed, here we assume that these galaxies do not have net rotation and are relaxed and isolated, whereas, in fact, these galaxies are in clusters or groups and show signs of interactions with nearby galaxies: NGC 1407 is at the centre of the Eridanus A group (Brough et al. 2006) and certainly feels the gravitational field of its neighbouring galaxy NGC 1400. NGC 4486, namely M87, is the central giant elliptical galaxy of the Virgo cluster (Strader et al. 2011) and NGC 5846 is the central and brightest galaxy of a galaxy group (Mahdavi et al. 2005). Moreover, the GC sample of NGC 5846 might be contaminated by the GCs of its neighbour NGC 5846A (Pota et al. 2013). The environment is expected to have a relevant effect on the intensity of the gravitational field in RG (Matsakos & Diaferio 2016); this effect clearly propagates into the actual value of the permittivity of the vacuum ϵ₀. Environmental effects can also be responsible for the poor modelling of the velocity dispersion profiles of the blue GCs of NGC 4486 and NGC 5846 shown in Fig. 8.

Assessing the gravitational interactions of neighbouring galaxies in RG is difficult and is beyond the purpose of our present work. The shape of the mass density distribution of an individual source determines the refraction angle of the lines of the gravitational field and thus how the direction and the intensity of the gravitational field depend on the vector r of the distance from the source (Matsakos & Diaferio 2016). Therefore, estimating the gravitational interactions among galaxies with arbitrary mass density distributions requires an accurate integration of the field equations with little room for simplyfing assumptions.

The tension between our vacuum permittivity ϵ₀ and the vacuum permittivity derived from the DMS sample might become more severe if we consider the DMS global value, namely the green dot in Fig. 9. Our ϵ₀ shows a 14.8σ tension from the DMS parameter, which is clearly driven by the small width of the posterior distribution of the DMS global value. On the contrary, the parameter Q is within 3.4σ from the DMS global parameter. More importantly, the critical density 𝒫_c, which represents the most relevant parameter of RG, is within 1.4σ from the DMS global parameter. The large discrepancy for ϵ₀ might be due to the approximate procedure adopted by Cesare et al. (2020) to estimate this DMS global value: with their MCMC analysis, Cesare et al. (2020) only explored the space of the three permittivity parameters while keeping fixed the mass-to-light ratio Υ and the disk-scale height h_z of each of the 30 individual galaxies. This approach is a shortcut to the appropriate, but computationally overwhelming, procedure of exploring the 63-dimensional parameter space of the full sample, namely the two parameters, Υ and h_z, of each of the 30 galaxies and the three permittivity parameters. The set of the permittivity global parameters returned by this shortcut still describes the kinematic profiles of each galaxy, but the agreement is poorer than the agreement obtained by modelling each galaxy individually. Therefore, these values of the permittivity parameters should be considered with caution.

7. Discussion and conclusion

We used the kinematic information of the stellar population and of the GCs in the outer regions of three E0 galaxies from the SLUGGS survey (Pota et al. 2013; Brodie et al. 2014; Forbes et al. 2017) to show that RG can describe the kinematics of spherical systems without requiring the existence of dark matter. Our results complement the ability of RG to describe the rotation curves and the vertical velocity dispersion profiles of 30 disk galaxies from the DMS (Bershady et al. 2010) shown in Cesare et al. (2020). In RG, dark matter is mimicked by the gravitational permittivity, a monotonic function of the local mass density that is expected to be universal. For the permittivity, according to Matsakos & Diaferio (2016) and Cesare et al. (2020), we adopted a smooth function that depends on three parameters: the permittivity of the vacuum ϵ₀, the critical density ρ_c, and the transition slope Q.

We find that the sets of parameters of the three E0 galaxies are consistent with each other within 2σ. Their values averaged over the three galaxies are {ϵ₀, Q, 𝒫_c} = $\{0.{089}_{-0.035}^{+0.038},0.{47}_{-0.21}^{+0.29},-24.{25}_{-0.20}^{+0.28}\}$, with 𝒫_c = log₁₀(ρ_c [g cm⁻³]). With this permittivity, RG requires stellar mass-to-light ratios in agreement with SPS models and tangential or radial orbits for the GCs, depending on the galaxy and the colour of the GCs. In particular, for NGC 1407, the mass-to-light ratio and the orbital anisotropy parameters are within 1σ from the values found by Pota et al. (2015), who modelled this galaxy with Newtonian gravity and a dark matter halo.

By modelling the kinematic properties of each individual DMS disk galaxy, Cesare et al. (2020) derived the mean values of the permittivity parameters {ϵ₀, Q, 𝒫_c}={0.56 ± 0.16, 0.92 ± 0.71, −25.30 ± 1.22}. The mean values of Q and 𝒫_c are within 1σ from the parameters we find here, whereas ϵ₀ is within 3σ. These combined results support the universality of the parameters of the RG permittivity.

Cesare et al. (2020) also derived a set of global parameters $\{{ϵ}_0,Q,{\mathcal{P}}_{\mathrm{c}}\}=\{0.{661}_{-0.007}^{+0.007},1.{79}_{-0.26}^{+0.14},-24.{54}_{-0.07}^{+0.08}\}$ with a MCMC analysis where the mass-to-light ratio and the disk thickness of each individual galaxy were kept fixed. This simplified approach returns an approximate estimate of the unique set of the permittivity parameters capable of describing the entire sample of the 30 DMS disk galaxies. This analysis of Cesare et al. (2020) is suggestive of the possible universality of the permittivity parameters, but their actual values should be taken with caution, because of the approximate procedure used to derive them. Nevertheless, the global DMS values of 𝒫_c and Q found by Cesare et al. (2020) are within 1.4σ and 3.4σ, respectively, from the values we find here and thus are still consistent with the expected universality. This consistency is particularly important for the critical density 𝒫_c, which sets the density scale of the transition between the Newtonian and the RG regimes. On the contrary, the global DMS ϵ₀ is 14.8σ discrepant from the mean value that we find here. This discrepancy is mostly driven by the width of the posterior distribution found by Cesare et al. (2020) which is small when compared to the uncertainties of the corresponding mean values (Fig. 9). If this width is not severely underestimated, this discrepancy might highlight at least three different possible problems: a too simplistic model of the ellipticals, an incorrect form of the permittivity, or a fundamental fault of RG.

The model of the ellipticals that we adopt here is oversimplified, because we treat the galaxies as isolated systems and completely neglect the net rotation of the galaxies. Indeed, NGC 1407 and NGC 5846 are within galaxy groups and NGC 4486 is the central galaxy of the Virgo cluster. Neglecting the galaxy environment can clearly overlook relevant interactions with neighbouring galaxies that could actually affect the parameters of the permittivity. Including these environmental effects requires modelling all the galaxies within the system at the same time, with each galaxy modelled according to its morphological type and shape. Properly modelling such a system is arduous and its non-linear effects are difficult to assess at the present stage. Thus, we cannot exclude that environmental effects might be part, if not all, of the cause of the possible tension on ϵ₀ that we find here. Alternatively, the explicit form of the permittivity ϵ(ρ) (Eq. (4)), which is chosen arbitrarily, might actually be inappropriate, especially in the very low-density regime: hints from the ultra-weak field limit of the covariant version of RG (Sanna et al. 2021) might suggest a better motivated permittivity that might alleviate this tension. Finally, all of the above might turn out to play a negligible role, and this possible tension on the permittivity of the vacuum might simply betray a fundamental fault of RG.

Acknowledgments

We sincerely thank the referee who identified a few errors in our original analysis and suggested relevant improvements. We also thank Xavier Hernández, Federico Lelli, Filippo Fraternali, and Heng Yu for their precious advice and detailed comments that improved the presentation of this work. We acknowledge the help of Alessandro Paggi and Valentina Missaglia in modelling the gas of NGC 5846 and the number density profiles of the GCs. We also thank Michele Cappellari, for his help in the derivation of the stellar root-mean-square velocity dispersion profiles of NGC 4486 and NGC 5846. We thank Compagnia di San Paolo (CSP) for funding the graduate-student fellowship of VC. We also acknowledge partial support from the INFN grant InDark and the Italian Ministry of Education, University and Research (MIUR) under the Departments of Excellence grant L.232/2016. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: Parameters of the stellar surface brightness profiles of NGC 4486 and NGC 5846

Tables A.1 and A.2 list the parameters of the model of the stellar surface brightness that we adopted for the galaxies NGC 4486 and NGC 5846. These parameters were determined by Scott et al. (2013) by fitting the measured surface brightness with the MGE model (Sect. 3.1.2).

Appendix B: Uncertainty on the root-mean-square velocity dispersions of GCs

We adopted the procedure of Danese et al. (1980) to estimate the uncertainties ΔV_rms(R) on the values of the root-mean-square velocity dispersion profile V_rms(R) that we computed with Eq. (22), where we considered the GCs in different radial bins. In the reference frame of the centre of mass of the galaxy, the line-of-sight velocity of the i-th GC is

$$\begin{array}{c}\hfill v_{||,i}=\frac{V_{\mathrm{rad},i}-V_{\mathrm{sys}}}{1+V_{\mathrm{sys}}/c},\end{array}$$(B.1)

where V_rad, i is the observed radial velocity of the GC, V_sys is the systemic velocity of the galaxy, and c = 3 × 10⁵ km s⁻¹ is the speed of light. The three E0 galaxies of our analysis have redshift z < 0.007 and the denominator of Eq. (B.1) can be approximated to 1. The variance of the line-of-sight velocities of the GCs in each radial bin of radius R is

$$\begin{array}{c}\hfill {\sigma }_{||}^2(R)=\frac{1}{N-1}{\displaystyle \sum _{i=1}^N{v}_{||,i}^2,}\end{array}$$(B.2)

where N is the number of GCs in the bin. This equation neglects the errors on V_rad and V_sys.

The quantity

$$\begin{array}{c}\hfill S^2(R)=(N-1)\frac{{\sigma }_{||}^2(R)}{V_{\mathrm{rms}}^2(R)}\end{array}$$(B.3)

is a ${\chi }_{\nu }^2$ random variate with ν = N − 1 degrees of freedom. We define $S_{-}^2$ and $S_{+}^2$ with the probability α of having S² in the range $(S_{-}^2,S_{+}^2)$. The value α = 0.68 defines the upper, indicated with the + sign, and lower, indicated with the − sign, 1σ uncertainties on the dispersion of the line-of-sight velocities derived with Eq. (B.2). The values $S_{-}^2$ and $S_{+}^2$ provide the following estimate of the uncertainty of V_rms(R) appearing in Eq. (22)

$$\begin{array}{c}\hfill \mathrm{\Delta }{V}_{\mathrm{rms},\pm }^2(R)={[{\left(\frac{\nu }{S_{\mp }^2(R)}\right)}^{1/2}-1]}^2{\sigma }_{||}^2(R)+\frac{\overline{{\delta }_{\ast }^2}}{N}(1+\frac{\overline{{\delta }_{\ast }^2}}{2{\sigma }_{||}^2(R)}),\end{array}$$(B.4)

where $\overline{{\delta }_{\ast }^2}$ is

$$\begin{array}{c}\hfill \overline{{\delta }_{\ast }^2}=\frac{\mathrm{\Delta }{V}_{\mathrm{sys}}^2}{{(1+V_{\mathrm{sys}}/c)}^2}\end{array}$$(B.5)

and ΔV_sys is the uncertainty on V_sys. In our analysis, V_sys ≪ c, and thus we can replace $\overline{{\delta }_{\ast }^2}$ with $\mathrm{\Delta }{V}_{\text{sys}}^2$. In Sect. 4.2.1, we used the symmetrised uncertainty

$$\begin{array}{c}\hfill \mathrm{\Delta }{V}_{\mathrm{rms}}(R)=\frac{\mathrm{\Delta }{V}_{\mathrm{rms},-}(R)+\mathrm{\Delta }{V}_{\mathrm{rms},+}(R)}{2}.\end{array}$$(B.6)

All Tables

All Figures

	Fig. 1. Surface number density profiles of blue (blue symbols and blue lines) and red (red symbols and red lines) GCs. Blue and red dots with error bars show the binned profiles; the solid lines are the models.
In the text

	Fig. 2. Three-dimensional mass density profile of the hot gas of NGC 5846. The solid circles with error bars show the discrete profile derived by Paggi et al. (2017); the curve shows our model.
In the text

Fig. 3. Root-mean-square velocity dispersion profiles of the stars (upper panels) and of the GCs (lower panels). Each column refers to the galaxy indicated in the upper panels. Upper middle and right panels: grey dots and black dots with error bars show the unbinned and binned profiles, respectively. Lower panels: red and blue dots with error bars refer to the blue and red GC populations, respectively. We note the different radial and velocity ranges of the panels.

In the text

Fig. 4. Cumulative mass (upper panels) and mass density (lower panels) profiles of the baryonic components of each galaxy: stars (black solid line), blue GCs (blue solid line), red GCs (red solid line), and hot X-ray emitting gas (magenta solid line); the green solid lines in the upper panels show the SMBH contributions. The purple solid lines show the sum of all the baryonic contributions. The stellar profiles assume Υ = 7.6 M⊙/L⊙ for NGC 1407, and Υ = 3.6 M⊙/L⊙ for NGC 4486 and NGC 5846. The grey shaded area around the stellar profiles shows the mass variation by adopting mass-to-light ratios in the range [4.0, 11.2] M⊙/L⊙, for NGC 1407, and [1.7, 5.5] M⊙/L⊙, for NGC 4486 and NGC 5846, in the B- and r-band, respectively. The profiles of blue and red GCs assume a GC mass MGC = 105 M⊙; the blue and red shaded areas show the mass variation by adopting MGC in the range [104, 106] M⊙. The magenta shaded areas show the uncertainties on the adopted gas mass density profiles; the gas profiles of NGC 4486 do not have a shaded area because the uncertainties on the profiles are unavailable. The green solid lines and shaded areas in the upper panels are the masses of the SMBH and their uncertainties. The purple shaded areas show the possible range of the total baryonic profiles. Each column refers to the galaxy indicated in the upper panels.

In the text

Fig. 5. Posterior distributions of the seven parameters derived from the two joint Markov chains for NGC 1407. The parameters are the stellar mass-to-light ratio, Υ, the velocity anisotropy parameters of the three tracers, ℬ*, ℬB, ℬR, and the three permittivity parameters, ϵ0, Q, and 𝒫c. The red squares with error bars show the medians and their uncertainites, set by the 15.9 and 84.1 percentiles of the posterior distributions. The yellow, cyan, and blue contours limit the 1σ, 2σ, and 3σ regions, respectively. The medians and their uncertainites are also reported above each column and in the one-dimensional posterior distributions in the top panels of each column as red solid and black dashed lines, respectively.

In the text

	Fig. 6. Same as Fig. 5, but for NGC 4486.
In the text

	Fig. 7. Same as Fig. 5, but for NGC 5846.
In the text

Fig. 8. Models of the root-mean-square velocity dispersion profiles emerging from our MCMC analysis in RG (solid lines). Upper panels: the green, blue, and red solid lines show the model profiles for the stars, blue GCs, and red GCs, respectively. Lower panels: zoom-in of the stellar profiles. The solid circles with error bars are the measured profiles from Fig. 3. The green, blue, and red dashed lines show the kinematic profiles of the stars, blue GCs, and red GCs computed with Eq. (23) where the values of the mass-to-light ratio Υ and the velocity anisotropy parameters β are those of the RG model but the gravitational field is Newtonian with no dark matter. The discrepancy with the RG curves highlights the gravitational boosting generated by the RG gravitational permittivity.

In the text

Fig. 9. Permittivity parameters estimated from the three E0 galaxies in our sample (orange squares with error bars) compared with the permittivity parameters estimated by Cesare et al. (2020) from the DMS disk galaxies. The purple squares with error bars show the means of the permittivity parameters found for the individual DMS galaxies. The light blue shaded areas show the posterior distributions of the three permittivity parameters found with an approximate procedure from the entire DMS sample, with the green dots indicating their median values and the yellow, red, and black contours indicating the 1σ, 2σ, and 3σ levels, respectively.

In the text