Galaxy haloes and Yukawa-like gravitational potentials

# Systematic biases on galaxy haloes parameters from Yukawa-like gravitational potentials

V.F. Cardone and S. Capozziello
Dipartimento di Scienze e Tecnologie dell’ Ambiente e del Territorio, Università degli Studi del Molise,
Contrada Fonte Lappone, 86090 - Pesche (IS), Italy
Dipartimento di Scienze Fisiche, Università degli Studi di Napoli Federico II”, Complesso Universitario
di Monte Sant’Angelo, Edificio N, via Cinthia, 80126 - Napoli, Italy
I.N.F.N. - Sezione di Napoli, Complesso Universitario di Monte Sant’Angelo, Edificio G, via Cinthia, 80126 - Napoli, Italy
Accepted xxx, Received yyy, in original form zzz
###### Abstract

A viable alternative to the dark energy as a solution of the cosmic speed up problem is represented by Extended Theories of Gravity. Should this be indeed the case, there will be an impact not only on cosmological scales, but also at any scale, from the Solar System to extragalactic ones. In particular, the gravitational potential can be different from the Newtonian one commonly adopted when computing the circular velocity fitted to spiral galaxies rotation curves. Phenomenologically modelling the modified point mass potential as the sum of a Newtonian and a Yukawa - like correction, we simulate observed rotation curves for a spiral galaxy described as the sum of an exponential disc and a Navarro - Frenk - White (NFW) dark matter halo. We then fit these curves assuming parameterized halo models (either with an inner cusp or a core) and using the Newtonian potential to estimate the theoretical rotation curve. Such a study allows us to investigate the bias on the disc and halo model parameters induced by the systematic error induced by forcing the gravity theory to be Newtonian when it is not. As a general result, we find that both the halo scale length and virial mass are significantly overestimated, while the dark matter mass fraction within the disc optical radius is typically underestimated. Moreover, should the Yukawa scale length be smaller than the disc half mass radius, then the logarithmic slope of the halo density profile would turn out to be shallower than the NFW one. Finally, cored models are able to fit quite well the simulated rotation curves, provided the disc mass is biased high in agreement with the results in literature, favoring cored haloes and maximal discs. Such results make us argue that the cusp/core controversy could actually be the outcome of an incorrect assumption about which theory of gravity must actually be used in computing the theoretical circular velocity.

###### keywords:
dark matter – gravitation – galaxies : kinematic and dynamics

## 1 Introduction

The picture of a spatially flat universe undergoing accelerated expansion is the nowadays accepted view of our cosmo. According to the successful concordance CDM model (Carroll et al., 1992; Sahni & Starobinski, 2000), there are two main ingredients in this scenario, namely dark matter (accounting for the clustering of the structures we observe) and the cosmological constant (dominating the energy budget and driving the cosmic speed up). The anisotropy spectrum of cosmic microwave background (de Bernardis et al. 2000; Komatsu et al. 2010), the galaxy power spectrum with the imprinted baryon acoustic oscillations (Eisenstein et al. 2005; Percival et al. 2010) and the Hubble diagram of Type Ia Supernovae (Kowalski et al. 2008; Hicken et al. 2009) represent an incomplete list of the wide amount of data this model is able to excellently reproduce in a single scenario.

From a theoretical point of view, however, the CDM is far to be satisfactory. First, the term is orders of magnitude larger than what expected from quantum field theory, while the ratio between its energy density and the matter one is coincidentally of order unity just today while it should have been much smaller or much larger than 1 over the rest of the universe history. Motivated by these unpleasing shortcomings, a plethora of models giving rise to a varying  - like term has been proposed mainly based on a scalar field evolving under the influence of its own self interaction potential. Needless to say, the absence of any candidate for this scalar field and the full arbitrariness in the choice of the potential are serious drawbacks of these models which therefore represent only a way to change the problems without actually solving them.

On galactic scales, the term gives a negligible contribution to the gravitational potential so that the classical Newtonian theory is usually adopted. While the evolution of the universe is driven by the cosmological constant, the formation of structure is mainly determined by the dark matter (DM) which provide the potential wells where baryons collapse to originate the visible component of galaxies. Numerical simulations allow to follow this process predicting the structure of DM haloes. Surprisingly, the theoretical expectations are not in agreement with observations on galactic scales. In particular, the density profile of DM haloes is expected to follow a double power - law with with , a characteristic length scale and giving the logarithmic slope in the inner regions. Although the debate on what the precise value of is remains open, what is granted is that should definitely be positive with for the most popular NFW model (Navarro et al., 1997). The DM haloes are thus described by cusped models or, in other words, the logarithmic slope never vanishes. On the contrary, rotation curves of low surface brightness galaxies (LSB) are definitely better fitted by cored models, i.e. in the inner regions, such as the pseudo - isothermal sphere, (Binney & Tremaine, 1987) or the Burkert model, (Burkert, 1995; Salucci & Burkert, 2000). As well reviewed in de Blok et al. (2010), cored models turn out to be statistically preferred over cusped ones from the fit to large samples of LSB rotation curves. Moreover, for those galaxies where a statistically acceptable fit for a cusp model is obtained, it turns out that the concentration parameter (defined later) is significantly larger than what expected from numerical simulations (see de Blok 2010 and refs. therein for further details).

From the above picture, it is clear that, although observationally successful on cosmological scales, the CDM model is far to be free of problems, especially if one also remember that there is up to now no laboratory final evidence for any of the many particles candidate to the role of DM. It is therefore worth going beyond the usual view and look at a radical revision of the underlying scheme. To this end, one can consider cosmic speed up as the first evidence of a breakdown of General Relativity (GR) as we know it. Rather than being due to a new actor on the scene, the accelerated expansion can be a consequence of gravity working in a different way than GR predicts. This consideration has attracted most interest towards braneworld - like models such as the five dimensions DGP models (Dvali et al., 2000; Lue & Starkman, 2003) or fourth order theories of gravity, where the GR Einstein - Hilbert Lagrangian is generalized by the introduction of a function of the scalar curvature (Capozziello, 2002; Nojiri & Odintsov, 2007; Capozziello & Francaviglia, 2008; Sotiriou & Faraoni, 2010; de Felice & Tsujikawa, 2010; Capozziello & Faraoni, 2010). It is worth stressing that, notwithstanding which is the correct modified gravity theory, should GR be incorrect on cosmological scales, one has to check whether the gravitational potential on galactic scale is. Most of modified theories indeed induce negligible changes to the gravitational potential on Solar System scale in order to pass the classical tests of gravity (Will, 1993), but this does not prevent to have significant deviations from the Newtonian potential on the much larger scale of galaxies where no direct experimental test is available.

Although modified gravity theories are investigated at cosmological scales as alternatives to dark energy, we are here interested in whether they can also impact the estimate of dark matter properties on galactic scales. Should indeed the gravitational potential be modified, the computation of the rotation curve must be done in this modified framework. On the contrary, one typically assumes that Newtonian mechanics holds and then constrains the halo model parameters by fitting the theoretical rotation curve to the observed one. We are here interested in investigating whether this incorrect procedure biases in a significant way the determination of the halo parameters and whether such a bias can explain the inconsistencies among theoretical expectations and observations. In a sense, we are wondering whether modified gravity could be a possible way to solve the cusp/core and similar problems of the DM scenario.

The plan of the paper is as follows. In Section 2, we present the modified gravitational potential used and detail the derivation of the rotation curve. It is worth noticing that Yukawa-like corrections emerges in any analytical -gravity model, except , where is the Ricci scalar adopted in the Hilbert-Einstein action. Section 3 describes how we estimate the bias induced on the halo model parameters by fitting data with the Newtonian potential while the underlying theory of gravity is non-Newtonian. The results of this analysis are discussed in Sections 4 and 5, while Section 6 summarizes our conclusions.

## 2 Yukawa-like gravitational potentials

As it is well known, Newtonian mechanics is the low energy limit of GR. Indeed, looking for a stationary and spherically symmetric solution of the Einstein equations gives the Schwarzschild metric whose component gives to the Newtonian gravitational potential in the weak field limit. Modifying GR leads to modified field equations hence to a different solution and potential in the low energy limit. As such, we should first choose a modified theory of gravity to finally get the modified potential.

An interesting example for its application to cosmology is provided by theories of gravity. Writing the metric in the weak field limit as

 ds2=−[1−2A(r)]dt2+[1+2B(r)]dr2+r2dΩ2 ,

and using the conformal equivalence with scalar - field theories, Faulkner et al. (2007) have shown that :

 A(r)=~A(~r)+ψ(~r)√6Mpl (1)

where tilted quantities are evaluated in the Einstein frame (with and ) and is the scalar field coupled to matter determined by :

 1~r2dd~r(~r2dψd~r)=∂Veff(ψ,~r)∂ψ . (2)

The effective potential is determined from both the functional expression for and the local mass density as :

 Veff=V(ψ)+χ−1/2¯ρ(~r) (3)

with given by Eq.(6) in Faulkner et al. (2007) and . For the particular case of a uniform sphere of mass and radius and a quadratic potential, one finally gets :

 ϕ(r)∝1r[1+(1/3)exp(−mψr)]

for the gravitational potential outside . For the general case, depending on the choice of , the chamaleon effect (see, e.g., Faulkner et al. 2007 and refs. therein) can take place leading to the same above potential but with a different factor instead of , with depending on both the source mass and the local density which is embedded in.

Motivated by this consideration, we will postulate that the gravitational potential generated by a point mass is :

 ϕ(r)=−Gm(1+δ)r[1+δexp(−rλ)] (4)

where depend on the parameters of the theory. It is worth noticing that a Yukawa - like correction has been invoked several times in the past. For instance, Sanders (1984) showed that the observed flat rotation curves of spiral galaxies may be well fitted by this model with no need for dark haloes provided and is adjusted on a case - by - case basis. As discussed above, Yukawa - like corrections have been obtained, as a general feature, in the framework of -theories of gravity111As hinted at above, for theories showing the chamaleon effect, should be a function of the source mass . We will therefore implicitly assume that all the stars in a galaxy have the same mass. Needless to say, this is far to be true, but, as far as the value of does not change too much with , the impact of this simplifying assumption can be neglected. Alternatively, one should introduce an average over the stellar mass function, but this latter quantity is largely unknown so that we prefer not to make any arbitrary assumption on its functional form and fully neglect the dependence of on . (Capozziello et al., 2009a) and successfully applied to clusters of galaxies setting (Capozziello et al., 2009b). In general, one can relate the length scale to the mass of the effective scalar field introduced by the Extended Theory of Gravity222Referring to -gravity, we prefer to deal with ”Extended Theories of Gravity” and not modified theories of gravity since these theories are nothing else but a straightforward extension of GR where the considered action is a generic function of the Ricci scalar (Capozziello & Francaviglia, 2008; Capozziello & Faraoni, 2010). . The larger is the mass, the smaller will be and the faster will be the exponential decay of the correction, i.e., the larger is the mass, the quicker is the recovering of the classical dynamics333Note that the factor has been explicitly introduced in order to recover the correct Newtonian potential for .. Eq.(4) then gives us the opportunity to investigate in a simple and unified way the impact of a large class of modified gravity theories, among these the Extended Theories, since other details do not have any impact on the galactic scales we are interested in.

Eq.(4) is our starting point for the computation of the rotation curve of an extended system. To this end, we first remember that, in the Newtonian gravity framework, the circular velocity in the equatorial plane is given by , with the total gravitational potential. Thanks to the superposition principle and the linearity of the point mass potential on the mass , this latter is computed by adding the contribution from infinitesimally small mass elements and then transforming the sum into an integral over the mass distribution. For a spherically symmetric body, one can simplify this procedure invoking the Gauss theorem to find out that the usual result with the total mass within . However, because of the Yukawa - like correction, the Gauss theorem does not apply anymore and hence we must generalize the derivation of the gravitational potential444The non-validity of the Gauss theorem is not a shortcoming since, as discussed in Capozziello et al. (2007), physical conservation laws are guaranteed by the Bianchi identities that must hold in any modified theories of gravity.. Alternatively, one can also remember that being the total gravitational force and a radial coordinate. This is the starting point adopted in Cardone et al. (2010) where a general expression has been derived for the case of a generic potential giving rise to a separable force, i.e. :

 Fp(μ,r)=GM⊙r2sfμ(μ)fr(η) (5)

with , and the Solar mass and a characteristic length of the problem. In our case, it is and :

 fr(η)=(1+ηηλ)exp(−η/ηλ)(1+δ)η2 (6)

with . Using cylindrical coordinates and the corresponding dimensionless variables (with ), the total force then reads :

 F(r) = Gρ0rs1+δ (7) × ∫∞0η′dη′∫∞−∞dζ′∫π0fr(Δ)~ρ(η′,θ′,ζ′)dθ′

with , a reference density, and we have defined

 Δ=[η2+η′2−2ηη′cos(θ−θ′)+(ζ−ζ′)2]1/2 . (8)

Since we will be interested in axisymmetric systems, we can set . Moreover, the systems of interest here are spiral galaxies which we will model as the sum of an infinitesimally thin disc and a spherical halo, so that a convenient choice for the scaling radius will be the disc scale length . Under these assumptions, the rotation curve may then be evaluated as :

 v2c(R) = Gρ0R2dη1+δ (9) × ∫∞0η′dη′∫∞−∞~ρ(η′,ζ′)dζ′∫π0fr(Δ0)dθ′

with

 Δ0=Δ(θ=ζ=0)=[η2+η′2−2ηη′cosθ′+ζ′2]1/2 . (10)

Inserting Eq.(6) into Eq.(9) gives rise to an integral which has to be evaluated numerically even for the spherically symmetric case. It is, however, clear that the total rotation curve may be splitted in the sum of the usual Newtonian one and a corrective term disappearing for , i.e. when the extended gravity has no deviations from GR on galactic scales.

## 3 Estimating the bias

Let us assume that a spiral galaxy can be modelled as the sum of an infinitesimally thin disc and a spherical halo and denote with the halo model parameters. We can then write the total rotation curve as :

 v2c(R,Md,pi) = v2dN(R,Md)+v2hN(R,pi) + v2dY(R,Md)+v2hY(R,pi)

where is the disc mass, the labels and denote disc and halo related quantities, while and refer to the Newtonian and Yukawa - like contributions. These latter terms fade away for so that the outer rotation curve is likely the same as the Newtonian one. On the other hand, in the inner region, the two curves may differ more or less depending on the value of . It is worth wondering whether such a difference may be compensated by adjusting the model parameters. That is to say, we are looking for a new set such that :

 v2cN(R,Md,p)=v2dN(R,M′d)+v2hN(R,p′) .

Formally, this problem could be solved explicitly writing down the above relation for values of with the number of halo parameters and then checking that the matching between the two curves is reasonably good (if not exact) along the full radial range. Actually, such a procedure is far from being ideal since it introduces a dependence of the results on the values of chosen. Moreover, we do not need to exactly match the two curves, but only find in such a way that the two curves trace each other within the typical observational uncertainties.

In order to find taking care of typical observations, we therefore adopt the procedure sketched below :

• Compute the theoretical circular velocity for input model parameters using the exact expression.

• Generate a simulated rotation curve by sampling the above and adding noise to the extracted points.

• Fit the simulated curve with the same disc + halo model, but using the Newtonian theory to compute .

• Compare the output parameters with the input ones as function of the normalized scale length of the modified gravity theory and other quantities of interest.

In the following, we describe in more details steps i. and ii., while Sections 4 and 5 are devoted to discussing the results from a large sample of simulated curves.

### 3.1 Modeling spiral galaxies

As yet said above, we will model a spiral galaxy as the sum of a thick disc and a spherical halo. For the disc, we adopt a double exponential disc so that the density profile reads :

 ρd(R,z)=Md4πR2dzdexp(−RRd−|z|zd) (11)

where are the disc mass, lengthscale and heightscale. The Newtonian rotation curve cannot be computed analytically, but, if (as we will assume), it is well approximated by the formula for the infinitesimally thin disc (Freeman, 1970; Binney & Tremaine, 1987) :

 v2dN(R)=(2GMd/Rd)y2[I0(y)K0(y)−I1(y)K1(y)] (12)

with and , the modified Bessel functions of order of the first and second kind, respectively.

While the observed photometry motivates the use of the exponential profile for the disc, the choice of the dark halo model is not trivial. Numerical simulations of structure formation are typically invoked as a direct evidence favouring the use of the NFW density law or its variants. However, the NFW model is the outcome of DM only simulations performed in a Newtonian framework, while here we are working in a modified gravity theory. In principle, one should therefore rely on the results of simulations which include both the effect of the different potential and the impact on the evolution of structure due to deviations from GR. To this end, one has first to specify which is the modified gravity theory one is considering, i.e., explicitly write down the gravity Lagrangian. In the case of -gravity, Schmidt et al. (2009) have shown that, provided satisfies the constraints from the Solar System, the halo density profiles of DM haloes from numerical simulations is still well approximated by the NFW model over the range of masses of interest here. Motivated by this result, we therefore assume that the DM density profile is the NFW one :

 ρh(r)=Mvir4πR3sg(Rvir/Rs)(rRs)−1(1+rRs)−2 (13)

with

 g(x)=ln(1+x)−x/(1+x) . (14)

In Eq.(13), and are the virial mass and radius. They are not independent being related by

 Rvir=(3Mvir4πΔth¯ρM)1/3

with the overdensity for spherical collapse and the mean matter density today. We follow Bryan & Norman (1998) for and set in accordance with Komatsu et al. (2010).

Because of the spherical symmetry, the Newtonian rotation curve may be easily evaluated as :

 v2hN(r)=GMh(r)r=GMvirRvirg(r/Rs)g(Rvir/Rs) . (15)

While the Newtonian contributions to the rotation curve may be computed in the usual way, the Yukawa - like terms in the potential give rise to two further terms in that have to be computed numerically. To this end, we must only insert Eqs.(11) and (13) into Eq.(9) with given by Eq.(6) subtracting the first term in parentheses. Some care must be taken in choosing the reference radius . Since the data typically probe a limited range in , a natural choice is to set . However, the reference density is not the density at . It is indeed more convenient to set for the disc and for the halo. With such a choice, for the halo, the dimensionless density profile entering Eq.(9) is given by the  - dependent part of Eq.(13) provided is replaced everywhere by . Finally, we remember the reader that, when computing the total rotation curve, the Newtonian terms, for both the disc and the halo, given by Eqs.(12) and (15) respectively, must be rescaled by the factor in order to recover the classical results in the GR limit ().

### 3.2 Simulating the rotation curve

In order to be useful, our approach should rely on simulated rotation curves that are as realistic as possible. By this, we mean that a.) they must refer to spiral galaxies with reasonable values of the model parameters and b.) the sampling and the noise should be the same as actual data. Point a.) is the easiest to address. As a first step, we randomly generate the disc scalelength and the halo virial mass from flat distributions over the ranges :

 0.5≤Rd/Rd,MW≤2.0  ,  11.5≤logMvir≤13.5  ,

with the disc scalelength for the Milky Way (Dehnen & Binney, 1998; Cardone & Sereno, 2005). In order to set the disc mass, we first define the halo mass fraction within the optical radius as :

 fDM=Mh(Ropt)Md+Mh(Ropt) (16)

with the optical radius and we have approximated the disc mass within with the total disc mass. We then randomly generate from a flat distribution in the range and (see, e.g., Williams et al. (2010) and refs. therein). The halo scalelength is computed as where the concentration is randomly generated from a Gaussian centred on :

 c=16.7(Mvir1011 h−1 M⊙)−0.125 (17)

and variance set to of the mean. Note that the above relation has been derived by Napolitano et al. (2005) for the mass range following the method detailed in Bullock et al. (2001) and updating the cosmological model. Finally, we need to set the modified potential parameters . Following the result obtained for -gravity (Capozziello et al., 2009b), we first set and run different simulations randomly generating from a flat distribution covering the wide range . In order to explore the impact of , we also consider the extremal case thus maximizing the contribution of the Yukawa terms.

Having thus generated a realistic galaxy model, we now need a rule for sampling it and adding noise in such a way that the simulated rotation curve is similar to an observed one. A unique choice is not possible since the details of any observation depend on the instrumental setup, the observing conditions and the galaxy surface brightness. After a visual examination of rotation curves samples in literature, we have adopted the strategy summarized below.

1. Take equally spaced points in the range . and replace each of them with with a randomly generated from the range (0.9, 1.1).

2. For each point in the sample, generate from a Gaussian distribution centred on the theoretical value and with variance set to .

3. Set the error on the  - th point as with randomly chosen in the range .

4. Generate two random numbers in the range and take the point if .

Typical rotation curves are well sampled up to the optical radius , while the sampling gets worse at larger radii. On the contrary, the percentage errors are of the same order in both regions, although it is possible that the innermost points have larger uncertainties because of deviations from ordered motions. For given input model parameters, we therefore generate two samples of points by first setting

 (Nsim,ηmin,ηmax,εc)=(30,0.1,3.1,0.25)

and then

 (Nsim,ηmin,ηmax,εc)=(20,3.0,10.0,0.20) .

We then add the two samples and rescale the errors in such a way that the standard equals 1 for the input rotation curve. Although such values are arbitrary, we have checked that the simulated rotation curves have a similar sampling and uncertainties of many dataset in literature (see Fig. 1) so that we will not try other possible combinations. In order to quantify the bias on the model parameters, we then simulate 1000 rotation curves and fit different models (as described in the following) to determine their best fit parameters. Note that we will not consider the errors on the fit quantities since they strongly depend on the uncertainties on the data points hence on the choice of the simulation parameters . Since we are interested in a statistical analysis of the full sample rather than on a detailed study of a particular case, we prefer to limit our attention to the best fit parameters only thus avoiding to deal with the details of the simulation procedure.

## 4 Bias on halo parameters

The sample of simulated rotation curves constructed by the procedure detailed above is the starting point of our analysis. Indeed, we can now fit each curve with a given (Newtonian) model and compare the output best fit parameters with the input ones. There is, however, a preliminary caveat to be addressed. When fitting a whatever model to a given dataset, one has to put down an objective criterium to deem the model as reliable or not. It is then common practice to compute the standard , i.e.

 χ2=N∑i=1[vc(Ri)−vc,th(Ri)σi]2 ,

and then evaluate the quality of the fit considering the value of with the number of degrees of freedom and () the number of data points (parameters). Formally, one can then conclude that the model is a good description of the data if with the threshold value depending on . Since, for our simulated curves , which is highly demanding selection criterium. Actually, this formal rule rigorously applies only if the data points are uncorrelated and the measurement uncertainties are Gaussian distributed. Both these assumptions break down for our simulated curves since we have rescaled the uncertainties in such a way that for the input model. By a visual examination of many fitted curves, we have checked that a reasonably good fit (with no systematic deviations from the outer or inner data and rms percentage deviations smaller than ) is obtained up to which we choose as our cut. We will therefore consider a model as successfully fitting the data if and with this latter criterium introduced to avoid unphysical solutions. However, we will also discuss the trends for the quantities of interest for a subsample obtained by using a stringent criterium, i.e. in order to explore the impact of the threshold value.

### 4.1 Navarro-Frenk-White models

In realistic situations, one has a set of data and a measurement of the galaxy surface brightness in a given band. It is then common to set the disc scalelength to the value inferred from photometry, while the disc mass can be inferred from the total luminosity provided an estimate of the stellar ratio is somewhat available (e.g., from stellar population synthesis models fitted to the galaxy colours). As a first step, we therefore assume that both are known and fit the simulated data for each rotation curve to determine the NFW model parameters555To this end, we use a straightforward Monte Carlo Markov Chain algorithm to minimize the with respect to the parameters and then infer the values of and the dark matter mass fraction . This approach is more stable than fitting directly for and allows to correctly recover the input model parameters when we fit rotation curves simulated in the Newtonian framework. only.

#### 4.1.1 δ=1/3 simulated samples

Let us first consider the results obtained setting when simulating the rotation curves. As yet said above, we are here trying to fit modified gravity rotation curves using theoretical models computed in a Newtonian framework so that the first point to address is whether this is indeed possible. To this end, it is sufficient to check what is the percentage of rotation curves passing the two selection criteria quoted before. We indeed find that of the simulated rotation curves are fitted with and (which we will refer to as the well fitted sample, hereafter WS), while this fraction decreases to if is demanded (defining what we will hereafter refer to as the best fitted sample, BS). Averaging over the sample values gives :

 ⟨~χ2⟩=1.27±0.24 (1.19±0.13)  for WS (BS)

while for the rms of the percentage deviations we get :

 ⟨rms(Δvc/vc)⟩=6.4%±0.8% (6.3%±0.7%) for WS (BS) .

Motivated by the high fraction of successful fits and the low values of and , we can therefore safely conclude that it is possible to fit a modified gravity rotation curve with a NFW Newtonian one obtaining both an acceptable value and reasonable model parameters.

Having checked that the model reasonably well fit the data, we can rely on the best fit model parameters and compare them with the input ones to estimate what is the bias induced by an incorrect assumption of the gravity law. To quantify this bias, we define , i.e. the ratio between the best fit and the input values of a given quantity . Needless to say, should be on average close to 1, we can conclude that assuming Newtonian gravity to fit modified gravity rotation curves does not induce a significant bias on the estimate of . The results, summarized in Table 1 for both the WS and BS samples, show that such a bias is indeed quite important and only mildly depend on the selection criteria adopted (so that we will hereafter refer to the WS results only unless otherwise stated).

Table 1 gives some statistics on the distribution of for the halo parameters . However, these values could be misleading since the histograms for and are strongly asymmetric with long tails towards the right. In other words, and could be much larger than their median value with being as high as 12 and values of as large as 15 so that both the halo scalelength and virial mass can be grossly overestimated. Since , one could naively expect that the concentration is overestimated too. On the contrary, the distribution of the values is almost symmetric around its mean clearly disfavouring , i.e. the concentration is underestimated. Actually, such a result can be understood remembering that so that finally leading to values smaller than unity. The emerging picture is therefore that of a halo having a larger mass than the input one, but also a larger scalelength. This can be qualitatively explained noting that the Yukawa - terms in the rotation curve increases the net circular velocity both in the inner and outer regions probed by the data. In order to adjust the fit in the halo dominated regions, one has to increase , but then has to be increased to in order to lower the density (and hence the contribution to the rotation curve) in the inner regions not to overcome the observed circular velocity. As a result, the concentration is underestimated and the same takes place for the dark matter mass fraction within the optical radius since the halo mass is now pushed outside the optical radius thus explaining why is smaller than unity. As Fig. 2 shows, is correlated with with the sign of the correlation depending on being typically larger or smaller than 1. Not surprisingly, the larger is the Yukawa scalelength , the closer is to 1, i.e., the smaller is the bias induced by the assumption of Newtonian gravity. However, there is a large scatter in these correlations since the relative importance of the correction also depends on the input halo parameters. Indeed, when , the halo contribute to the modified rotation curve is maximized so that the Yukawa corrections are more important and the bias of the fitted parameters is larger. As a consequence, one can therefore expect a correlation between and which is what we indeed find looking at the results in Table 1.

#### 4.1.2 δ=1.0 simulated samples

Let us now consider the results for the rotation curves sample simulated under the extreme assumption , i.e. when the contribution of the Yukawa term to the point mass potential (4) is the same order of magnitude as the Newtonian one. Note that, different from the previous case, is, as far as we know it, an unmotivated assumption, but we consider it to investigate the impact of on the results.

The first striking result is that the fraction of well fitted rotation curves dramatically drops to applying the WS selection criteria ( and ), while no curves pass the cut thus showing that it is likely not possible to reproduce the simulated curves with NFW model if Newtonian gravity is assumed. Such a result can be understood noting that, for , the modified gravity rotation curve is much larger than the Newtonian one so that one has to greatly increase the halo virial mass to fit the outer rotation curve. In order to compensate the corresponding increase in the inner region, one should set as large as possible, but this also tend to lower the circular velocity for . Finding a compromise between these two opposite trends is quite difficult thus leading to unacceptably large values.

It is worth investigating what is the bias induced on the halo model parameters for the successfully fitted cases. The trends of the bias parameters are the same as in the previous case, but the typical values are now much larger. In particular, can be as high as 40 thus leading to grossly underestimates of both the concentration and the dark matter mass fraction . Not surprisingly, the only way to reduce such large biases is to increase to values larger than our upper limit . However, in this case, the modified gravity potential differs from the Newtonian one only on group and cluster scales where the rotation curve are no more measured so that our analysis becomes meaningless.

### 4.2 Generalized Navarro-Frenk-White models

Up to now, we have considered the NFW profile for the halo model, but such a choice does not allow us to investigate whether the mismatch between modified and Newtonian gravity can have an impact on the inner logarithmic slope. In order to explore this issue, we therefore consider the generalized NFW model (gNFW, Jing & Suto 2000) :

 ρh(r)=Mvir4πR3sh(Rvir/Rs)(rRs)−α(1+rRs)−(3−α) (18)

with

 h(x)=∫x0ξ2−αdξ(1+ξ)3−α . (19)

The parameter depends on the behaviour of the rotation curve in the inner regions where the disc contribute can be larger than the halo one. It is therefore worth investigating which is the impact of the uncertainties on in order to see how the bias on depend on it.

#### 4.2.1 gNFW models with known disc mass

As a first step, we assume that is set to the input value and only look at the bias on the parameters . As a preliminary remark, it is mandatory explaining how the concentration is defined for gNFW models. To this end, one may simply note that, for the NFW density profile, is the radius at which the logarithmic slope equals the isothermal value . Therefore, the concentration for the gNFW model may be easily defined as :

 c=RvirR−2=Rvir(2−α)Rs (20)

with . Note that, for , the gNFW model reduces to the NFW one and Eq.(20) gives the standard definition for the concentration of NFW haloes.

Applying the selection criteria used above to the sample with , we find that the WS sample contains now of the simulated curves, while this fraction lowers to for the BS sample. Averaging over the samples gives :

 ⟨~χ2⟩=1.28±0.19 (1.23±0.12)  for WS (BS) ,
 ⟨rms(Δvc/vc)⟩=6.4%±0.8% (6.3%±0.8%) for WS (BS) .

These values are almost identical to those obtained for the fits of the NFW model described before and allows us to confidently conclude that it is possible to reproduce the simulated rotation curves with a gNFW model and Newtonian gravity. In a sense, this is not surprising since the gNFW model has one more parameter than the NFW one so that there is a much larger freedom to adjust the shape of the theoretical rotation curve to fit the observed one. As a result, the percentage of well fitted rotation curves slightly increases, but the quality estimator ( and the rms value of ) stay almost unchanged.

The bias induced on the model parameters can be still estimated considering the quantity defined before and summarized in Table 2 for the present case666Note that, for all the simulated curves, it is since we have assumed a NFW model in the simulation process. Therefore, for the gNFW models.. Comparing to the results in Table 1, it is apparent that the bias induced on the halo parameters in common, i.e. , is almost the same both for the WS and BS samples. The same qualitative discussion can be repeated here so that we will not consider anymore this issue. It is, on the contrary, more interesting to look at the bias on which asks for some caution. As can be seen from the left panel in Fig. 3, the distribution of values is quite large and asymmetric, while the right panel is a clear evidence that its mean value strongly depend on . Indeed, should we have only fitted models with , would have been smaller than 1, i.e. the inner slope would have been underestimated and shallower models be preferred. Should the modified gravity scalelength be smaller than the typical disc one, fitting rotation curves assuming the validity of Newtonian gravity would have us led incorrectly to believe that the inner slope of the density profile is much smaller than the NFW one. Such a result, therefore, suggests that the cusp/core controversy is just a fake problem originated by an incorrect assumption of what the underlying gravity theory actually is.

As a final remark, it is worth noting that a second difference among the results in Tables 1 and 2 is represented by the markedly different values of the Spearman correlation coefficients. Actually, this is a consequence of having added one more parameter which introduces a degeneracy between the fitted quantities hence partially washing out the correlations of with the input halo parameters and the modified gravity scalelength. In particular, and are clearly correlated since both set the shape of the inner rotation curve. Although with a large scatter, is typically smaller (larger) than 1 for () so that actually correlates with , but not linearly thus giving rise to a small Spearman correlation coefficient (which looks for linear correlations). Since has to be adjusted according to the changes in , has to follow thus losing its correlation with .

Such results are strengthened if we consider the simulated curves for , but now relying on a much smaller statistics. Indeed, the WS sample is now made out of only of the full set of simulated curves, while only of them are included in the BS sample. These fractions are larger than in the NFW fits considered before, but are still so small to make us conclude that such extreme modified gravity curves can not be reproduced in a Newtonian framework using the gNFW model. For the few surviving curves, we can repeat the same discussion above with the only caveat that now the values deviate more from 1 than in the case. Moreover, the slope can be underestimated also when depending on the input halo parameters. This is a still further evidence in favour of the suggested interpretation of the cusp/core controversy as the outcome of a systematic error in the assumed gravity theory.

#### 4.2.2 gNFW models with free disc mass

Up to now, we have set the disc mass to the input value implicitly assuming that one is able to perfectly estimate this quantity from the galaxy colors and a stellar population synthesis code. Actually, this is not always the case and, on the contrary, it is usual to include in the set of parameters to be determined by fitting the rotation curve data. We have therefore repeated the above analysis for the gNFW model leaving free to be adjusted by the fit.

Considering there is one more fit parameter, it is not surprising that the fraction of successfully fitted rotation curves (for ) increases to for the WS (BS) samples and also the quality of the fit is improved being :

 ⟨~χ2⟩=1.19±0.20 (1.14±0.13)  for WS (BS) ,
 ⟨rms(Δvc/vc)⟩=6.1%±0.6% (6.0%±0.8%) for WS (BS) .

Table 3 summarizes the results on the bias estimator for the different parameters involved which allow us to draw some interesting considerations. First, we note that the halo model parameters are now better recovered than in previous cases. For instance, although the corresponding distributions have long tails up to , most of the values of and are smaller than 2 leading to a modal value (i.e., the most peak of the histogram) close to 1, i.e. are not biased anymore for most of the WS and BS sample curves. A similar discussion also apply to the concentration and the dark matter mass fraction whose histograms are quite symmetric and centred close to the no bias value, . Concerning , one must still keep in mind that its mean value depend on the cut on typically being (i.e., the inner slope is shallower than the input one) when . Note, however, that this effect is now less pronounced than before thus explaining why the distribution of values is peaked close to 1 with a not too large width. On the contrary, the disc mass turns out to be biased low with being smaller than the input value for most of the cases. Such an unexpected result can be qualitatively understood considering that the Yukawa terms make the modified gravity rotation curve larger than the Newtonian one in the inner and outer regions. As in the case of the NFW fit, one must increase the halo virial mass to recover this additional contribution in the halo dominated regions which increases also the halo inner circular velocity. In the NFW case, one has then to increase to prevent overestimating the simulated curve, while the same effect is now accomplished by lowering the disc mass. This explains why the halo parameters turn out to be less biased than in the NFW fit case. It is, however, worth stressing that, as shown in the left panel of Fig. 4, the distribution is actually bimodal. As a consequence, when the disc mass is not biased low, one has again to increase thus explaining the long tails towards values. Alternatively, one can leave unchanged, but make the density profile shallower hence giving rise to the tail towards small in the right panel of Fig. 4. Model degeneracies also explain why the correlation coefficients among and the input model parameters are typically quite small. We, however, caution the reader that a small value of Spearman correlation does not necessarily imply that doe not depend on the corresponding quantity (see, e.g., the dependence on for the gNFW fits with fixed disc mass). We have therefore checked whether this is the case or not by directly looking at the vs plots (with the  - th input parameter) finding that the scatter of the points in each plot clearly indicates a complicated dependence a weak dependence on all the input parameters but .

Finally, we have repeated the above analysis setting to generate the simulated rotation curves. Differently from the previous cases considered up to now, we now find that it is possible to fit the gNFW model with free disc mass to the simulated rotation curves. In a sense, the additional freedom we have now to adjust the disc mass makes it possible to recover the inner rotation curve better than in the previous models thus leading to a successful fit for () of the galaxies using the selection criteria for sample WS (BS). The quality of the fit is also remarkably good :

 ⟨~χ2⟩=1.32±0.30 (1.18±0.15)  for WS (BS) ,
 ⟨rms(Δvc/vc)⟩=6.5%±1.0% (6.2%±0.7%) for WS (BS) .

Concerning the values of quantities, they are almost the same as those in Table 3 for the parameters which are therefore almost unbiased. On the contrary, for the other parameters (and the WS sample), we get :

 ⟨R(x)⟩=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩0.54±0.20  for  x=Md1.04±0.28  for  x=α0.85±1.18  for  x=Mvir ,
 [R(x)]rms=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩0.58  for  x=Md1.08  for  x=α1.5  for  x=Mvir .

As it is clear from the large standard deviations and rms values, the distribution of and are strongly asymmetric and, as a result, the one for presents a non negligible amount of points with . This result can be explained noting that setting maximizes the contribution of the Yukawa terms so that we now have to mimic a similar effect by adjusting the model parameters in the Newtonian fitting. One possibility is to leave unchanged with respect to the input value , but increasing the halo mass by a significant factor thus leading to large values of . But, in this case, one has also to decrease the disc contribution to the inner rotation curve to compensate for the additional dark matter. On the contrary, one can leave almost unchanged , but redistribute the dark mass pushing it towards the inner region. To this aim, one should make the profile steeper, i.e. increase , while must be smaller to not overcome the inner circular velocity. This strategy will lead to and , while again it is . The final histograms of for is then the outcome of a mixture of these two possible strategies whose relative contribution depends on the details of the individual rotation curves.

## 5 Bias on the halo density profile

The use of the gNFW model to fit the simulated rotation curve may also be read as a first step towards completely mismatching the halo density profile. In a sense, the gNFW model departs from the input NFW one only because of the possibly different inner logarithmic slope. It is, however, common in data analysis to use also models that differs from the NFW one both in the inner and outer regions. Moreover, the gNFW is still a cuspy model, while the cusp/core controversy comes from the observation that cored models better fit the observed rotation curves. To this end, we now drop the implicit assumption that the trial model tracks or generalizes the NFW profile and investigate whether completely different models in the Newtonian framework are in accordance with our modified gravity simulated rotation curves.

### 5.1 The pseudo - isothermal sphere

A classical example of a cored model is represented by the pseudo - isothermal (hereafter, PI) model whose density profile reads :

 ρh(r)=Mvir4πR3shpi(Rvir/Rs(1+rRs2)−1 (21)

with

 h(x)=x−arctanx . (22)

As it is clear, the PI density law approaches a constant value for , while falls off as in the outer regions. As such, the PI model is radically different from the NFW and gNFW ones so that it is interesting to see whether can fit or not the data. Following common practice, we leave the disc mass as an unknown and estimate .

Applying the corresponding selection cuts, we find that of the simulated galaxies fall into the WS sample, while the BS one contains of the full sample. Such high fraction and the low values of both the reduced

 ⟨~χ2⟩=1.49±0.35 (1.24±0.15)  for WS (BS) ,

and of the rms of percentage residuals

 ⟨rms(Δvc/vc)⟩=7.0%±0.9% (6.9%±0.9%) for WS (BS)

make us confident that the Newtonian circular velocity of the PI model can indeed mimic the rotation curve predicted by modified gravity. This is in agreement with what is indeed found in the literature where galaxies rotation curves data are typically well fitted by PI haloes.

Having used now a different halo model to fit the data, it is not possible to straightforwardly compare the output parameters with the input ones. First, although we use the same symbol, the scalelength of the PI model is not defined by the condition as for the NFW model, but rather represents the core radius. As such, it is meaningless to compare the input with the output since they refer to a different characteristic of the density profile. Moreover, a concentration for the PI model is not defined so that we can not compare with the input one. On the contrary, the total disc mass , the virial mass and the dark matter mass fraction refer to global properties of the disc + halo model so that it makes sense to compare them with the input ones. Table 4 then gives the bias on the global parameters both the WS and BS samples (having set ). It is clear that all these three quantities are severely biased with being overestimated and grossly underestimated. Such a behaviour can be explained noting that, once the core radius has been set, the only way to increase the Newtonian rotation curve is to increase the global masses thus leading to asymmetric distributions both peaked in . It is worth noting that, if one ignores that the underlying gravity theory is not Newtonian, a large disc mass will be interpreted as the presence of a maximal disc which is indeed the case in many spiral galaxies (Palunas & Williams, 2000). Being both biased high, one could find the result an unexpected nonsense. Actually, one has first to note that also the disc mass is overestimated which automatically reduces the dark matter mass fraction. Moreover, for the PI model, the mass profile increases almost linearly, while, for the input NFW profile, grows approximately in a logarithmic way. Being the input almost the same as the output core radius, we therefore get