Dark Matter in Fractional Gravity. I. Astrophysical Tests on Galactic Scales

A multitude of astrophysical and cosmological probes have firmly established that baryons constitute only some 15% of the total matter content in the universe, the rest being in the form of a dark matter (DM) component. A nonexhaustive list of such evidence includes the kinematics of spiral galaxies (e.g., Rubin et al. 1980; Persic et al. 1996; see review by Salucci 2019 and references therein); dynamical modeling of ellipticals (e.g., Bacon et al. 2001; Cappellari et al. 2013; Mendel et al. 2020); cosmic microwave background and Big Bang nucleosynthesis constraints (e.g., Bennett et al. 2003; Ade et al. 2014; Aver et al. 2015; Aghanim et al. 2020); baryon acoustic oscillation measurements (e.g., Eisenstein et al. 2005; Beutler et al. 2011; Zhao et al. 2022); type Ia cosmography (e.g., Perlmutter et al. 1999; Scolnic et al. 2018; Brout et al. 2022); cosmic shear galaxy surveys (e.g., Heymans et al. 2013; Secco et al. 2022; Amon et al. 2022); X-ray, Sunyaev–Zel'dovich, and strong/weak lensing observations and number counts of galaxy clusters (e.g., White et al. 1993; Allen et al. 2011; Umetsu et al. 2020; Garrel et al. 2022; Mantz et al. 2022); and observations of the "Bullet Cluster" (e.g., Markevitch et al. 2004; Clowe et al. 2006; Paraficz et al. 2016). Incidentally, note that most of the baryons in the universe are locked in a hot/warm intergalactic medium, with only a small fraction (less than about 20%) residing as cold gas and stars in galaxies (e.g., Tumlinson et al. 2017 and references therein).

Still, no firm detection of DM particles has been made so far, despite the big efforts carried on with colliders (see Mitsou 2013; Argyropoulos et al. 2021) or with direct (see Aprile et al. 2018; Bernabei et al. 2020) and indirect (see Ackermann et al. 2017; Zornoza 2021) searches in the sky. The standard picture envisages DM to be constituted of weakly interacting particles with GeV masses (see Bertone & Hooper 2018) that are nonrelativistic at the epoch of decoupling (hence, they are dubbed "cold" dark matter, or CDM) and feature negligible free-streaming velocities (i.e., they do not diffuse out of perturbations before collapse). As a consequence, bound CDM structures called halos grow sequentially in time and hierarchically in mass by stochastically merging together (e.g., Frenk & White 2012; Lapi & Danese 2020, Lapi et al. 2022a).

Although on cosmological scales, the CDM hypothesis is remarkably consistent with the data, on (sub)galactic scales, it faces some severe challenges. For example, with respect to the predictions of gravity-only N-body simulations, the shape of the inner density profiles inferred from the rotation curve (RC) in DM-dominated dwarfs is too flat (e.g., Flores & Primack 1994; Gentile et al. 2004; de Blok et al. 2008; Oh et al. 2015; Salucci et al. 2021), and the number and dynamical properties of observed Milky Way satellites differ from those of subhalos (see Boylan-Kolchin et al. 2012; Bullock & Boylan-Kolchin 2017). Moreover, the emergence of tight empirical relationships between the properties of the dark and luminous components in disk-dominated galaxies, such as the radial acceleration relation (RAR; see McGaugh et al. 2016; Lelli et al. 2017; Di Paolo et al. 2019), the universal core surface density (see Donato et al. 2009), and the scaling between the core radius and the disk scale length (see Donato et al. 2004), are puzzling and seem to be indicative of a new dark sector and/or nongravitational coupling between DM particles and baryons (see Salucci et al. 2020). Some of the above findings can, in principle, be explained in CDM by invoking physical processes causing a transfer of energy and angular momentum from the baryons to DM particles, such as dynamical friction (see El-Zant et al. 2001; Tonini et al. 2006) or feedback effects from stars and active galactic nuclei (see Pontzen & Governato 2014; Freundlich et al. 2020), but in recent years, it has become progressively clearer that a fine-tuning of these solutions is required to explain in detail the current observations.

This has triggered the consideration of alternative, and perhaps more fascinating, explanations that rely on nonstandard particle candidates (see Bertone et al. 2004; Adhikari et al. 2017; Salucci et al. 2021); the most widespread possibilities include keV-scale warm DM particles (see Bode et al. 2001; Lovell et al. 2014), fuzzy or particle-wave DM (see Hu et al. 2000; Hui et al. 2017), self-interacting DM (see Vogelsberger et al. 2016), and dark photon DM (McDermott & Witte 2020; Bolton et al. 2022). As a consequence of free streaming, quantum pressure effects, and/or dark sector interactions, all of these scenarios produce a matter power spectrum suppressed on small scales, fewer (sub)structures, and flatter inner density profiles within halos relative to CDM. However, all of these possibilities, despite being still viable to some extent, have been strongly constrained by various astrophysical probes, most noticeably the Lyα forest (e.g., Irsic et al. 2017), high-redshift galaxy counts (e.g., Shirasaki et al. 2021), cosmic reionization (e.g., Lapi et al. 2022b), gravitational lensing (e.g., Ritondale et al. 2018), integrated 21 cm data (e.g., Carucci et al. 2015), γ-ray emission (e.g., Grand & White 2022), fossil records of the Local Group (e.g., Weisz & Boylan-Kolchin 2017), dwarf galaxy profiles and scaling relations (e.g., Burkert 2020), Milky Way satellite galaxies (e.g., Newton et al. 2021), cosmic star formation rate density (e.g., Gandolfi et al. 2022b), or a combination of these.

All that has motivated, on different grounds, the emergence of modified gravity theories. Perhaps the most famous one is the modified Newtonian dynamics (MOND) that was originally proposed by Milgrom (1983) and further investigated in a rich literature (see Bekenstein 2004; Bruneton & Esposito-Farese 2007; Famaey & McGaugh 2012). As the name suggests, MOND aims to explain galactic dynamics through a modification of Newtonian gravity (or, more generally, Newton's second law) that comes into action at accelerations well below a definite universal threshold; in its original formulation, DM is not included, and baryons are the only source of the gravitational field. In fact, it has been claimed that MOND (or theories reducing to it in the weak-field limit) can properly fit galactic RCs (de Blok & McGaugh 1998; Sanders & McGaugh 2002) and provide a satisfying description of the RAR (e.g., Li et al. 2018). On the other hand, the performances of MOND are highly debated when moving to galaxy clusters and cosmological scales (see Aguirre et al. 2001; Angus et al. 2006; Scott et al. 2001; McGaugh 2015; Nieuwenhuizen 2017; Boran et al. 2018). Other more recent proposals include, e.g., an emergent (entropic) theory of gravity (see Verlinde 2017; Yoon et al. 2023) or a dynamical nonminimal coupling between the DM component and the gravitational metric (see Gandolfi et al. 2021, 2022a).

In this paper, we explore yet another possibility that has some features in common with modified gravity models, though it is fundamentally different. Specifically, we propose that DM in galaxies originates fractional gravity; this means that the gravitational potential associated with a given DM density distribution is determined by a modified Poisson equation that includes fractional derivatives (i.e., derivatives of noninteger type) meant to describe nonlocal effects. Fractional calculus has found applications in the descriptions of several physical phenomena that subtend nonlocal behavior in space and time (e.g., material science, rheology, seismology, fluid dynamics, nuclear physics, medicine, finance, etc.). In astrophysics and cosmology, the related investigations are still scanty, though their number and variety have been progressively increasing in recent years (e.g., Calcagni 2010, 2013; Varieschi 2020; Giusti 2020; Giusti et al. 2020; Varieschi 2021; Calcagni 2021; Varieschi 2022; Calcagni & Varieschi 2022; García-Aspeitia et al. 2022; Varieschi 2023; see also Section 4).

Here in particular, we aim at testing whether fractional gravity as originated by the DM component can solve the aforementioned problems of the standard paradigm in galaxies while retaining its successes on cosmological scales. We stress that in our fractional gravity framework, the law of inertia is unchanged, while only the Poisson equation is modified. Note that fractional gravity is not necessarily meant to be an ab initio theory but may constitute an effective description for a whole class of models that imply the development of nonlocal effects in the gravitational behavior of the DM component. As for the baryons, we assume that they originate standard gravity, but feel the overall gravitational potential of the system.

We will show that in fractional gravity, the dynamics of a test particle in the potential sourced by a standard Navarro–Frenk–White (NFW) density distribution (that we consider as a representative rendition of the density distribution for halos of noninteracting DM) is substantially altered with respect to the Newtonian case. Specifically, such a framework can (i) provide accurate fits to the stacked RCs of spiral galaxies with different properties, including dwarfs and high/low surface brightness systems; (ii) reproduce to reasonable accuracy the observed shape and scatter of the RAR over an extended range of galaxy accelerations; and (iii) properly account for the universal surface density and the scaling relation between the core radius of the DM component and the disk scale length.

The plan of the paper is as follows. In Section 2, we introduce the fractional gravity framework and describe how this can be applied to fit the stacked RCs of spiral galaxies and infer related physical parameters via a Bayesian Markov Chain Monte Carlo (MCMC) framework. In Section 3, we present and discuss our results, with particular focus on the RC of dwarf galaxies, the halo to (stellar) disk mass relation, the RAR, and the universal core surface density. In Section 4, we explore an alternative MONDian viewpoint to fractional gravity without DM. In Section 5, we summarize our findings and outline future perspectives and applications of the fractional gravity framework. Appendix A contains a primer on fractional calculus, including the fractional Laplacian and Poisson equation. In Appendix B, we provide original analytic solutions of the gravitational potential in fractional gravity for various spherically symmetric density distributions.

Throughout this work, we adopt the standard flat ΛCDM cosmology (see Aghanim et al. 2020) with rounded parameter values: matter density Ω_M = 0.3, dark energy density Ω_Λ = 0.7, baryon density Ω_b = 0.05, and Hubble constant H₀ = 100 h km s⁻¹ Mpc⁻¹ with h = 0.7. Unless otherwise specified, G ≈ 6.67 × 10⁻⁸ cm³ g⁻¹ s⁻² indicates the standard gravitational constant.

As anticipated in Section 1, in this paper, we propose that the kinematics observed in galactic structures may be explained by assuming that the DM component originates fractional gravity. In such a framework, although the standard law of inertia continues to hold, the gravitational potential associated with a given DM distribution is determined by a modified Poisson equation of the fractional type. The latter has an inherently nonlocal nature and can profoundly alter the dynamics of a test particle with respect to the Newtonian case (i.e., based on the standard Poisson equation).

In this section, we first present our new solution for the potential sourced by the NFW density profile in fractional gravity and then describe how this can be applied to analyze stacked RCs of spiral galaxies via a Bayesian MCMC technique, thus inferring relevant galaxy properties and at the same time constraining the parameters ruling the strength and length scale of fractional gravity.

2.1. DM in Fractional Gravity

The density distribution of virialized halos for collisionless DM is routinely described via the NFW profile (see Navarro et al. 1997),

$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{\rho }_{s}\,{r}_{s}^{3}}{r{\left(r+{r}_{s}\right)}^{2}},\end{eqnarray} \tag{ 1 }$

where r_s is a scale radius, and ρ_s is a characteristic density. The logarithmic density slope $\gamma \equiv d\mathrm{ln}\rho /d\mathrm{ln}r$ takes on values around −3 in the outskirts, −2 at r ≈ r_s , and −1 in the inner region. The associated cumulative mass is given by

$\begin{eqnarray}\begin{array}{rcl}M(\lt r) & = & 4\pi \,{\int }_{0}^{r}{dr}^{\prime} \,{r}^{{\prime} 2}\,\rho (r^{\prime} )\\ & = & {M}_{s}\left[\mathrm{ln}\left(1+\displaystyle \frac{r}{{r}_{s}}\right)-\displaystyle \frac{r/{r}_{s}}{1+r/{r}_{s}}\right],\end{array}\end{eqnarray} \tag{ 2 }$

where ${M}_{s}\equiv 4\pi \,{\rho }_{s}\,{r}_{s}^{3}$.

In Newtonian gravity, the potential Φ(r) associated with a given density distribution ρ(r) is computed from the Poisson equation supplemented with appropriate boundary conditions (usually taken as a vanishing potential at infinity),

$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\Phi }}({\boldsymbol{r}})=4\pi G\,\rho ({\boldsymbol{r}}),\end{eqnarray} \tag{ 3 }$

where Δ is the Laplacian operator. For the spherically symmetric NFW profile, it is easily found that

$\begin{eqnarray}&&{\rm{\Phi }}(r)=-\displaystyle \frac{{{GM}}_{s}}{r}\,\mathrm{log}\left(1+\displaystyle \frac{r}{{r}_{s}}\right).\end{eqnarray} \tag{ 4 }$

The rotational velocity, or, in other words, the RC, v²(r) = r ∣dΦ/dr∣, of a test mass in such a potential is just

$\begin{eqnarray}&&{v}^{2}(r)=\displaystyle \frac{{{GM}}_{s}}{{r}_{s}}\left[\displaystyle \frac{\mathrm{ln}(1+r/{r}_{s})}{r/{r}_{s}}-\displaystyle \frac{1}{1+r/{r}_{s}}\right];\end{eqnarray} \tag{ 5 }$

this also equals v²(r) = G M(<r)/r as a consequence of the Birkhoff theorem.

In fractional gravity, the potential is instead derived from the modified Poisson equation (see Giusti 2020; more details can be found in Appendix A),

$\begin{eqnarray}&&{\left(-{\rm{\Delta }}\right)}^{s}\,{\rm{\Phi }}({\boldsymbol{r}})=-4\pi G\,{{\ell }}^{2-2s}\,\rho ({\boldsymbol{r}}),\end{eqnarray} \tag{ 6 }$

where ${\left(-{\rm{\Delta }}\right)}^{s}$ is the fractional Laplacian, s ∈ [1, 3/2] is the fractional index (this range of values for s is required to avoid divergences; see Appendix A), and ℓ is a fractional length scale (more on the meaning of these parameters at the end of this section). We have solved the fractional Poisson equation sourced by the NFW density distribution, as detailed in Appendix B. For s ∈ [1, 3/2), we find the solution

$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{s}(r) & = & -\displaystyle \frac{{{GM}}_{s}}{{r}_{s}}\,\displaystyle \frac{1}{{2}^{2s}\,\sqrt{\pi }}{\left(\displaystyle \frac{{\ell }}{{r}_{s}}\right)}^{2-2s}\,\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{3}{2}-s\right)}{{\rm{\Gamma }}(s+1)}\,\displaystyle \frac{{r}_{s}}{r}\\ & & \times \left\{\displaystyle \frac{2\pi s}{\sin (2\pi s)}\left[{\left(1+\displaystyle \frac{r}{{r}_{s}}\right)}^{2s-2}-{\left(1-\displaystyle \frac{r}{{r}_{s}}\right)}^{2s-2}\right]\right.\\ & & +\displaystyle \frac{{\left(r/{r}_{s}\right)}^{2s}}{1-{\left(r/{r}_{s}\right)}^{2}}\left[\left(1+\displaystyle \frac{r}{{r}_{s}}\right){}_{2}{F}_{1}\left(1,1,2s+1,\displaystyle \frac{r}{{r}_{s}}\right)\right.\\ & & \left.\left.+\left(1-\displaystyle \frac{r}{{r}_{s}}\right){}_{2}{F}_{1}\left(1,1,2s+1,-\displaystyle \frac{r}{{r}_{s}}\right)-\displaystyle \frac{4s}{2s-1}\right]\right\}\end{array}\end{eqnarray} \tag{ 7 }$

in terms of the Euler gamma function Γ and the ordinary hypergeometric function ₂ F₁; it is straightforward to verify that Φ_s=1(r) coincides with the Newtonian expression. For the limiting case s = 3/2, we instead find the solution

$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{s=3/2}(r) & = & -\displaystyle \frac{G\,{M}_{s}}{{\ell }}\,\displaystyle \frac{1}{\pi }\,\displaystyle \frac{{r}_{s}}{r}\left\{2\,\displaystyle \frac{r}{{r}_{s}}\left[\mathrm{log}\left(\displaystyle \frac{r}{{r}_{s}}\right)-1\right]\right.\\ & & -\left(1+\displaystyle \frac{r}{{r}_{s}}\right)\mathrm{log}\left(\displaystyle \frac{r}{{r}_{s}}\right)\mathrm{log}\left(1+\displaystyle \frac{r}{{r}_{s}}\right)\\ & & +\left(\displaystyle \frac{r}{{r}_{s}}-1\right){\mathrm{Li}}_{2}\left(1-\displaystyle \frac{r}{{r}_{s}}\right)\\ & & \left.-\left(1+\displaystyle \frac{r}{{r}_{s}}\right){\mathrm{Li}}_{2}\left(-\displaystyle \frac{r}{{r}_{s}}\right)+\displaystyle \frac{{\pi }^{2}}{6}\right\}\end{array}\end{eqnarray} \tag{ 8 }$

in terms of the dilogarithm function Li₂. The corresponding RC computed after ${v}_{s}^{2}(r)=r\,| d{{\rm{\Phi }}}_{s}/{dr}|$, which remains true because the standard law of inertia continues to hold in fractional gravity, is illustrated in Figure 1 (top panel) for different values of the fractional index s. For s increasing above 1 (corresponding to Newtonian gravity), we find that the RC is steeper in the inner region, behaving as v_s (r) ∝ r^s−1/2 for s < 3/2 and ${v}_{s=3/2}(r)\propto r\,\sqrt{-\mathrm{log}r/{r}_{s}}$ for s = 3/2; curiously, for the limiting value s = 3/2, the RC tends to closely mirror a (Newtonian) solid body. Contrariwise, in the outskirts, the RC in fractional gravity is appreciably flatter than the Newtonian case, behaving as v_s (r) ∝ r^s−3/2 for s < 3/2 and ${v}_{s=3/2}(r)\propto \sqrt{\mathrm{log}r/{r}_{s}}$ for s = 3/2.

Zoom In Zoom Out Reset image size

In fractional gravity, the Birkhoff theorem does not hold, but one can insist in writing ${v}_{s}^{2}(r)=G\,{M}_{\mathrm{eff},s}(\lt r)/r$ in terms of an effective mass M_eff,s (<r) and then differentiate this expression to obtain an effective density profile via ρ_eff,s (r) = (1/4π r²) × dM_eff,s /dr. This is actually the density behavior that one would infer by looking at the RC and interpreting the result in terms of Newtonian gravity. We illustrate such an effective density profile in Figure 1 (bottom panel) for different values of the fractional index s. With s increasing from unity (Newtonian case), the effective density profile progressively flattens. Specifically, in the inner region, a core-like behavior tends to be enforced, while in the outskirts, the effective profile resembles an isothermal sphere. To have a grasp on the overall effect, consider the s = 3/2 case, where the effective density has a particularly simple analytic expression, given by

$\begin{eqnarray}&&{\rho }_{\mathrm{eff},s=3/2}(r)=\displaystyle \frac{2\,{\rho }_{s}\,{r}_{s}}{\pi \,{\ell }}\,\displaystyle \frac{\mathrm{ln}(r/{r}_{s})}{{\left(r/{r}_{s}\right)}^{2}-1};\end{eqnarray} \tag{ 9 }$

interestingly, the logarithmic slope of this profile is −2 in the outskirts, −1 at r ≈ r_s , and ≲−0.5 for r ≲ r_s /10. All in all, the effective density is practically indistinguishable from the run of a cored pseudoisothermal sphere. This slope around r ∼ 0.1r_s is remarkably close to the observational value estimated from the kinematic of individual dwarf galaxies (e.g., Oh et al. 2015; Salucci et al. 2021).

The physical reason behind this behavior is that the fractional Laplacian is inherently nonlocal; the effect is markedly evident in Fourier space, where it is easily found that ${ \mathcal F }\{{\rho }_{\mathrm{eff},s}\}/{ \mathcal F }\{\rho \}\propto {(\lambda /{\ell })}^{2s-2}$ in terms of the Fourier transform ${ \mathcal F }$ and the wavelength λ = 2π/k associated with the Fourier mode with wavenumber k. This means that fractional gravity causes, with respect to the Newtonian case, a transfer of power from modes with λ ≲ ℓ, which are suppressed, to modes with λ ≳ ℓ, which are instead enhanced. Therefore, the parameter s measures the strength of the nonlocality, while the length scale ℓ can be interpreted as the typical size below which gravitational effects are somewhat reduced and above which they are instead amplified by nonlocality (around r ≈ ℓ, the dynamics is almost unaffected and indistinguishable from the Newtonian case).

2.2. Bayesian Analysis of Stacked RCs

The outcomes of the fitting procedure are illustrated in Figure 2 for the dwarf sample, Figure 3 for HSB spirals, and Figure 4 for LSB spirals. In each panel, the best fit (solid lines) and the 1σ credible intervals sampled from the posterior (shaded areas) are shown for the halo (green), disk (blue), and total (red) RC models; for reference, the results in Newtonian gravity are also displayed (dashed lines). The MCMC posterior distributions in the fractional gravity framework for some representative bins are shown in the corner plot of Figure 5. Finally, Table 1 reports the marginalized posterior estimates of the fitting parameters, the reduced ${\chi }_{r}^{2}$ of the fits, and the difference in the Bayesian inference criterion (BIC) for the model comparison between the Newtonian and fractional gravity fits. ⁶

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Marginalized Posterior Estimates (Mean and 1σ Confidence Intervals) for the Fractional and Newtonian Gravity

	Fractional Gravity					Newtonian Gravity
Bin	s	$\mathrm{log}{\ell }/{r}_{s}$	$\mathrm{log}{M}_{{\rm{H}}}$ (M⊙)	$\mathrm{log}{M}_{d}$ (M⊙)	${\chi }_{r}^{2}$	$\mathrm{log}{M}_{{\rm{H}}}$ (M⊙)	$\mathrm{log}{M}_{d}$ (M⊙)	${\chi }_{r}^{2}$	ΔBIC
Dwarfs	${1.41}_{-0.04}^{+0.05}$	$-{0.72}_{-0.04}^{+0.03}$	${9.18}_{-0.26}^{+0.22}$	${6.59}_{-0.38}^{+0.38}$	4.91	${10.52}_{-0.01}^{+0.03}$	${6.25}_{-0.57}^{+0.57}$	25.5	−250
HSB1	${1.38}_{-0.06}^{+0.07}$	$-{1.07}_{-0.09}^{+0.11}$	${9.68}_{-0.36}^{+0.31}$	${7.38}_{-0.89}^{+0.79}$	1.37	${11.69}_{-0.03}^{+0.03}$	${6.78}_{-0.51}^{+0.51}$	4.36	−50
HSB2	${1.35}_{-0.06}^{+0.06}$	$-{1.39}_{-0.17}^{+0.17}$	${9.85}_{-0.36}^{+0.40}$	${8.03}_{-0.67}^{+1.3}$	0.94	${12.34}_{-0.06}^{+0.06}$	${8.81}_{-0.07}^{+0.20}$	3.54	−42
HSB3	${1.30}_{-0.04}^{+0.03}$	$-{1.11}_{-0.13}^{+0.21}$	${10.70}_{-0.16}^{+0.11}$	${9.71}_{-0.03}^{+0.10}$	1.98	${12.24}_{-0.06}^{+0.04}$	${9.62}_{-0.01}^{+0.06}$	4.87	−51
HSB4	${1.27}_{-0.05}^{+0.05}$	$-{1.17}_{-0.17}^{+0.27}$	${10.98}_{-0.18}^{+0.14}$	${9.87}_{-0.11}^{+0.21}$	1.70	${12.40}_{-0.08}^{+0.03}$	${9.88}_{-0.06}^{+0.09}$	3.53	−30
HSB5	${1.03}_{-0.03}^{+0.10}$	$-{0.72}_{-0.89}^{+0.89}$	${12.16}_{-0.20}^{+0.31}$	${9.99}_{+0.18}^{+0.22}$	17.7	${12.41}_{-0.06}^{+0.04}$	${10.17}_{-0.01}^{+0.05}$	15.8	+5
HSB6	${1.11}_{-0.04}^{+0.05}$	$-{1.59}_{-0.37}^{+0.37}$	${11.74}_{-0.54}^{+0.27}$	${10.34}_{-0.09}^{+0.10}$	0.79	${12.60}_{-0.06}^{+0.04}$	${10.27}_{-0.03}^{+0.06}$	1.13	−2
HSB7	${1.08}_{-0.07}^{+0.07}$	$-{0.68}_{-0.78}^{+0.78}$	${12.28}_{-0.31}^{+0.37}$	${10.61}_{+0.04}^{+0.03}$	0.52	${12.45}_{-0.07}^{+0.05}$	${10.54}_{-0.02}^{+0.05}$	0.47	+5
HSB8	${1.02}_{-0.05}^{+0.05}$	$-{0.64}_{-0.86}^{+0.86}$	${12.16}_{-0.07}^{+0.15}$	${10.82}_{+0.02}^{+0.05}$	3.74	${12.24}_{-0.09}^{+0.06}$	${10.79}_{-0.04}^{+0.06}$	3.33	+5
HSB9	${1.09}_{-0.09}^{+0.09}$	$-{0.40}_{-0.76}^{+0.76}$	${12.41}_{-0.37}^{+0.27}$	${10.97}_{-0.04}^{+0.06}$	1.42	${12.48}_{-0.09}^{+0.05}$	${10.88}_{-0.06}^{+0.08}$	1.26	+5
HSB10	${1.11}_{-0.09}^{+0.09}$	$-{0.29}_{-0.78}^{+0.78}$	${12.45}_{-0.46}^{+0.03}$	${11.19}_{-0.04}^{+0.06}$	2.25	${12.26}_{-0.08}^{+0.05}$	${11.11}_{-0.05}^{+0.07}$	1.99	+5
HSB11	${1.12}_{-0.09}^{+0.09}$	$-{0.02}_{-0.65}^{+0.65}$	${12.56}_{-0.34}^{+0.10}$	${11.43}_{-0.03}^{+0.04}$	0.66	${12.39}_{-0.09}^{+0.07}$	${11.35}_{-0.03}^{+0.06}$	0.59	+6
LSB1	${1.36}_{-0.06}^{+0.11}$	$-{0.31}_{-0.19}^{+0.07}$	${9.93}_{-0.55}^{+0.42}$	${7.57}_{-0.6}^{+1.1}$	1.84	${10.46}_{-0.01}^{+0.04}$	${6.67}_{-0.66}^{+0.66}$	12.3	−103
LSB2	${1.08}_{-0.09}^{+0.09}$	$-{1.14}_{-0.79}^{+0.79}$	${10.83}_{-0.03}^{+0.67}$	${8.17}_{-0.15}^{+0.90}$	0.65	${11.29}_{-0.02}^{+0.06}$	${8.39}_{-0.06}^{+0.47}$	0.52	+5
LSB3	${1.23}_{-0.10}^{+0.06}$	$-{0.17}_{-0.49}^{+0.54}$	${11.44}_{-0.69}^{+0.60}$	${9.81}_{-0.03}^{+0.27}$	0.25	${11.64}_{-0.04}^{+0.04}$	${9.29}_{-0.05}^{+0.17}$	2.37	−16
LSB4	${1.03}_{-0.06}^{+0.06}$	$-{0.80}_{-0.84}^{+0.84}$	${11.39}_{-0.04}^{+0.18}$	${10.35}_{+0.07}^{+0.13}$	2.42	${11.48}_{-0.03}^{+0.12}$	${10.43}_{-0.02}^{+0.04}$	1.74	+4
LSB5	${1.05}_{-0.09}^{+0.09}$	$-{0.70}_{-0.91}^{+0.91}$	${11.64}_{-0.28}^{+0.36}$	${10.73}_{-0.36}^{+0.39}$	8.9	${11.99}_{-0.06}^{+0.02}$	${11.04}_{-0.05}^{+0.08}$	8.8	−9

Note. The columns report the values of the fractional parameter s, fractional scale length ℓ, halo mass M_H, disk mass M_d , and reduced ${\chi }_{r}^{2}$. The last column refers to the difference in values of the BIC between the two models. Note that for values s ≲ 1.1, the estimate of ℓ is highly uncertain, since the fractional model features a rotation velocity approximately independent of such parameters and tends to become very similar to Newtonian gravity.

Download table as:  ASCIITypeset image

The fractional gravity framework performs substantially better than the Newtonian case in several bins, especially those relative to small and intermediate-mass systems, where the DM component dominates or appreciably contributes to the total RC out to r_opt. This is manifest both from the quality of the fit in terms of a substantial improvement in the reduced ${\chi }_{r}^{2}$ and in terms of the strongly negative values of the ΔBIC (meaning that fractional gravity is favored, in a Bayesian sense, over the Newtonian case). For example, one can consider the dwarf RC shown in Figure 2, where the contribution of baryons to the dynamics is known to be negligible; hence, the overall RC is determined solely by the DM component. In the Newtonian case, the cuspy behavior of the NFW density distribution (with standard values of the concentration) struggles in reproducing the inner part of the RC; historically, this has been among the first evidence of the cusp–core problem. However, the dynamics in fractional gravity is substantially modified in such a way that the RC sourced by the NFW density profile can actually reproduce the observed inner RC to a remarkable accuracy; this is because, as discussed in Section 2, the gravitational potential sourced by the NFW density in fractional gravity mirrors that sourced by a cored density distribution (e.g., a pseudoisothermal sphere) in Newtonian gravity. In more massive galaxies, where the contribution of baryons to the inner RC out to r_opt becomes increasingly relevant and eventually dominates over the DM, the evidence for fractional gravity in the DM gets weaker; in such systems, the fits in fractional and Newtonian gravity turn out to be very similar (though with an appreciable difference in the estimated total mass; see below). However, we offer the caveat that the contribution of the DM to the RC in massive spirals emerges at radii r ≳ r_opt, where the data are still subject to appreciable uncertainties; more precise observations over an extended radial range would be necessary before drawing a definitive conclusion on the impact of fractional gravity in massive spirals.

In Figure 6, we illustrate the disk (stellar) versus halo mass relation constructed with the best-fit estimates in fractional (magenta circles) and Newtonian (orange squares) gravity. For reference, the dotted–dashed and dashed lines are the relations inferred by Lapi et al. (2018) and Karukes & Salucci (2017) by fitting galaxy RCs of different masses with the cored Burkert profile that are also in agreement with abundance-matching determinations (Moster et al. 2013; Lapi et al. 2018). The inferred masses in the fractional gravity framework are in excellent agreement with the Lapi et al. (2018) and Karukes & Salucci (2017) results. Contrariwise, there is instead a clear tendency for the NFW fits in Newtonian gravity to yield higher DM masses at a fixed stellar mass (especially so for small galaxies) that is somewhat inconsistent with the aforementioned empirical determinations; in the literature, it has been pointed out that more reasonable halo masses can be obtained by leaving the concentration parameter free to vary in the fit, but at the price of obtaining concentration values in disagreement with the concentration versus mass relation measured in cosmological simulations (e.g., Dehghani et al. 2020; Gandolfi et al. 2022a). We fit the stellar mass versus halo mass relation in fractional gravity via an orthogonal distance regression (ODR) algorithm that takes into account error bars in both axes. Using a linear (solid line) shape $\mathrm{log}{M}_{d}[{M}_{\odot }]=a+b(\mathrm{log}{M}_{{\rm{H}}}[{M}_{\odot }]-11)$, we obtain the best-fit parameters a = 9.31 ± 0.11 and b = 1.34 ± 0.09 and a reduced ${\chi }_{r}^{2}\approx 0.37;$ a nonlinear fit (dotted line) inspired by Equation (2) in Moster et al. (2013) yields only a modest improvement in the reduced ${\chi }_{r}^{2}\approx 0.31$ and large uncertainties on the fitted parameters.

Zoom In Zoom Out Reset image size

In Figure 7, we illustrate the behavior of the fractional gravity index s and length scale ℓ on halo mass. The former features a decreasing behavior as a function of halo mass, passing from values around s ≈ 1.3–1.4 in dwarf galaxies with M_H ≲ 10¹⁰ M_⊙, to s ≈ 1.2–1.3 in intermediate-mass galaxies, to s ≲ 1.1 in massive galaxies with M_H ≳ some 10¹¹ M_⊙. This implies a strongly fractional behavior of the DM component in small and intermediate-mass systems, whereas deviations from Newtonian gravity in high-mass systems seem to be more limited (at least given the present quality of the data in the outer regions, where the DM component can be probed in massive galaxies; see above). We fit the s–M_H relation with a linear shape $s=a+b(\mathrm{log}{M}_{{\rm{H}}}[{M}_{\odot }]-11)$ via an ODR algorithm to obtain the best-fit parameters a = 1.21 ± 0.02 and b = −0.13 ± 0.01 and a reduced ${\chi }_{r}^{2}\approx 0.58;$ a nonlinear fit (dotted line) $s=\displaystyle \frac{5}{4}+\tfrac{1}{4}\,\tanh (c(\mathrm{log}{M}_{{\rm{H}}}[{M}_{\odot }]-d)]$ with asymptotic values of s = 1 and 1.5 at small and large masses yields the best-fit parameters c = −0.68 ± 0.09 and d = 10.71 ± 0.12 and a modest improvement in the reduced ${\chi }_{r}^{2}\approx 0.52$.

Zoom In Zoom Out Reset image size

As to the fractional length scale, there is a tendency of increasing with halo mass. However, it may be stressed that the determination of ℓ is robust only in small and intermediate-mass galaxies with strong fractional gravity behavior, or, in other words, where the index s is substantially larger than unity. Contrariwise, in massive galaxies, where s gets close to 1, the fractional Poisson equation (Equation (6)) becomes almost independent on ℓ (the quantity ℓ^2s-2 appears on the right-hand side), so that any robust inference on this parameter is difficult. We fit the ℓ−M_H relation with a linear shape ${\ell }[\mathrm{kpc}]=a+b(\mathrm{log}{M}_{{\rm{H}}}[{M}_{\odot }]-11)$ via an ODR algorithm to obtain the best-fit parameters a = 0.27 ± 0.13 and b = 0.57 ± 0.11 and a reduced ${\chi }_{r}^{2}\approx 1.25;$ a nonlinear fit (dotted line) in terms of a double power-law shape leaves the reduced ${\chi }_{r}^{2}\approx 1.32$ unchanged and yields large errors on the fitted parameters. Interestingly, the best-fit linear relation implies a scaling of ℓ very similar to the NFW scale radius r_s , in such a way that ℓ/r_s ≈ 0.25 ± 0.15 is approximately independent of halo mass.

The trends of the parameters s and ℓ with halo mass may be suggestive that the effects of fractional gravity on progressively larger scales become weaker; this would retain the successes of the standard DM paradigm in Newtonian gravity on cosmological scales (see Section 1). However, given the still large uncertainties affecting the kinematic analysis of massive galaxies, it would be interesting to validate our results by looking at larger systems like clusters of galaxies. In fact, the overall gravitational potential of the cluster is dominated by the DM component, while most of the baryons are in the form of a hot intracluster medium in hydrostatic equilibrium within the DM potential well. In fractional gravity, one may expect the thermodynamic properties of the intracluster medium (e.g., density, temperature, pressure, entropy runs) to be substantially modified, so comparing with those empirically inferred from X-ray, Sunyaev–Zel'dovich effect, and weak lensing data could be extremely informative on fractional gravity parameters at scales much larger than galaxies. We plan to follow this program in a forthcoming paper.

3.1. Radial Acceleration Relation

The RAR, a relationship between the total g_tot versus the baryonic acceleration g_bar, has been proposed by McGaugh et al. (2016) via the analysis of individual RCs of the SPARC sample (Lelli et al. 2016); note that the acceleration of any component can be simply computed from the RC as g ≡ v²(r)/r. The RAR is thought to subsume and generalize a plethora of well-known dynamical laws of galaxies (see discussion by Lelli et al. 2017), although its tightness and physical interpretation are still debated (e.g., Di Paolo et al. 2019; Salucci et al. 2020). Here we just aim to show the RAR that emerges from the dynamics of galaxies in the fractional gravity framework, using as inputs the outcomes from the previous analysis of stacked RC data.

Toward this purpose, note that the RAR is a local scaling law that combines data at different radii in galaxies with different masses, which feature different contributions of the stellar/gas disk and DM. Therefore, we approach the problem via a semiempirical method. We first build up mock RCs of galaxies with different halo masses by exploiting well-known empirical relationships between halo mass and baryonic properties; then we sample the RCs to derive the total and baryonic accelerations and construct the RAR. Such a procedure was described in Gandolfi et al. (2022b; see their Section 4), to which we refer the reader for details. The only difference with respect to that work is that the RC of the DM component is modeled in fractional gravity (with an NFW density profile).

The RAR is illustrated in Figure 8. The observational determinations are shown in gray; data for spiral galaxies (binned) are from McGaugh et al. (2016; squares), Local Group dwarf spheroidals are from Lelli et al. (2017; triangles), Coma Cluster ultradiffuse galaxies are from Freundlich et al. (2022; reversed triangles), Centaurus A dwarf spheroidals are from Muller et al. (2021; diamonds), and dwarf LSBs with g_bar(0.4 ≲ r/R_opt ≲ 1) < −11 are from Di Paolo et al. (2019; stars). For reference, the dotted black line displays the one-to-one relation g_tot = g_bar.

Zoom In Zoom Out Reset image size

The big green circles illustrate the outcome for s = 1 (independent of ℓ and corresponding to Newtonian gravity); as pointed out by Gandolfi et al. (2022a), the observed RAR is reasonably reproduced at large baryonic acceleration g_bar, but it is considerably overpredicted at small g_bar. The big red and blue circles show the outcomes in fractional gravity for s = 1.2 and 1.4, respectively, and ℓ/r_s = 0.25 (typical values from the stacked RC analysis); the shaded area displays instead how the predicted RAR changes when ℓ/r_s is varied in the range 0.1–0.5 (extremal values from the stacked RC analysis). Finally, the solid black line is the outcome by assuming the scaling of s and ℓ/r_s with halo mass M_H from the stacked RC analysis (see Figure 7). All in all, the RAR in fractional gravity has a shape and an overall scatter that are pleasingly consistent with the data.

These results can be understood as follows. The high-acceleration regime of the RAR is mainly dominated by the contribution at small radii in high-mass galaxies, where the total gravitational acceleration is dominated by baryons, to imply g_tot ≈ g_bar. Contrariwise, at small g_bar, the behavior of the RAR is dictated by small/intermediate radii in intermediate- and low-mass galaxies, where the baryon acceleration is dominated by the stellar disk, while the total gravitational acceleration is contributed by both the disk and the halo; for the NFW density profile in Newtonian gravity, the halo typically dominates over the disk, and the RAR deviates upward with respect to the high-acceleration regime. However, in fractional gravity, such an upward deviation is mitigated because the inner RC corresponding to the NFW density cusp is somewhat steepened, ${v}_{\mathrm{DM}}\propto {r}^{s-1/2}\,{({\ell }/{r}_{s})}^{1-s}$; hence, the total acceleration toward the center is appreciably smaller than in Newtonian dynamics, especially for s approaching the limiting value 1.5 and for higher ℓ/r_s ratios. We conclude by noticing that this behavior is different from intrinsically cored models in Newtonian dynamics that, as demonstrated by Gandolfi et al. (2022b; see their Figure 7), tend to enforce a low-acceleration behavior similar to (or only slightly different from) the high-acceleration regime.

3.2. Universal Core Surface Density and Core Radius versus Disk Scale Length

It has been well established (see Salucci & Burkert 2000; Donato et al. 2009; Gentile et al. 2009; Burkert 2015) that, when fitting galaxy RCs with a cored profile like Burkert or a pseudoisothermal sphere, the product of the core density ρ₀ and core radius r₀, which is a kind of core surface density, is an approximately universal constant among different galaxies; typical values of the latter have been estimated around ${\rho }_{0}\,{r}_{0}\approx {140}_{-30}^{+80}\,{M}_{\odot }$ pc⁻² by Donato et al. (2009) and around ${\rho }_{0}\,{r}_{0}\approx {75}_{-45}^{+85}\,{M}_{\odot }$ pc⁻² by Burkert (2015). This universality in the core surface density is a somewhat unexpected property that poses a serious challenge to any theoretical model of core formation (e.g., Deng et al. 2018; Burkert 2020). We now aim to illustrate that halos in the fractional gravity framework are broadly consistent with such a remarkable scaling law.

To fairly compare with the literature universal core surface density relation, we need to identify an equivalent core density and radius for a DM halo described by fractional gravity. To this purpose, we start from the best-fit rendition of the stacked RCs of Section 3 in the fractional gravity framework, refit each of these RCs with an equivalent Burkert or pseudoisothermal sphere RC via a least-squares minimization algorithm, and finally extract the core density and radius. The outcome of this procedure is illustrated in Figure 9 (top panel), where the core surface density is plotted against the core radius; magenta circles refer to the pseudoisothermal sphere refitting, and cyan squares refer to the Burkert profile refitting. All in all, the measured core surface densities in the fractional gravity framework are broadly consistent with the universal relation and its overall normalization.

Zoom In Zoom Out Reset image size

Another interesting relation known to hold (see Donato et al. 2004; Salucci et al. 2007; Karukes & Salucci 2017; Lapi et al. 2018; Salucci et al. 2020) from dwarf to massive galaxies involves the core radius r₀ and disk scale length r_d . Its origin is uncertain, with some of the above authors arguing for evidence in favor of a nongravitational interaction between baryons and DM. To check such a relationship in fractional gravity, we extract r₀ as described above for every stacked RC bin and plot it as a function of the average r_d = r_opt/3.2 in the bin. The outcome is illustrated in Figure 9 (bottom panel). We find a pleasing agreement with the determinations by Donato et al. (2004) and Salucci et al. (2020). In fractional gravity, the r₀ versus r_d relation can be traced back to nonlocal effects that ultimately enforce the scaling with halo mass of the core radius and the disk mass, the latter being, in turn, trivially related to the disk scale length.

As recalled in Section 1, the evidence for the presence of DM in the universe goes much beyond galaxy kinematics; thus, our fiducial approach in this work has been to explore whether fractional gravity associated with the DM component could possibly alleviate some of the small-scale issues emerging from the analysis of galaxy RCs, especially dwarfs. However, for the sake of completeness and to ease the comparison with other literature studies, it is worth exploring an alternative MONDian viewpoint, where no DM component is considered and fractional gravity is associated with the baryons. The underlying idea is that, similar to the original version of the MOND or other modified gravity theories, DM phenomenology could be emulated by baryons in fractional gravity.

We have repeated the analysis of Section 2.2, adopting this MONDian viewpoint. Specifically, we have assumed that no DM component is present, and we have fitted the stacked RCs using the circular velocity of a razor-thin exponential disk model in fractional gravity (provided in Appendix B; we recall that for the systems considered in this work, the contribution to the stacked RCs from any bulge or gas component can be safely neglected). In this setup, the free parameters of the model are the disk mass M_d , fractional index s, and fractional length scale ℓ. As for ℓ, we have tried to either let it vary as a free fitting parameter or set it to the value ${\ell }\sim (2/\pi )\times \sqrt{G\,{M}_{d}/{a}_{0}}$, where a₀ ≈ 1.2 × 10⁻¹⁰ m s⁻², that can be naively expected in a MONDian interpretation of the fractional gravity theory (see the end of Appendix A); we find that the fits to the stacked RCs in these two cases are practically indistinguishable, so we focus on the latter one. The results of the fitting procedure on a subset of representative RCs of our sample are summarized in Figure 10. In each panel, we illustrate the best-fit RC in the MONDian approach with no DM and report for comparison the one from our fiducial setup with DM; we also indicate the corresponding best-fit estimates of the fractional index s and disk mass M_d .

Zoom In Zoom Out Reset image size

An important point to stress is that the MONDian viewpoint shares with our fiducial setup a similar decreasing trend of s, from values of s ≈ 1.5 in small galaxies toward s ≈ 1 in the largest systems. This is essentially imposed by the behavior of the measured RCs in the inner region, which behave like a solid body in dwarfs while steeply rising in massive galaxies. However, it is seen that the MONDian viewpoint struggles somewhat in reproducing the outer behavior of the RCs, especially in small galaxies. This is because for an s appreciably larger than 1, the RC of a razor-thin disk in fractional gravity (and this is true for other baryonic profiles as well; see Appendix B) tends to saturate to a constant in the outer regions, while the measured RCs in small galaxies are quite steeply rising even in the outermost sampled point. Moreover, the fit of the MONDian viewpoint in the outer regions does not substantially improve in moving toward massive galaxies; this is because s progressively decreases (as required to fit the inner part; see above), so the outer RC in fractional gravity tends to recover the steeply decreasing trend of the Newtonian case, while the measured RCs retain a rather flat or modestly decreasing behavior.

One can speculate that perhaps the MONDian viewpoint can be reconciled with the observed RCs if the fractional index s is allowed to vary with the radius to originate a sort of variable-order fractional theory (see Giusti 2020). However, some caveats are in order. First, nontrivial numerical techniques for solving the fractional Poisson equation must be developed in such a case, since the analytic treatment becomes impossible. Second, it is not clear whether the variable-order theory is conceptually well posed, since fractional gravity is inherently nonlocal; hence, a simple dependence s(r/ℓ) on the local radial coordinate does not seem easily justified. Third, to describe the radial dependence of s, other parameters must be introduced that plainly make the theory subject to additional degeneracies and less elegant or predictive.

Another issue with the MONDian viewpoint is that, lacking the DM component, the correct normalization of the RC (even in the inner region) can be obtained only with disk masses systematically larger than in our fiducial setup. This may constitute a problem, especially in dwarf spheroidals, since the visible mass that can be inferred from photometry (I band) could be appreciably smaller. In other words, the mass-to-light ratios that can be derived from the best-fit stellar mass in the MONDian viewpoint tend to be substantially larger than the typical values M_⋆/L_I ≲ 1.5 expected on the basis of stellar population synthesis models (e.g., Portinari et al. 2004). To highlight this point, in each panel of Figure 10, we have reported the average I-band luminosity in the bin and the maximum disk mass that can be inferred from such a photometry according to the relationship between luminosity and disk-to-total mass at the optical radius estimated by Salucci et al. (2008; see also Persic & Salucci 1990). There is clearly a tendency for the best-fit stellar masses in the MONDian viewpoint to exceed such values.

We also note that the values of the fractional length scale derived in the MONDian viewpoint from the disk masses according to the relation ${\ell }\sim (2/\pi )\times \sqrt{G\,{M}_{d}/{a}_{0}}$ range from 0.3 kpc in dwarfs to about 10 kpc in massive galaxies. For comparison, in our fiducial approach with DM, ℓ ranges from 0.3 kpc in dwarfs to a few tens of kiloparsecs in massive galaxies (see right panel in Figure 7). Thus, the overall dependence of ℓ on the total mass in the two approaches is rather similar within the respective fitting uncertainties.

In the context of the MONDian viewpoint, it is worth mentioning literature studies where fractional gravity is invoked as an alternative to DM in galaxies. A relevant example is the theory of Newtonian fractional-dimensional gravity by Varieschi (2020, 2021, 2022, 2023), which introduces a generalized law of Newtonian gravity in a spatial dimension D < 3, representing the local effective Hausdorff dimension of the matter distribution. Another approach by Calcagni (2010, 2013, 2021; see also Calcagni & Varieschi 2022) relies on multifractional spacetimes with variable Hausdorff and spectral dimensions directly inspired from quantum gravity theories. At variance with our fiducial setup, both of these theories adopt a MONDian viewpoint where DM is not present, and the galaxy kinematics is interpreted as a pure geometrical effect. Another important difference with our approach (even in the MONDian viewpoint with no DM) is that the above theories are formulated in terms of integer-order local operators (though generalizations are possible), whereas the fractional Laplacian exploited in this study is inherently nonlocal. Recently, both Varieschi's and Calcagni's theories have been tested against a few prototypical examples of individual galaxy RCs. The former decently reproduces the measured RCs, though at the price of introducing a radially dependent Hausdorff dimension. The latter can easily accommodate the observed rising or flattening behavior of the observed RCs in the intermediate and outer regions, while it struggles somewhat in fitting the inner behavior. All in all, such theories could be promising, but an extended analysis of different mass systems, and especially dwarf galaxies with poor baryonic content, would be interesting.

In this work, we have explored the possibility that the dark matter (DM) component in galaxies may originate fractional gravity; in such a framework, albeit the standard law of inertia continues to hold, the gravitational potential associated with a given DM density distribution is determined by a modified Poisson equation including fractional derivatives (i.e., derivatives of noninteger type) that accounts for nonlocal effects.

We have shown that in fractional gravity, the dynamics of a test particle in the potential sourced by an NFW density distribution is substantially altered with respect to the Newtonian case (i.e., based on the standard Poisson equation), mirroring what in Newtonian gravity would be originated by cored density profiles (like a pseudoisothermal sphere). The fractional gravity framework can be fully characterized by two parameters: the index s of the fractional Poisson equation ruling the strength of the nonlocal effects and a fractional length scale ℓ marking the size below which gravitational effects are somewhat reduced and above which they are amplified by nonlocality.

We have then tested such a model against the stacked RCs of local spiral galaxies with different properties (e.g., optical radius/velocity) and morphology (HSB and LSB spirals, dwarfs). Our main results can be summarized as follows.

• 1.  
  We have shown that the fractional gravity framework performs substantially better than the Newtonian case in fitting stacked RCs, especially for dwarfs and intermediate-mass galaxies where the DM component dominates or appreciably contributes to the kinematics in the inner regions.
• 2.  
  Interestingly, the strength of fractional gravity effects gets progressively weaker in more massive systems dominated by baryons toward the center, though precise measurements of the outer RCs in massive spirals are needed to robustly confirm this finding. If true, this trend will imply that our fractional gravity framework can substantially alleviate the small-scale issues of the standard DM paradigm in Newtonian gravity while saving its successes on large cosmological scales.
• 3.  
  We have derived the relationship between (stellar) disk and halo mass in the fractional gravity framework, finding remarkable agreement with literature determinations based on empirical RC modeling via cored profiles in Newtonian gravity.
• 4.  
  We have computed the radial acceleration relation (RAR) in the fractional gravity framework, finding that its behavior at low baryonic acceleration appreciably deviates from that in Newtonian gravity to yield an overall shape and scatter in pleasing agreement with the current observational determinations.
• 5.  
  We have highlighted that the fractional gravity originates RCs that are consistent with the observed universal core surface density behavior and the observed scaling of the core radius versus the disk scale length. Such relationships, which pose a serious challenge to other models of modified gravity and/or nonstandard particle candidates, are naturally originated by the nonlocal nature inherent to the fractional gravity framework.
• 6.  
  We have explored an alternative MONDian viewpoint on fractional gravity where the DM component is not present, and we have discussed its performance and criticalities in reproducing the observed galaxy kinematics.

All in all, the fractional gravity framework effectively describes nonlocal effects in the DM component as a source of gravity, but how are these originated? Although the investigation of such an issue from the theoretical side is far beyond the scope of the present paper, in the vein of a qualitative discussion, we can envisage three possibilities: (i) nonlocality may stem from microscopic properties of the DM particles, e.g., some form of quantum entanglement, nonstandard interactions, or coupling with the gravitational metric; (ii) nonlocal dynamics may be enforced at the mesoscopic level, related to the fluid, coarse-grained description of the DM particles' collective behavior in a finite volume; and (iii) nonlocality may be sourced on macroscopic scales by the response of a complex (e.g., clumpy and inhomogeneous) DM distribution to the long-range action of gravitational forces against DM and baryonic particles. That said, it should be noted that fractional gravity is not necessarily meant to be an ab initio theory but may constitute an effective description for a whole class of models that imply the development of nonlocal effects in the DM component when sourcing the gravitational field.

In a future perspective, it will be interesting to investigate the performance of the fractional gravity framework in fitting the outermost RCs of z ∼ 1 spiral galaxies (see Sharma et al. 2022) to more robustly check whether the strength and length scale of the fractional dynamics show signs of evolution with redshift. In addition, our analysis of galaxy kinematics suggests that fractional gravity effects tend to fade away when moving to progressively massive systems; to confirm such a trend, a direct inspection of fractional gravity behavior on larger cosmological scales would be welcome. A natural test bed would be provided by galaxy clusters (e.g., Gandolfi et al. 2023); specifically, we will investigate to what extent the thermodynamic profiles of the intracluster medium are affected by the modified DM gravitational potential. We also plan to revise the interpretation of strong and weak gravitational lensing events that may be substantially altered in the fractional gravity framework. Finally, a more theoretical development will involve applying the fractional gravity framework to nonstandard DM particle candidates and checking whether nonlocal effects may solve the issues affecting these scenarios; e.g., halos in fuzzy DM have an inner core, but the universal core surface density relation cannot be reproduced in Newtonian gravity (see Burkert 2020). Ultimately, it would be relevant to trace the physical origin of fractional gravity, especially in terms of the mechanism determining the index s and length scale ℓ in different systems and in originating the scaling of these quantities between themselves and with halo mass (so actually reducing the number of effective parameters in the theory). Finally, full N-body simulations could be exploited to study the emergence of the fractional gravitational dynamics in evolving cosmological structures.

We acknowledge M. Giulietti and M. Massardi for helpful discussions. We thank the anonymous referee for constructive comments. A.L. has been supported by the EU H2020-MSCA-ITN-2019 Project 860744 "BiD4BESt: Big Data applications for black hole Evolution STudies" and PRIN MIUR 2017 prot. 20173ML3WW, "Opening the ALMA window on the cosmic evolution of gas, stars and supermassive black holes."

Fractional calculus is the field of mathematics dealing with differentiation and integrations of noninteger order. Nowadays, it has found many applications in material science, rheology, seismology, transport phenomena, nuclear physics, medicine, finance, etc. Given that the concepts and techniques related to fractional calculus are not very common among the astrophysics community, for the reader's convenience, we present a short (nonexhaustive) primer in this Appendix. In particular, we focus on the definitions of the fractional integrals and derivatives, fractional Laplacian, and fractional Poisson equation that are often invoked and employed in the main text.

A.1. Fractional Integrals and Derivatives

In one dimension, the simplest definition of a fractional integral can be viewed as a natural extension of the well-known formula (usually attributed to Cauchy)

$\begin{eqnarray}&&{{ \mathcal I }}^{n}f(x)=\displaystyle \frac{1}{(n-1)!}\,{\int }_{0}^{x}{dy}\,f(y){\left(x-y\right)}^{n-1},\end{eqnarray} \tag{ A1 }$

where n is an integer, and n ! = 1 × 2 × ... × n is the factorial function. The above relates the n-fold primitive of a function f(x) to a convolution integral, as it can be easily verified by repeated integration by parts, provided that f satisfies f(x) = 0 for x < 0. One is immediately led to extend the above formula for any positive real index s by means of the Euler gamma function ${\rm{\Gamma }}(s)\equiv {\int }_{0}^{\infty }{dx}\,{x}^{s-1}\,{e}^{-x}$. In fact, since Γ(n) = (n − 1)! for an integer n, one can define the fractional integral of order s via the expression

$\begin{eqnarray}&&{{ \mathcal I }}^{s}f(x)=\displaystyle \frac{1}{{\rm{\Gamma }}(s)}\,{\int }_{0}^{x}{dy}\,f(y){\left(x-y\right)}^{s-1},\end{eqnarray} \tag{ A2 }$

prescribing ${{ \mathcal I }}^{0}f(x)=f(x)$. The fractional integral enjoys the (semigroup) property ${{ \mathcal I }}^{s}\,{{ \mathcal I }}^{p}={{ \mathcal I }}^{s+p}$, which implies commutativity ${{ \mathcal I }}^{s}\,{{ \mathcal I }}^{p}={{ \mathcal I }}^{p}\,{{ \mathcal I }}^{s}$ and the straightforward action on power function

$\begin{eqnarray}&&{{ \mathcal I }}^{s}{x}^{a}=\displaystyle \frac{{\rm{\Gamma }}(a+1)}{{\rm{\Gamma }}(a+s+1)}\,{x}^{a+s}.\end{eqnarray} \tag{ A3 }$

It is natural to think that the fractional derivative could be the inverse operation of the above fractional integral; however, this idea must be implemented with some care. In particular, note that even for an integer n, the differentiation ${{ \mathcal D }}^{n}$ and integration ${{ \mathcal I }}^{n}$ operators in general do not commute, and ${{ \mathcal D }}^{n}$ is only a left inverse of ${{ \mathcal I }}^{n};$ In fact, ${{ \mathcal D }}^{n}\,{{ \mathcal I }}^{n}f(x)=f(x)$ holds, but the converse is not necessarily true due to the appearance of numerous integration constants. Therefore, given an integer m, one defines the fractional derivative ${{ \mathcal D }}^{s}$ of order s ∈ (m − 1, m) in such a way that ${{ \mathcal D }}^{s}\equiv {{ \mathcal D }}^{m}\,{{ \mathcal I }}^{m-s}$, or, more explicitly,

$\begin{eqnarray}&&{{ \mathcal D }}^{s}f(x)=\displaystyle \frac{{d}^{m}}{{{dx}}^{m}}\left[\displaystyle \frac{1}{{\rm{\Gamma }}(m-s)}\,{\int }_{0}^{x}{dy}\,f(y){\left(x-y\right)}^{m-s-1}\right],\end{eqnarray} \tag{ A4 }$

with ${{ \mathcal D }}^{0}f(x)=f(x)$. Using this definition and the semigroup property, one has ${{ \mathcal D }}^{s}\,{{ \mathcal I }}^{s}f(x)={{ \mathcal D }}^{m}\,{{ \mathcal I }}^{m-s}\,{{ \mathcal I }}^{s}f(x)\,={{ \mathcal D }}^{m}\,{{ \mathcal I }}^{m}f(x)=f(x)$, so that ${{ \mathcal D }}^{s}$ is truly the left inverse of ${{ \mathcal I }}^{s}$. It is also trivial to verify the action on power functions

$\begin{eqnarray}&&{{ \mathcal D }}^{s}{x}^{a}=\displaystyle \frac{{\rm{\Gamma }}(a+1)}{{\rm{\Gamma }}(a-s+1)}\,{x}^{a-s}.\end{eqnarray} \tag{ A5 }$

The definition of a fractional derivative implies an inherently nonlocal behavior, since ${ \mathcal D }f(x)$ does not depend only on the values of the dependent variable infinitesimally close to x (as in standard calculus) but rather on a range of values even far away from x. More details on the basis of fractional calculus, including an alternative definition of fractional derivatives (e.g., the Caputo one), can be found in the books by Oldham & Spanier (1974) and Samko et al. (1993).

A.2. Fractional Laplacian

The standard Laplacian Δ is a linear, elliptic, second-order, local differential operator whose applications in physics include wave mechanics, gravity, electrodynamics, hydrodynamics, biophysics, and probability theory. However, some physical phenomena have been experimentally shown to exhibit a nonlocal behavior that has motivated the extension of the Laplacian to a fractional operator ${\left(-{\rm{\Delta }}\right)}^{s};$ among the many possible definitions, here we rely on the simplest, "spectral" one in Fourier space. In an n-dimensional space, the traditional Laplacian satisfies

$\begin{eqnarray}&&{ \mathcal F }\{-{\rm{\Delta }}f\}({\boldsymbol{k}})={k}^{2}\,{ \mathcal F }\{f\}({\boldsymbol{k}}),\end{eqnarray} \tag{ A6 }$

where ${ \mathcal F }\{\cdot \}$ denotes the Fourier transform and hereafter k ≡ ∣ k ∣. The fractional Laplacian of order s ∈ (0, n/2) is defined as the linear operator satisfying

$\begin{eqnarray}&&{ \mathcal F }\{{\left(-{\rm{\Delta }}\right)}^{s}f\}({\boldsymbol{k}})={k}^{2s}\,{ \mathcal F }\{f\}({\boldsymbol{k}}).\end{eqnarray} \tag{ A7 }$

Taking the inverse Fourier transform ${{ \mathcal F }}^{-1}\{\cdot \}$ of both sides, one obtains an explicit form for the fractional Laplacian in real space as

$\begin{eqnarray}&&{\left(-{\rm{\Delta }}\right)}^{s}\,f({\boldsymbol{x}})={{ \mathcal F }}^{-1}\{{k}^{2s}\,{ \mathcal F }\{f\}({\boldsymbol{k}})\}({\boldsymbol{x}}).\end{eqnarray} \tag{ A8 }$

In analogy with the one-dimensional fractional integral, one may wonder whether there is an operator that acts like the inverse of the fractional Laplacian. This is provided by the Riesz operator ${{ \mathcal R }}_{s}$, defined by

$\begin{eqnarray}&&{{ \mathcal R }}_{s}\,f({\boldsymbol{x}})\equiv \displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{n-s}{2}\right)}{{2}^{s}\,{\pi }^{n/2}\,{\rm{\Gamma }}\left(\tfrac{s}{2}\right)}\,\int {d}^{n}{\boldsymbol{x}}^{\prime} \,\displaystyle \frac{f({\boldsymbol{x}}^{\prime} )}{| {\boldsymbol{x}}-{\boldsymbol{x}}^{\prime} {| }^{n-s}};\end{eqnarray} \tag{ A9 }$

the above can be shown to satisfy

$\begin{eqnarray}&&{\left(-{\rm{\Delta }}\right)}^{s}\,{{ \mathcal R }}_{2s}\,f({\boldsymbol{x}})={{ \mathcal R }}_{2s}\,{\left(-{\rm{\Delta }}\right)}^{s}\,f({\boldsymbol{x}})=f({\boldsymbol{x}})\end{eqnarray} \tag{ A10 }$

so that one can naively set ${{ \mathcal R }}_{2s}={\left(-{\rm{\Delta }}\right)}^{-s}$.

More details on the definition and properties of the fractional Laplacian can be found in the classic book by Stein (1971) and the modern article by Lischke et al. (2020).

A.3. Fractional Poisson Equation

In astrophysics, the main relevance of the Laplacian arises in the context of Newtonian gravity; specifically, in three spatial dimensions, the Poisson equation links the gravitational potential Φ to a given DM density ρ as

$\begin{eqnarray}&&{\rm{\Delta }}\,{\rm{\Phi }}({\boldsymbol{r}})=4\pi G\,\rho ({\boldsymbol{r}}).\end{eqnarray} \tag{ A11 }$

The above constitutes a partial differential equation that, when solved with appropriate boundary conditions, yields the gravitational potential anywhere in space. It is worth stressing that the Poisson equation is local in nature, meaning that the potential at a given space location depends only on the density distribution at the same point.

The Green function ${ \mathcal G }$ of the Poisson equation is the solution of the equation ${\rm{\Delta }}\,{ \mathcal G }({\boldsymbol{r}})=4\pi G\,{\delta }_{{\rm{D}}}^{3}({\boldsymbol{r}})$ in terms of the Dirac delta function δ_D( · ); by Fourier transform, this has the well-known solution ${ \mathcal G }({\boldsymbol{r}})=1/r$. The relevance of ${ \mathcal G }$ is that the solution of the Poisson equation vanishing at infinity for a generic density distribution ρ( r ) can be written as the convolution integral

$\begin{eqnarray}&&{\rm{\Phi }}({\boldsymbol{r}})=\int {d}^{3}{\boldsymbol{r}}^{\prime} \,{ \mathcal G }({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} )\rho ({\boldsymbol{r}}^{\prime} )=-G\,\int {d}^{3}{\boldsymbol{r}}^{\prime} \,\displaystyle \frac{\rho ({\boldsymbol{r}}^{\prime} )}{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | };\end{eqnarray} \tag{ A12 }$

when applied to a pointlike distribution at the origin $\rho ({\boldsymbol{r}})=m\,{\delta }_{{\rm{D}}}^{3}({\boldsymbol{r}})$, this yields the famous potential Φ( r ) = −Gm/r.

To introduce some nonlocal behavior (which may be needed to solve small-scale problems in galactic dynamics, as we advocate in the main text), one can generalize the above to the fractional case. This is accomplished by the natural ansatz (see Giusti 2020),

$\begin{eqnarray}&&{\left(-{\rm{\Delta }}\right)}^{s}\,{\rm{\Phi }}({\boldsymbol{r}})=-4\pi G\,{{\ell }}^{2-2s}\,\rho ({\boldsymbol{r}}),\end{eqnarray} \tag{ A13 }$

where ℓ is a constant with physical dimensions of a length; the choice of −Δ is dictated by the need to have a positively defined operator so as to apply the properties discussed in the previous section.

The associated fractional Green function ${{ \mathcal G }}_{s}$ may be found as the solution of the equation ${\left(-{\rm{\Delta }}\right)}^{s}\,{{ \mathcal G }}_{s}({\boldsymbol{r}})=-4\pi G\,{{\ell }}^{2-2s}\,{\delta }_{{\rm{D}}}^{3}({\boldsymbol{r}})$. Applying the Riesz operator for s ∈ (0, 3/2), one finds that

$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal G }}_{s}({\boldsymbol{r}}) & = & -4\pi G\,{{\ell }}^{2-2s}\,{{ \mathcal R }}_{2s}\,{\delta }_{\ {\rm{D}}}^{3}({\boldsymbol{r}})\\ & = & -\displaystyle \frac{G\,{{\ell }}^{2-2s}\,{\rm{\Gamma }}\left(\tfrac{3}{2}-s\right)}{{4}^{s-1}\,\sqrt{\pi }\,{\rm{\Gamma }}(s)}\,\displaystyle \frac{1}{{r}^{3-2s}};\end{array}\end{eqnarray} \tag{ A14 }$

it is straightforward to check that ${{ \mathcal G }}_{s=1}({\boldsymbol{r}})$ collapses to the Newtonian Green function. The case of s = 3/2 is a bit tricky and may be solved more conveniently in Fourier space; there the Poisson equation reads ${ \mathcal F }\{{{ \mathcal G }}_{s}\}({\boldsymbol{k}})=-4\pi G/{\ell }\,{k}^{3}$, and thus ${{ \mathcal G }}_{s}=-(4\pi G/{\ell })\times {{ \mathcal F }}^{-1}\{1/{k}^{3}\}$. Given that $-{k}^{2}\,{ \mathcal F }\{\mathrm{ln}(r/{\ell })\}={ \mathcal F }\{{\rm{\Delta }}\mathrm{log}(r/{\ell })\}={ \mathcal F }\{1/{r}^{2}\}=2{\pi }^{2}/k$ (where the last equality is meant to hold in the principal value sense), one finds ${ \mathcal F }\{\mathrm{ln}(r/{\ell })\}=-2{\pi }^{2}/{k}^{3}$, and hence

$\begin{eqnarray}&&{{ \mathcal G }}_{s=3/2}({\boldsymbol{r}})=\displaystyle \frac{2G}{\pi \,{\ell }}\,\mathrm{ln}\left(\displaystyle \frac{r}{{\ell }}\right).\end{eqnarray} \tag{ A15 }$

The potential for a point mass m located at the origin goes like Φ_s ( r ) ∝ Gm ℓ^2−2s /r^3−2s for s ∈ (0, 3/2) and ${{\rm{\Phi }}}_{s=3/2}({\boldsymbol{r}})\propto (G\,m/{\ell })\times \mathrm{ln}(r/{\ell })$ for s = 3/2, thus decreasing faster than the Newtonian case (s = 1) for s ∈ (0, 1) and slower for s ∈ (1, 3/2].

Correspondingly, the rotational velocity v²(r) = r ∣dΦ/dr∣ of a test mass on a circular orbit in such a potential goes like ${v}_{s}{(r)\propto \sqrt{{Gm}/{\ell }}\times ({\ell }/r)}^{3/2-s}$. Curiously, for s = 3/2, it is found that ${v}_{s=3/2}(r)=\sqrt{2\,G\,m/\pi \,{\ell }}$ is constant, an occurrence that has triggered some connection with MOND theories (see Giusti 2020); in fact, since at low accelerations, MOND implies that v⁴ = Gm a₀ in terms of a universal constant a₀ ≈ 1.2 × 10⁻¹⁰ m s⁻², the two expressions are consistent provided that ${\ell }=(2/\pi )\times \sqrt{G\,m/{a}_{0}}$ applies. We are not adopting this MONDian viewpoint as our fiducial setup in the present paper, but we explore it somewhat in Section 4.

In this Appendix, we provide original analytic solutions Φ_s (r) of the fractional Poisson equation for various literature density profiles in spherical symmetry. To derive these, it is convenient to use the coordinates ${\boldsymbol{r}}^{\prime} =(r^{\prime} ,\theta ,\phi )$, placing the point r where the potential is calculated at the north pole; using Equation (A14) for s ∈ [1, 3/2), one has

$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{s}(r) & = & \int {d}^{3}{\boldsymbol{r}}^{\prime} \,{{ \mathcal G }}_{s}({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} )\rho (r^{\prime} )\\ & = & -\displaystyle \frac{G\,{{\ell }}^{2-2s}\,{\rm{\Gamma }}\left(\tfrac{3}{2}-s\right)}{{4}^{s-1}\,\sqrt{\pi }\,{\rm{\Gamma }}(s)}{\int }_{0}^{2\pi }d\phi \,{\int }_{0}^{\infty }{dr}^{\prime} \,{r}^{{\prime} 2}\,\rho (r^{\prime} )\\ & & \times \,{\int }_{0}^{\pi }d\theta \,\displaystyle \frac{\sin \theta }{{\left({r}^{2}+{r}^{{\prime} 2}-2\,r\,r^{\prime} \,\cos \theta \right)}^{3/2-s}}.\end{array}\end{eqnarray} \tag{ B1 }$

A straightforward calculation yields the solution

$\begin{eqnarray}&&{{\rm{\Phi }}}_{s}(r)=-\displaystyle \frac{\sqrt{\pi }\,G\,{{\ell }}^{2-2s}\,{\rm{\Gamma }}\left(\tfrac{3}{2}-s\right)}{{4}^{s-3/2}(2s-1){\rm{\Gamma }}(s)}\,\displaystyle \frac{{{ \mathcal J }}_{s}^{\rho }(r)}{r},\end{eqnarray} \tag{ B2 }$

where

$\begin{eqnarray}&&{{ \mathcal J }}_{s}^{\rho }(r)={\int }_{0}^{\infty }{dr}^{\prime} \,r^{\prime} \,\rho (r^{\prime} )\left[{\left(r+r^{\prime} \right)}^{2s-1}-| r-r^{\prime} {| }^{2s-1}\right].\end{eqnarray} \tag{ B3 }$

In the case of s = 3/2, one has to instead use Equation (A15), obtaining

$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{s=3/2}(r) & = & \int {d}^{3}{\boldsymbol{r}}^{\prime} \,{{ \mathcal G }}_{s=3/2}({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} )\rho (r^{\prime} )\\ & = & \displaystyle \frac{2G}{\pi {\ell }}\,{\int }_{0}^{2\pi }d\phi \,{\int }_{0}^{\infty }{dr}^{\prime} \,{r}^{{\prime} 2}\,\rho (r^{\prime} )\\ & & \times \,{\int }_{0}^{\pi }d\theta \,\sin \theta \,\mathrm{ln}\left(\displaystyle \frac{\sqrt{{r}^{2}+{r}^{{\prime} 2}-2\,r\,r^{\prime} \,\cos \theta }}{{\ell }}\right),\end{array}\end{eqnarray} \tag{ B4 }$

which yields the integral representation

$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{s=3/2}(r)=\displaystyle \frac{2G\,{\ell }}{r}\,{\int }_{0}^{\infty }{dr}^{\prime} \,r^{\prime} \,\rho (r^{\prime} )\\ \,\times \,\left[{\left(\displaystyle \frac{r+r^{\prime} }{{\ell }}\right)}^{2}\,\mathrm{ln}\left(\displaystyle \frac{r+r^{\prime} }{{\ell }}\right)-{\left(\displaystyle \frac{r-r^{\prime} }{{\ell }}\right)}^{2}\,\mathrm{ln}\left(\displaystyle \frac{r-r^{\prime} }{{\ell }}\right)-2\,\displaystyle \frac{r\,r^{\prime} }{{{\ell }}^{2}}\right].\end{array}\end{eqnarray} \tag{ B5 }$

Another method of solving the fractional Poisson equation for spherically symmetric distribution is to work in the Fourier domain, where the density can be represented as

$\begin{eqnarray}&&{ \mathcal F }\{\rho \}(k)=\displaystyle \frac{4\pi }{k}\,{\int }_{0}^{\infty }{dr}\,r\,\rho (r)\sin ({kr}).\end{eqnarray} \tag{ B6 }$

Then, by inverse Fourier transform, the fractional potential for s ∈ [1, 3/2] is given by

$\begin{eqnarray}&&{{\rm{\Phi }}}_{s}(r)=-\displaystyle \frac{2G\,{{\ell }}^{2-2s}}{\pi \,r}\,{\int }_{0}^{\infty }{dk}\,\displaystyle \frac{\sin ({kr})}{{k}^{2s-1}}\,{ \mathcal F }\{\rho \}(k);\end{eqnarray} \tag{ B7 }$

in particular, such a Fourier method turns out to be very useful to solve the s = 3/2 case.

Below, we will derive the solutions for the following density profiles: NFW, Hernquist, Plummer, and razor-thin exponential disk. Note that in the main text, we mainly consider the NFW and exponential disk density profiles; however, the other solutions are presented for completeness and possible exploitation from the interested scientific community.

B.1. NFW Profile

The NFW (Navarro et al. 1997) density profile, which is routinely exploited to model the DM distribution within virialized halos, is defined as

$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{\rho }_{s}\,{r}_{s}^{3}}{r{\left(r+{r}_{s}\right)}^{2}}\end{eqnarray} \tag{ B8 }$

in terms of a scale density ρ_s and scale radius r_s ; note that the cumulative mass is logarithmically divergent at large distances, so a finite halo boundary must be introduced (usually taken to be the virial radius). The Newtonian potential can be easily computed:

$\begin{eqnarray}&&{{\rm{\Phi }}}_{s=1}(r)=-\displaystyle \frac{4\pi G\,{\rho }_{s}\,{r}_{s}^{3}}{r}\,\mathrm{log}\left(1+\displaystyle \frac{r}{{r}_{s}}\right).\end{eqnarray} \tag{ B9 }$

In the fractional case, one has to solve Equation (B2) in terms of the integral

$\begin{eqnarray}&&{{ \mathcal J }}_{s}^{\rho }(r)={\rho }_{s}\,{r}_{s}^{3}\,{\int }_{0}^{\infty }{dr}^{\prime} \,\displaystyle \frac{{\left(r+r^{\prime} \right)}^{2s-1}-| r-r^{\prime} {| }^{2s-1}}{{\left(r^{\prime} +{r}_{s}\right)}^{2}},\end{eqnarray} \tag{ B10 }$

which converges for s ∈ [1, 3/2) to

$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{s}(r) & = & -\sqrt{\pi }\,G\,{\rho }_{s}\,{r}_{s}^{2}{\left(\displaystyle \frac{2\,{\ell }}{{r}_{s}}\right)}^{2-2s}\,\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{3}{2}-s\right)}{{\rm{\Gamma }}(s+1)}\,\displaystyle \frac{{r}_{s}}{r}\\ & & \times \left\{\displaystyle \frac{2\pi s}{\sin (2\pi s)}\left[{\left(1+\displaystyle \frac{r}{{r}_{s}}\right)}^{2s-2}-{\left(1-\displaystyle \frac{r}{{r}_{s}}\right)}^{2s-2}\right]+\right.\\ & & +\displaystyle \frac{{\left(r/{r}_{s}\right)}^{2s}}{1-{\left(r/{r}_{s}\right)}^{2}}\left[{\left(1+\displaystyle \frac{r}{{r}_{s}}\right)}_{2}{F}_{1}\left(1,1,2s+1,\displaystyle \frac{r}{{r}_{s}}\right)\right.\\ & & +\left.\left.{\left(1-\displaystyle \frac{r}{{r}_{s}}\right)}_{2}{F}_{1}\left(1,1,2s+1,-\displaystyle \frac{r}{{r}_{s}}\right)-\displaystyle \frac{4s}{2s-1}\right]\right\}.\end{array}\end{eqnarray} \tag{ B11 }$

in the above, ₂ F₁ is the ordinary hypergeometric function, which is defined by the power series ${}_{2}{F}_{1}{(a,b,c;x)\equiv {\sum }_{k=0}^{\infty }{(a)}_{k}{(b)}_{k}\,{x}^{k}/(c)}_{k}\,k!$ in terms of the Pochammer symbols (q)₀ = 1 and (q)_k = q(q + 1)...(q + k − 1) for any positive integer k.

The case of s = 3/2 is more easily treated in the Fourier domain and requires regularizing the integrals appearing in Equation (B7); the result reads

$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{s=3/2}(r) & = & -\displaystyle \frac{4\,G\,{\rho }_{s}\,{r}_{s}^{3}}{{\ell }}\,\displaystyle \frac{{r}_{s}}{r}\left\{2\,\displaystyle \frac{r}{{r}_{s}}\left[\mathrm{log}\left(\displaystyle \frac{r}{{r}_{s}}\right)-1\right]\right.\\ & & -\left(1+\displaystyle \frac{r}{{r}_{s}}\right)\mathrm{log}\left(\displaystyle \frac{r}{{r}_{s}}\right)\mathrm{log}\left(1+\displaystyle \frac{r}{{r}_{s}}\right)+\\ & & +\left(\displaystyle \frac{r}{{r}_{s}}-1\right){\mathrm{Li}}_{2}\left(1-\displaystyle \frac{r}{{r}_{s}}\right)\\ & & \left.-\left(1+\displaystyle \frac{r}{{r}_{s}}\right){\mathrm{Li}}_{2}\left(-\displaystyle \frac{r}{{r}_{s}}\right)+\displaystyle \frac{{\pi }^{2}}{6}\right\}\end{array}\end{eqnarray} \tag{ B12 }$

in terms of the dilogarithm function ${\mathrm{Li}}_{2}(x)\equiv {\sum }_{k\,=\,1}^{\infty }{x}^{k}/{k}^{2}$.

The fractional NFW potential is extensively exploited in the main text (see Section 2); the corresponding (suitably normalized) RC v² = r ∣dΦ/dr∣ is illustrated in Figure 1.

B.2. Hernquist Profile

The Hernquist (1990) density profile, which is often used to model the stellar distribution within galactic bulges, is defined as

$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{M}_{\infty }}{2\pi }\,\displaystyle \frac{a}{r{\left(r+a\right)}^{3}}\end{eqnarray} \tag{ B13 }$

in terms of a scale radius a and the total mass M_∞. The Newtonian potential can be easily computed:

$\begin{eqnarray}&&{{\rm{\Phi }}}_{s=1}(r)=-\displaystyle \frac{{{GM}}_{\infty }}{a}\,\displaystyle \frac{1}{1+r/a}.\end{eqnarray} \tag{ B14 }$

In the fractional case, one has to solve Equation (B2) in terms of the integral

$\begin{eqnarray}&&{{ \mathcal J }}_{s}^{\rho }(r)=\displaystyle \frac{{M}_{\infty }\,a}{2\pi }\,{\int }_{0}^{\infty }{dr}^{\prime} \,\displaystyle \frac{{\left(r+r^{\prime} \right)}^{2s-1}-| r-r^{\prime} {| }^{2s-1}}{{\left(r^{\prime} +a\right)}^{3}},\end{eqnarray} \tag{ B15 }$

which converges for s ∈ [1, 3/2) to

$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{s}(r)=-\displaystyle \frac{{{GM}}_{\infty }}{a}\,\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{3}{2}-s\right)}{{4}^{s-1}\,\sqrt{\pi }(2s-1){\rm{\Gamma }}(s)}{\left(\displaystyle \frac{{\ell }}{a}\right)}^{2-2s}\\ \ \times \left\{\displaystyle \frac{(2s-3)(a/r)+(2s-1)(r/a)}{{\left(1-{r}^{2}/{a}^{2}\right)}^{2}}{\left(\displaystyle \frac{r}{a}\right)}^{2s}\right.\\ \ +\displaystyle \frac{\pi \,(2s-1)(s-1)}{\sin (2\pi s)}\,\displaystyle \frac{a}{r}\left[{\left(1-\displaystyle \frac{r}{a}\right)}^{2s-3}-{\left(1+\displaystyle \frac{r}{a}\right)}^{2s-3}\right]\\ \ -\,\displaystyle \frac{(2s-1)(s-1)}{2\,s\,{\left(1-{r}^{2}/{a}^{2}\right)}^{2}}{\left(\displaystyle \frac{r}{a}\right)}^{2s-1}\\ \ \times \,\left[{\left(1+\displaystyle \frac{r}{a}\right)}^{2}{\,}_{2}{F}_{1}\left(1,1,2s+1,\displaystyle \frac{r}{a}\right)\right.\\ \ \left.\left.+\,{\left(1-\displaystyle \frac{r}{a}\right)}^{2}{\ }_{2}{F}_{1}\left(1,1,2s+1,-\displaystyle \frac{r}{a}\right)\right]\right\}.\end{array}\end{eqnarray} \tag{ B16 }$

The case of s = 3/2 is more easily treated in the Fourier domain and requires regularizing the integrals appearing in Equation (B7); the result reads

$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{s=3/2}(r)=\displaystyle \frac{{{GM}}_{\infty }}{{\ell }}\,\displaystyle \frac{2}{\pi }\left\{\mathrm{ln}\left(\displaystyle \frac{r}{a}\right)+\displaystyle \frac{{\pi }^{2}}{12}\,\displaystyle \frac{a}{r}-\displaystyle \frac{1}{2}\,\displaystyle \frac{a}{r}\right.\\ \ \left.\times \ \left[\mathrm{ln}\left(\displaystyle \frac{r}{a}\right)\mathrm{ln}\left(1+\displaystyle \frac{r}{a}\right)+{\mathrm{Li}}_{2}\left(-\displaystyle \frac{r}{a}\right)+{\mathrm{Li}}_{2}\left(1-\displaystyle \frac{r}{a}\right)\right]\right\}.\end{array}\end{eqnarray} \tag{ B17 }$

The (suitably normalized) RC v² = r ∣dΦ/dr∣ corresponding to the Hernquist profile is illustrated in Figure B1 (left panel).

Zoom In Zoom Out Reset image size

B.3. Plummer Profile

The Plummer (1911) density profile, which is often exploited to model the stellar distribution within globular clusters, is defined as

$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{3\,{M}_{\infty }}{4\pi \,{a}^{3}}{\left(1+\displaystyle \frac{{r}^{2}}{{a}^{2}}\right)}^{-5/2}\end{eqnarray} \tag{ B18 }$

in terms of a scale radius a and the total mass M_∞. The Newtonian potential can be easily computed:

$\begin{eqnarray}&&{{\rm{\Phi }}}_{s=1}(r)=-\displaystyle \frac{{{GM}}_{\infty }}{a}\,\displaystyle \frac{1}{\sqrt{1+{r}^{2}/{a}^{2}}}.\end{eqnarray} \tag{ B19 }$

The fractional Poisson equation is most easily solved in Fourier space via Equation (B7); for s ∈ [1, 3/2), we obtain the expression

$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{s}(r) & = & \displaystyle \frac{{{GM}}_{\infty }}{a}\,\displaystyle \frac{{\rm{\Gamma }}(\tfrac{3}{2}-s){\rm{\Gamma }}(\tfrac{5}{2}-s)}{{2}^{2s-1}\,\pi (s-1)(s-2)}{\left(\displaystyle \frac{{\ell }}{a}\right)}^{2-2s}\\ & & \times \left\{{\left(1+\displaystyle \frac{{a}^{2}}{{r}^{2}}\right)}_{2}{F}_{1}\left(\displaystyle \frac{3}{2}-s,\displaystyle \frac{5}{2}-s,-\displaystyle \frac{1}{2},-\displaystyle \frac{{r}^{2}}{{a}^{2}}\right)\right.\\ & & -\left.\left[4(2-s)+\displaystyle \frac{{a}^{2}}{{r}^{2}}\right]{\,}_{2}{F}_{1}\left(\displaystyle \frac{3}{2}-s,\displaystyle \frac{5}{2}-s,+\displaystyle \frac{1}{2},-\displaystyle \frac{{r}^{2}}{{a}^{2}}\right)\right\}.\end{array}\end{eqnarray} \tag{ B20 }$

The case of s = 3/2, found after appropriate regularization, reads

$\begin{eqnarray}&&{{\rm{\Phi }}}_{s=3/2}(r)=\displaystyle \frac{{{GM}}_{\infty }}{{\ell }}\,\displaystyle \frac{2}{\pi }\,\sqrt{1+\displaystyle \frac{{a}^{2}}{{r}^{2}}}\,\mathrm{arcsinh}\left(\displaystyle \frac{r}{a}\right).\end{eqnarray} \tag{ B21 }$

The (suitably normalized) RC v² = r ∣dΦ/dr∣ corresponding to the Plummer profile is illustrated in Figure B1 (middle panel).

B.4. Razor-thin Exponential Disk

The density profile of a razor-thin exponential disk in cylindrical coordinates is defined as

$\begin{eqnarray}&&\rho (R,z)=\displaystyle \frac{{M}_{\infty }}{2\pi \,{r}_{d}^{2}}\,{e}^{-R/{r}_{d}}\,{\delta }_{{\rm{D}}}(z).\end{eqnarray} \tag{ B22 }$

The Newtonian potential on the plane of the disk (z = 0) reads

$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{s=1}(R,0) & = & -\displaystyle \frac{{{GM}}_{\infty }}{{r}_{d}}\,\displaystyle \frac{R}{2\,{r}_{d}}\left[{I}_{0}\left(\displaystyle \frac{R}{2\,{r}_{d}}\right){K}_{1}\left(\displaystyle \frac{R}{2\,{r}_{d}}\right)\right.\\ & & \left.-{I}_{1}\left(\displaystyle \frac{R}{2\,{r}_{d}}\right){K}_{0}\left(\displaystyle \frac{R}{2\,{r}_{d}}\right)\right].\end{array}\end{eqnarray} \tag{ B23 }$

The solution of the fractional Poisson equation may be most easily found in Fourier space after Equation (B7). For an index s ∈ [1, 3/2), we find

$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{s}(R,0)=-\displaystyle \frac{{{GM}}_{\infty }}{{r}_{d}}\,\displaystyle \frac{{\rm{\Gamma }}\left(s-1/2\right)}{2\,\sqrt{\pi }\,{\rm{\Gamma }}(s)}{\left(\displaystyle \frac{{\ell }}{{r}_{d}}\right)}^{2-2s}\\ \ \times \ \left[{2}^{2s}\,\displaystyle \frac{{\rm{\Gamma }}(2-2s){\rm{\Gamma }}(s)}{{\rm{\Gamma }}(1-s)}{\ }_{1}{F}_{2}\left(\displaystyle \frac{3}{2}-s;1,1-s,\displaystyle \frac{{R}^{2}}{4\,{r}_{d}^{2}}\right)\right.\\ \ -\ \left.\displaystyle \frac{\pi }{\sin (\pi s){{\rm{\Gamma }}}^{2}(s+1)}{\left(\displaystyle \frac{R}{2\,{r}_{d}}\right)}^{2s}{\ }_{1}{F}_{2}\left(\displaystyle \frac{3}{2};1+s,1+s,\displaystyle \frac{{R}^{2}}{4\,{r}_{d}^{2}}\right)\right].\end{array}\end{eqnarray} \tag{ B24 }$

in the above, ₁ F₂ is the generalized hypergeometric function, which is defined by the power series ${}_{1}{F}_{2}{(a,b,c;x)\equiv {\sum }_{k=0}^{\infty }{(a)}_{k}\,{x}^{k}/{(b)}_{k}(c)}_{k}\,k!$ in terms of the Pochammer symbols (q)₀ = 1 and (q)_k = q(q + 1)...(q + k − 1) for any positive integer k.

In the case of s = 3/2, after regularization of the integrals, we obtain

$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{s=3/2}(R,0) & = & -\displaystyle \frac{G\,{M}_{\infty }}{{\ell }}\,\displaystyle \frac{2}{\pi }\left[\mathrm{Ei}\left(-\displaystyle \frac{R}{{r}_{d}}\right)\right.\\ & & -\left.\,{e}^{-R/{r}_{d}}-\mathrm{ln}\left(\displaystyle \frac{R}{{r}_{d}}\right)\right]\end{array}\end{eqnarray} \tag{ B25 }$

in terms of the exponential integral function $\mathrm{Ei}(x)\equiv {\int }_{-\infty }^{x}{dt}\,{e}^{t}/t$.

The (suitably normalized) RC v² = R ∣dΦ/dR∣ corresponding to the razor-thin exponential disk profile is illustrated in Figure B1 (right panel).