Origin of the fundamental plane of elliptical galaxies in the Coma Cluster without fine-tuning

TL;DR

This paper explains the tilt of the fundamental plane of elliptical galaxies in the Coma Cluster using modified Newtonian dynamics (MOND), without requiring fine-tuning, by analyzing galaxy data with specific mass models and assumptions.

Contribution

It demonstrates that the fundamental plane tilt can be naturally explained by MOND with a simple interpolating function and realistic assumptions, avoiding fine-tuning.

Findings

The tilt is explained by MOND with a simple interpolating function.

A well-motivated anisotropic velocity distribution is assumed.

The analysis supports the use of Kroupa or Salpeter stellar mass-to-light ratios.

Abstract

After thirty years of the discovery of the fundamental plane, explanations to the tilt of the fundamental plane with respect to the virial plane still suffer from the need of fine-tuning. In this paper, we try to explore the origin of this tilt from the perspective of modified Newtonian dynamics (MOND) by applying the 16 Coma galaxies available in Thomas et al.[1]. Based on the mass models that can reproduce de Vaucouleurs' law closely, we find that the tilt of the traditional fundamental plane is naturally explained by the simple form of the MONDian interpolating function, if we assume a well motivated choice of anisotropic velocity distribution, and adopt the Kroupa or Salpeter stellar mass-to-light ratio. Our analysis does not necessarily rule out a varying stellar mass-to-light ratio.

Equations

$${\sigma }_{\mathrm{ap}}^2=\frac{G\mathcal{M}}{R_h}\mathcal{H}(\tilde{R}){\Gamma }_h^M(\tilde{R},\beta ,{\rm Y}{I}_e)$$

$$\displaystyle\Gamma^{M}_{h}=2\int^{\tilde{\it{R}}}_{0}\tilde{\it{R}}^{\prime}d\tilde{\it{R}}^{\prime}\int^{\infty}_{\tilde{\it{R}}^{\prime}}\frac{\tilde{\it{r}}d\tilde{\it{r}}}{\sqrt{\tilde{\it{r}}^{2}-\tilde{\it{R}}^{\prime 2}}}\Big{(}1-\beta(\tilde{\it{r}})\frac{\tilde{{\it{R}}}^{\prime 2}}{\tilde{\it{r}}^{2}}\Big{)}$$

$$\displaystyle\int^{\infty}_{\tilde{\it{r}}}\frac{\tilde{\nu}(\Upsilon{\,}I_{e},\tilde{\it{r}}^{\prime})}{\tilde{\it{r}}^{\prime}(\tilde{\it{r}}^{\prime}+1)^{5}}\exp\Big{(}\int^{\tilde{\it{r}}^{\prime}}_{\tilde{\it{r}}}\frac{2\beta}{{\tilde{\it{r}}^{\prime\prime}}}{d\tilde{\it{r}}^{\prime\prime}}\Big{)}{d\tilde{\it{r}}^{\prime}}$$

$${\sigma }_{\mathrm{ap}}^2=2\pi G\times \eta (\tilde{R},\beta )\times {\rm Y}\times R_e\times I_e\times {\Xi }^M(\tilde{R},\beta ,{\rm Y}{I}_e),$$

$$\eta (\tilde{R},\beta )=R_h/R_e\times \mathcal{H}(\tilde{R}){\Gamma }_h^N(\tilde{R},\beta )$$

$${\Xi }^M={\Gamma }_h^M/{\Gamma }_h^N$$

$$\displaystyle\Gamma^{M}_{j}=2\int^{\tilde{\it{R}_{j}}}_{0}\tilde{\it{R}_{j}}^{\prime}d\tilde{\it{R}_{j}}^{\prime}\int^{\infty}_{\tilde{\it{R}_{j}}^{\prime}}\frac{\tilde{\it{r}_{j}}d\tilde{\it{r}_{j}}}{\sqrt{\tilde{\it{r}_{j}}^{2}-\tilde{\it{R}_{j}}^{\prime 2}}}\Big{(}1-\beta(\tilde{\it{r}_{j}})\frac{\tilde{{\it{R}_{j}}}^{\prime 2}}{\tilde{\it{r}_{j}}^{2}}\Big{)}$$

$$\displaystyle\int^{\infty}_{\tilde{\it{r}_{j}}}\frac{\tilde{\nu}(\Upsilon{\,}I_{e},\tilde{\it{r}_{j}}^{\prime})}{\tilde{\it{r}_{j}}^{\prime 3}(\tilde{\it{r}_{j}}^{\prime}+1)^{3}}\exp\Big{(}\int^{\tilde{\it{r}_{j}}^{\prime}}_{\tilde{\it{r}_{j}}}\frac{2\beta}{{\tilde{\it{r}_{j}}^{\prime\prime}}}{d\tilde{\it{r}_{j}}^{\prime\prime}}\Big{)}{d\tilde{\it{r}_{j}}^{\prime}},$$

$${\beta }_e(r)=\frac{r^2}{(r^2+R_e^2)}.$$

$$\log \frac{R_e}{\mathrm{kpc}}=a\log \frac{{\sigma }_{\mathrm{ap}}}{\mathrm{km}{\mathrm{s}}^{-1}}+b\log \left(\frac{I_e}{{\mathrm{L}}_{\odot }{\mathrm{pc}}^{-2}}{\Xi }^M\right)+c,$$

$$c=-\log \frac{⟨{\rm Y}⟩}{{\mathrm{M}}_{\odot }/{\mathrm{L}}_{\odot }}-\log \left(2\pi \frac{G}{(\mathrm{km}{\mathrm{s}}^{-1}{)}^2\mathrm{kpc}{\mathrm{M}}_{\odot }^{-1}}\eta \right).$$

$$\log \frac{R_e}{\mathrm{kpc}}=a\log \frac{{\sigma }_{\mathrm{ap}}}{\mathrm{km}{\mathrm{s}}^{-1}}+b\log \left(\frac{{\Sigma }_e}{{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}}{\Xi }^M\right)+c,$$

$$c=-\log \left(2\pi \frac{G}{(\mathrm{km}{\mathrm{s}}^{-1}{)}^2\mathrm{kpc}{\mathrm{M}}_{\odot }^{-1}}\eta \right).$$

$$\tilde{\nu }(x_N)=\frac{1}{2}\left(1-\frac{g_e}{g_N}+\sqrt{4(\frac{a_0}{g_N}+\frac{g_e}{g_N})+(1-\frac{g_e}{g_N}{)}^2}\right),$$

Full text

I Introduction

The first scaling law of elliptical galaxies, the Faber-Jackson relation, was found by (Faber and Jackson, 1976) in 1976, and it was later incorporated into a fundamental plane (FP) between the size, surface brightness, and velocity dispersion of elliptical galaxies (Djorgovski and Davis, 1987; Dressler et al., 1987; Jørgensen et al., 1996). However, this empirical FP does not comply with the dimensional-analysis of the virial theorem in the standard Newtonian dynamics while a constant total mass-to-light ratio ($M/L$) is assumed (e.g. Gerhard et al., 2001; Treu et al., 2005). This deviation of the observed FP to the theoretical prediction is called the tilt of the FP. The tilt is usually attributed to a varying total $M/L$ that increase with mass. In the past decade, thanks to the independent, accurate measurement of total mass by stellar dynamics (Cappellari et al., 2006; Thomas et al., 2011; Cappellari et al., 2013) or strong lensing (Bolton et al., 2007, 2008; Auger et al., 2010), a genuine variation of total $M/L$ has been confirmed. Not only the total $M/L$ derived from these two different techniques are consistent, but also both techniques show that other than the FP, there is a less scattering total mass plane (MP) that agrees well with the prediction of the virial theorem (Bolton et al., 2008; Thomas et al., 2011; Cappellari et al., 2013). Based on this independent estimation of total mass, it has also been shown that the total $M/L$ is highly correlated to the residuals of the FP, so a variation of $M/L$ among galaxies is responsible for the tilt of the FP (Bolton et al., 2007).

Despite of the clear evidences of the systematic variation of total $M/L$ among galaxies, the origin of this variation, so the tilt of the FP, is an ongoing debate (Desmond and Wechsler, 2016; D’Onofrio et al., 2016). Explanations include structural non-homology in elliptical galaxies (Trujillo et al., 2004), variations in the distribution of dark matter (Renzini and Ciotti, 1993), and the systematic variation of stellar population (Forbes et al., 1998). Of these, non-homology not only requires fine-tuning between stellar $M/L$ and Sersic index, but also just plays a minor role in the tilt of the FP at best (Bolton et al., 2008; Cappellari et al., 2013). On the other hand, the total $M/L$ used to derive the MP in Thomas et al. (2011) obviously differs from the stellar $M/L$ estimated from the population synthesis, which implies that the tilt of the FP is not a pure stellar population effect. Instead, either a non-universal stellar initial mass function (IMF) or variations of dark matter fraction among systems shall account for, at least partly, the tilt of the FP. However, the explanations from IMF and dark matter both require strong fine-tuning (Ciotti et al., 1996).

Since the nature of dark matter can contribute to the tilt of the FP, it is only one step forward to ask whether or not the tilt of the FP can also be attributed to other alternatives to dark matter. Modified Newtonian dynamics (MOND), as an alternative paradigm to dark matter (Milgrom, 1983a), is very successful at galactic systems, not only the overwhelmingly convincing baryonic Tully-Fisher relation (McGaugh, 2005, 2011, 2012) and rotation-curve analysis in spiral galaxies (Sanders and McGaugh, 2002; Famaey and McGaugh, 2012), but also analysis of galactic strong lensing (Chiu et al., 2006, 2011; Zhao et al., 2006; Tian et al., 2013), galaxy-galaxy lensing (Milgrom, 2013), dynamical analysis of planetary nebulae (Tian and Ko, 2016), the velocity dispersions of the dwarf satellites of Andromeda (McGaugh and Milgrom, 2013a, b; Pawlowski and McGaugh, 2014), gas ring in elliptical galaxies (Milgrom, 2012), and newly found relation between the baryonic and dynamical central surface densities of disc galaxies (Lelli et al., 2016; Milgrom, 2016), and radial acceleration in rotating galaxies (McGaugh et al., 2016). However, although the Tully-Fisher relation of spiral galaxies is naturally explained in MOND (Milgrom, 1983b), the origin of the FP of elliptical galaxies in MOND is much less clear. Several papers have attempted to explain the FP in MOND (Sanders, 2000; Sanders and Land, 2008; Cardone et al., 2011), but none of them have explicitly applied observational data to compare exponents of the observed MOND ‘fundamental plane’ with their theoretical counterpart. To the best of our knowledge, this paper is the first attempt to do so.

In this paper, we apply $16$ Coma elliptical galaxies in Thomas et al. (2011) to construct the FP and MP in MOND, where stellar mass is tantamount to total mass, because the stellar $M/L$ estimated from stellar population synthesis is available for every member of these galaxies. In addition, since some of these galaxies have more than $50\%$ of dark matter within one effective radius, $\rm{R}_{\rm{e}}$ (Thomas et al., 2011), we expect a substantial MONDian effect will be detected in these data.

The structure of the paper is organized as below. In Sec. II, we discuss the dynamics of elliptical galaxies in MOND, and show how to obtain the theoretical MOND FP in the following section. We then construct the observational MOND FPs from the $16$ Coma elliptical galaxies, and compare them with their theoretical counterparts in Sec. IV. Finally, we give a brief discussion and conclusion in Sec. V.

II Dynamics of elliptical galaxies in MOND

For spherical systems, the Jeans equation in the context of MOND only differs from its counterpart in the Newtonian dynamics in the radial gravitational force, $g=\tilde{\nu}(x_{N})g_{N}$ (Sanders, 2000; Cardone et al., 2011; Chae and Gong, 2015; Tian and Ko, 2016), where $x_{N}$ is the ratio between the Newtonian gravity, $g_{N}$, and the MOND acceleration constant, $a_{0}$ ($=1.2\times 10^{-8}\rm{cm\,s}^{-2}$), and where $\tilde{\nu}$ is the interpolating function of MOND (e.g. Milgrom, 1983a; Chiu et al., 2011; Ko, 2016). Following exactly the same way as in the Newtonian dynamics (e.g. Binney and Mamon, 1982; Agnello et al., 2014), we can solve for the aperture velocity dispersion, $\sigma_{\rm{ap}}$, from the modified Jeans equation in MOND with certain matter density and anisotropy profile of velocity distribution, $\beta$. For the matter density, we assume that mass always follows light, so the Hernquist model (Hernquist, 1990), which approximates de Vaucouleurs’ law closely, shall describe all mass in elliptical galaxies. Accordingly, we obtain

[TABLE]

where $G$ is Newton’s constant, $\mathcal{M}$ is total mass of galaxy, $\Upsilon$ is stellar mass-to-light ratio, $I_{e}$ is surface brightness within the effective radius, and $R_{h}$, a characteristic length in the Hernquist model, is approximately related to the effective radius by 0.551$$R_{e}. In Eq. (1), we also denote $\mathcal{H}(\tilde{R})$ as a ratio between total luminosity and surface brightness within $\tilde{R}$, which itself is a radius normalised to $R_{h}$, and use

[TABLE]

to describe the component that will be influenced by MOND. For any fixed aperture radius, $\Gamma^{M}_{h}$ will in general depend on anisotropy, surface brightness within the effective radius, mass-to-light ratio, and the MONDian interpolating function. On the contrary, when the MONDian effect vanishes, $\Gamma^{M}_{h}$ will only depend on $\beta$. We denote this special case of $\Gamma^{M}_{h}$ as $\Gamma^{N}_{h}(\tilde{R},\beta)$. According to the notations above, Eq. (1) can be rewritten as

[TABLE]

where

[TABLE]

is defined as a structural parameter that is solely decided by mass model and anisotropy, and

[TABLE]

is a MONDian parameter that describes how much $\sigma_{\rm{ap}}$ will deviate from its Newtonian counterpart. Basically, Eq. (3) plays a role similar to the virial theorem in the traditional analysis of the FP.

We would like to emphasize that although the derivation of Eq. (3) is based on the Hernquist model, we can easily extend it to other mass models, and will obtain similar results. For example, for the Jaffe model (Jaffe, 1983), which is another analytical mass profile that successfully reproduces de Vaucouleurs’ law, we can obtain an equation similar to Eq. (3) except for replacing $\Gamma^{M}_{h}$ in the structural parameter, $\eta$, and the MONDian parameter, $\Xi^{M}$, with

[TABLE]

where the radii are normalised to the characteristic length in the Jaffe model, $R_{j}=1.3106R_{e}$, and for substituting $R_{h}$ and $\mathcal{H}(\tilde{R})$ in the structural parameter with $R_{j}$ and $\mathcal{H}(\tilde{R}_{j})$.

III The MONDian Fundamental Plane

III.1 Anisotropy, Mass-to-light ratios, and MONDian effects

In this paper, we only consider three most well-known forms of the MONDian interpolating function, that is, the standard, simple, and Bekenstein forms (e.g. Famaey and McGaugh, 2012). As shown in Fig. 1, under the same conditions, the Bekenstein form will yield the strongest MONDian effect to boost the velocity of dispersion; on the contrary, the standard form will lead to the weakest. Figure 1 also shows that for a given surface brightness, the larger the stellar $M/L$, the weaker the MONDian effect; on the other hand, like spiral galaxies (McGaugh and de Blok, 1998), for a fixed stellar $M/L$, the lower surface brightness will have stronger MONDian effects. In addition, Fig. 1 demonstrates the values of $\Xi^{M}$ of the 16 Coma galaxies considered in this paper for two different types of $\Upsilon$: the Kroupa ($\Upsilon_{\rm{Krou}}$) and Salpeter ($\Upsilon_{\rm{Salp}}$) mass-to-light ratio in Thomas et al. (2011). The values of $I_{e}\Upsilon_{\rm{Krou}}$ and $I_{e}\Upsilon_{\rm{Salp}}$ of these 16 galaxies range from around $100$ to $6000$ $\rm{M}_{\odot}\rm{pc}^{-2}$. Finally, we compare the influences of anisotropy on $\Xi^{M}$. We consider anisotropy profile $\beta=r^{2}/(r^{2}+r_{a}^{2})$ (van Albada, 1982; Sanders, 2000; Milgrom and Sanders, 2003; Agnello et al., 2014), which fits well the velocity distributions of the systems formed by dissipationless collapse (van Albada, 1982; Sanders, 2000), but set $r_{a}=R_{e}$ to avoid an ad hoc free parameter, and obtain

[TABLE]

To check the sensitivity of $r_{a}$, we also consider $r_{a}=3R_{e}$, which is used by Milgrom and Sanders (2003) to fit the line-of-sight velocity dispersion and the planetary nebulae observations. This alternation will only cause less than $5\%$ changes of $\Xi^{M}$. Based on $\beta_{e}$, we find the stark differences of the influences of anisotropy profile between the Newtonian dynamics and MOND, which have never been noticed before. In the Newtonian dynamics, the contribution of $\beta_{e}$ to the boost of $\sigma_{\rm{ap}}$ is a constant, and only $10\%$ more than its counterpart in the case of isotropy. However, in MOND, the impact of introducing $\beta_{e}$ into Eq. (2) is not only much larger, but also not a constant: it depends on $\Upsilon$, $I_{e}$, and $\tilde{\nu}$.

III.2 The tilt of the Fundamental Plane

While we assume a universal stellar $M/L$ in elliptical galaxies, the existence of a theoretical MOND FP between $R_{e}$, ${\sigma_{\rm{ap}}}$, and $I_{e}{\Xi^{M}}$ can be easily recognized from Eq. (3), which leads to

[TABLE]

with $a=2$, $b=-1$, and

[TABLE]

Here, we take the mean of $\Upsilon$, $\langle\Upsilon\rangle$, to guarantee that $c$ is a constant because even for a universal stellar $M/L$, there must exist some degree of random variations of $\Upsilon$ for individual galaxies. Equation (8) is identical to its Newtonian counterpart except for the additional parameter of the MONDian effects, $\Xi^{M}$, which implicitly includes the MOND acceleration constant, $a_{0}=1.2\times 10^{-8}\rm{cm\,s}^{-2}$. Unlike the Newtonian FP derived from the virial theorem, however, we do not use a free parameter (i.e. a structure coefficient) to ambiguously represent galactic structure (e.g. Bolton et al., 2007; Thomas et al., 2011; Cappellari et al., 2013); instead, the influence of galactic structure on the MOND FP is analytically decided by mass models of galaxies. Since we assume a Hernquist model throughout our derivation, according to Eq. (4), the structure parameter, $\eta$, will only depend on $\beta$ for any fixed radius, $\tilde{R}$. This fact also implies that homology of galaxies (galactic structure is independent of the galaxy mass) is implicitly assumed in Eq. (8).

Following Bolton et al. (2007), we might also incorporate variations of stellar $M/L$ among galaxies into the MOND FP by combining $\Upsilon$ with $I_{e}$ in Eq. (8), so that surface brightness will be replaced with surface matter density, and the MOND FP will become the MOND ‘mass plane’,

[TABLE]

where ${\Sigma}_{e}$ is the surface matter density within the effective radius, and the coefficients $a$ and $b$ remain the same as in Eq. (8), but

[TABLE]

This is the first time that the MOND FP and MOND MP are written down explicitly and lucidly. According to Eqs (8)-(11), the tilt of the FP is attributed to a single parameter, $\Xi^{M}$.

Although ideally the MOND FP and MOND MP are identical for a universal stellar $M/L$, in practice, these two planes are not exactly the same. Firstly, we need $\langle\Upsilon\rangle$ in Eq. (9), but individual $\Upsilon$ to calculate $\Sigma_{e}$ in Eq. (10). Secondly, even if there is no random variation around $\langle\Upsilon\rangle$, the MOND MP will still differ from the MOND FP in the exponent of $\Upsilon$, which together with $I_{e}$ is implicitly bound up with the exponent of $\Sigma_{e}$ in the MOND MP, but is always equal to $-1$ in the MOND FP. These subtle differences are non-trivial when we try to find and compare best-fitting coefficients of the MOND FP and the MOND MP.

IV Data analysis

Similar to Bolton et al. (2007) and Thomas et al. (2011), we perform regression analysis of the 16 Coma galaxies in order to find the best-fitting coefficients of the observed MOND FP and MOND MP for various $\Upsilon$, $\beta$, and $\tilde{\nu}$ outlined above. The error bars of these coefficients are estimated from one standard deviation of the mean of bootstrap resampling. Unlike the traditional regression analysis of the FP (e.g. Bolton et al., 2007; Thomas et al., 2011; Cappellari et al., 2013), however, the intercepts of the observed MOND FP and MOND MP in our analysis are not decided by the best-fitting value, but calculated from Eq. (9) and Eq. (11) accordingly, where $\Upsilon$ in Eq. (9) is taken from the mean value of $\Upsilon_{\rm{Salp}}$ or $\Upsilon_{\rm{Krou}}$ of the 16 Coma galaxies in Thomas et al. (2011). By adopting the no-intercept regression analysis, we will have fewer free parameters to fit.

Firstly, we study the FP and MP of the 16 galaxies without any MONDian effect (left end of Fig. 2). Since the contributions of $\beta_{e}$ to the these planes are independent of $\sigma_{\rm{ap}}$, $R_{e}$, and $I_{e}$, we do not expect any obvious deviation of our results from previous studies. The consistency of our fit with that of (Bernardi et al., 2003; Thomas et al., 2011) within the statistical uncertainties reinforces this expectation.

Among various assumptions of the dynamics of elliptical galaxies in MOND, we find that the assumption of isotropic velocity dispersion of galaxies ($\beta=0$) will lead to great difficulties in removing the tilts of the empirical MOND FP and MOND MP with respect to their theoretical counterparts (the left end of Fig. 2). Indeed, we can hardly improve the tilts in MOND for $\beta=0$ whatever the interpolating functions we choose. In contrast, while taking the anisotropy profile $\beta_{e}$, we find that the observed MOND FPs and MOND MPs both agree well with their theoretical counterparts along with a universal stellar $M/L$ and simple form of $\tilde{\nu}$. Under these conditions, the best-fitting coefficients of the MOND FP for $\Upsilon_{\rm{Krou}}$ are $a=1.92\pm 0.15$ and $b=-0.88\pm 0.12$ with ${rms}=0.16$, and that of the MOND MP are $a=2.10\pm 0.16$ and $b=-1.03\pm 0.11$ with ${rms}=0.15$ (Figure 2). Similarly, if we replace $\Upsilon_{\rm{Krou}}$ by $\Upsilon_{\rm{Salp}}$ but leave the other conditions unchanged, we will obtain the best-fitting coefficients, $a=1.96\pm 0.13$ and $b=-0.87\pm 0.11$, of the MOND FP as well as $a=2.04\pm 0.16$ and $b=-0.96\pm 0.10$ for the MOND MP; both planes have ${rms}=0.15$. It is worth noting that we can hardly detect any change of the slopes if we replace ${R}_{e}$ with $3R_{e}$ in $\beta_{e}$. In addition, the larger rms errors of the MOND FPs and MOND MPs of $\beta_{e}$, compared to that of $\beta=0$, might be largely attributed to the fact that although $\beta_{e}$ contributes a lot to the slopes of MOND FPs and MOND MPs, it is at best an averaged approximation to an ensemble of real anisotropies of the 16 coma galaxies. Unlike the simple form, although the standard form with $\beta_{e}$ yields a very good value of the coefficient $b$, its another best-fitting coefficient $a$ is not so good; similarly, in the case of the Bekenstein form, the best-fitting coefficient $b$ is not as good as $a$.

In clusters of galaxies, the MONDian effect in elliptical galaxies might be affected by the external fields of the host clusters. To investigate the influence of this external field effect on our sample, we follow Bekenstein and Milgrom (1984); Famaey et al. (2007a), and express this effect in terms of the modification $\tilde{\nu}-$function. For the simple form, the $\tilde{\nu}-$function becomes

[TABLE]

where $g_{e}$ denotes the external field. In Fig. 3, we show how the external field effect will affect the slope of the MOND MP by considering $g_{e}=0.1a_{0}$, $0.5a_{0}$, and $a_{0}$. Although overall the external effect has only small impact on the best-fitting coefficients of the MOND MPs in our sample, we do find that a stronger external field will move the plane a little bit closer to the Newtonian one. In reality, however, not all galaxies in the sample will be affected by an external field as strong as $g_{e}=a_{0}$. To consider a more accurate case, we thus roughly model the external field based on the Fig. 3 of richtler_etal11, and classify the 16 galaxies into three groups. We assign $g_{e}=a_{0}$ to the two galaxies (GMP $3414$ and $3510$) that are located at the inner part of the cluster, assign $g_{e}=0.5a_{0}$ to GMP $3792$, $2417$, and $2440$, whose projected distances from the center of the cluster are between $0.56$ and $0.67$ Mpc, and assign $g_{e}=0.1a_{0}$ to the other 11 galaxies which are at least 1.45 Mpc away from the center. The MOND MP of this more realistic modeling of the external field effect is consistent with the theoretical prediction of Eq. (10).

Finally, to complete our analysis, we also consider the Jaffe model as an alternative to the mass profile in our analysis (lower panel of Fig. 3). Although the influence of the external field effect on the slope of the MOND MP is a little bit stronger in the Jaffe model than that in the Hernquist model, the MOND MP based on our modeling of the external field still agrees well with the theoretical prediction. In general, we do not find any obvious disparity between the results of the Hernquist and the Jaffe models.

V Discussion and conclusion

This paper, for the first time since the existence of the FP was confirmed in 1987, successfully explain the origin of the FP without any fine-tuning. We explicitly show that the tilt of the fundamental plane can be naturally explained in the framework of MOND by a combined effect of anisotropy of velocity dispersion and modification of gravity. All the parameters involved in our model are either well-motivated, such as the anisotropy profile $\beta_{e}$ (van Albada, 1982; Sanders, 2000), and a universal stellar $M/L$ (Kroupa, 2001), or consistent with other independent observations (Famaey and Binney, 2006; Zhao and Famaey, 2006; Famaey et al., 2007b), such as the simple form of the interpolating function. This is striking, because although the MONDian interpolating function might arguably be contrived to replace the requirement of dark matter in spiral galaxies, it is never the case for the slope of the FP. Everything just fits together naturally.

We shall emphasize that the success of $\Upsilon_{\rm{Krou}}$ or $\Upsilon_{\rm{Salp}}$ at removing the tilt of the observed FP does not mean that our results rebut a non-universal IMF in MOND. Indeed, under the same conditions, $\Upsilon_{*,\rm{dyn}}$, a non-universal $M/L$ in Thomas et al. (2011), will also match the empirical MOND MP (but not the MFP) with its theoretical counterpart. Although $\Upsilon_{*,\rm{dyn}}$ contains the contribution of dark matter and only provides an upper limit to stellar mass, this result does support that a slight variant of stellar $M/L$ among galaxies is not ruled out. Hence, our finding does not necessarily contradict with Sanders (2000), but disagrees with Tortora et al. (2014) at that MOND must require variant IMF to fit observations.

Over the past years, much attention has been paid to elliptical galaxies in MOND (Cardone et al., 2011; Chiu et al., 2011; Ferreras et al., 2012; Milgrom, 2012; Chae and Gong, 2015; Tian and Ko, 2016), partly because it is important to learn whether or not MOND is as successful in elliptical galaxies as in spiral ones. Unlike cluster of galaxies, where massive neutrino might accumulate and contribute to invisible matter (Sanders, 2007; Milgrom, 2008), any requirement for dark matter in galaxies will be devastating to the MONDian paradigm. This paper offers an independent piece of evidence to support the success of MOND, either as a fundamental theory or some effective theory, at galactic scales.

Acknowledgements

We would like to thank Mordehai Milgrom and David Koo for helpful comments. MCC would also like to thank CMK for the hospitality during his stay in NCU. CGS and MCC are partly supported by the Chinese National Nature Science Foundation No. 11433003, the Shanghai Science Foundation No. 13JC14044400 and National Key Project 973 No. 2014CB845704. CMK is supported in part by the Taiwan Ministry of Science and Technology grants MOST 104-2923-M-008-001-MY3 and MOST 105-2112-M-008-011-MY3.