Dark Matter in Fractional Gravity II: Tests in Galaxy Clusters

[abridged] Recently, in Benetti et al. (Astrophys. J. 2023, 949, 65), we suggested that the dark matter (DM) component in galaxies may originate fractional gravity. In such a framework, the DM component exists, but the gravitational potential associated to its density distribution is determined by a modified Poisson equation including fractional derivatives, which are meant to describe nonlocal effects. In Benetti et al., we showed that fractional gravity worked very well for reproducing the kinematics of disk-dominated galaxies, especially dwarfs; there is also preliminary evidence that the strength of fractional effects tends to weaken toward more massive systems. Here, we aim to test fractional gravity in galaxy clusters, with a twofold aim: (i) perform an independent sanity check that it can accurately describe such large and massive structures; (ii) derive a clear-cut trend for its strength in systems with different DM masses. To this purpose, we forward model the density and pressure distributions of the intracluster medium (ICM), working out the hydrostatic equilibrium equation in fractional gravity. Then, we perform a Bayesian analysis of the X-COP galaxy cluster sample and infer constraints on the fractional gravity parameters, for individual clusters as well as stacked clusters. We find that fractional gravity performs remarkably well in modeling the ICM profiles for the X-COP sample. We also confirm the weakening of the fractional gravity effects toward more massive systems and derive the overall scaling of the fractional gravity parameters from dwarf galaxies to massive clusters, spanning six orders of magnitude in DM mass. Such an overall trend implies that fractional gravity can substantially alleviate the small-scale issues of the standard DM paradigm, while remaining successful on large cosmological scales.


Introduction
Galaxy clusters constitute the largest bound structures in the Universe, with dark matter (DM) masses M ∼ 10 14-15 M ⊙ and sizes extending out to R ∼ a few Mpcs. Most of the baryons are in the form of a hot diffuse gas, referred to as the intracluster medium (ICM), with a mass ratio over the DM very close to the cosmic fraction Ω b /Ω M ≈ 0.16 [2].
The density n(r) and temperature T(r) distributions of the ICM throughout the cluster can be probed thanks to the copious X-ray powers L X ∝ n 2 √ T R 3 ∼ 10 44-46 erg s −1 emitted by the ICM via thermal Bremsstrahlung and high-excitation lines [3,4]. The inferred high average temperatures k B T ∼ several keVs and low average number densities n ∼ 10 −3 cm −3 make the ICM the best plasma in the Universe ever, with thermal to electrostatic energy ratios k B T/e 2 n 1/3 ∼ 10 12 .
In addition, the pressure distribution p(r) can be probed thanks to the Sunyaev-Zel'dovich (SZ; [5,6]) effect, arising when the hot ICM electrons Compton upscatter the CMB photons crossing the cluster, tilting the latter's black-body spectrum toward high energies. In the microwave band, such a tilt mimics a diminution of the CMB temperature proportional to the Comptonization parameter y ∝ dℓ p(r), which encompasses the line-of-sight integral of the pressure profile. Combining X-ray and SZ data allows one to reconstruct the ICM thermodynamic profiles throughout most of the cluster volume, from the center to a few times R 500 or even beyond the virial boundary 1 .
In massive and sufficiently relaxed clusters, the ICM is expected to settle in hydrostatic equilibrium within the overall gravitational potential well mainly provided by the DM component. Under this assumption, the gas density profile reconstructed from X-rays and the gas pressure profile from SZ data can be combined to probe the shape of the DM gravitational potential and check whether this is consistent with the DM density run extracted from N-body simulations in the ΛCDM cosmology. This is the rationale of many investigations aimed at exploiting galaxy clusters to probe modified gravity scenarios [7][8][9][10][11][12][13], which have been developed to solve cosmological problems such as the origin of dark energy [14][15][16], and/or to alleviate small-scales issues of the standard cold DM paradigm [17][18][19][20][21][22]. In the latter vein, a prototypical example of such theories is the modified Newtonian dynamics (MOND) framework, which was originally designed 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 was not included, and baryons were the only source of the gravitational field. Although MOND can properly fit galactic RCs [23,24], its performances at the scales of galaxy clusters are somewhat debated [25,26].
More connected with the present work, in the last few years, various authors have put forward the idea that fractional calculus (i.e., the field of mathematics dealing with differentiation and integration of noninteger order) could be exploited to formulate modified gravity theories [27][28][29][30][31][32][33][34]. A relevant example is the theory of Newtonian fractional-dimensional gravity by [29,30], which introduces a generalized law of Newtonian gravity in a spatial dimension smaller than three, representing the local effective Hausdorff dimension of the matter distribution. Another approach by [27,28] relies on multifractional spacetimes with variable Hausdorff and spectral dimensions directly inspired from quantum gravity theories. The framework by [31,32] directly modifies the Laplacian operator in the Poisson equation to alter the dynamics followed by a test particle in a given gravitational well; a similar route is followed by [34], using fractional Fourier derivatives. All these theories adopt a MONDian viewpoint where DM is not present, and the galaxy kinematics is interpreted as a pure geometrical effect.
Recently, in [1], we suggested that the DM component itself may originate fractional gravity. In such a framework, the DM component exists, but the gravitational potential associated to its density distribution is determined by a modified Poisson equation including fractional derivatives (i.e., derivatives of noninteger type), which are meant to describe nonlocal effects; as such, this scenario is substantially different from the above theories where baryonic matter emulates DM-like effects via modifications of gravity. In [1], we showed that DM in fractional gravity worked very well for reproducing the kinematics of disk-dominated galaxies, especially dwarfs. In addition, we found preliminary evidence that the strength of fractional effects tends to weaken toward more massive systems; however, the latter finding is still subject to large uncertainties since the rotation curves of massive spirals were not probed out to radii large enough for the DM contribution to clearly emerge.
In the present work, we aim to extend our previous investigation to much larger scales and test fractional gravity in galaxy clusters. Our aim is twofold: (i) perform an independent sanity check that it can accurately describe the distributions of the ICM in clusters; (ii) derive a clear-cut trend for the strength of its effects over an extended DM mass range, from dwarf galaxies to galaxy clusters. To this purpose, we forward model the density and pressure distributions of the ICM, working out the hydrostatic equilibrium equation in fractional gravity. Such theoretical framework is then compared with data from the XMM-Newton Cluster Outskirts Project (X-COP 2 ; [35][36][37][38]), which consists of 12 clusters with well-observed X-ray and SZ data, providing density and pressure profiles over an extended radial range of ∼0.2-2 Mpc. We then perform a Bayesian analysis of the thermodynamic profiles of the X-COP sample and infer constraints on the fractional gravity parameters, for individual clusters and also for clusters stacked together.
The structure of the paper is straightforward: in Section 2, we describe our methods and analysis; in Section 3, we present and discuss our results; in Section 4, we summarize our findings and highlight future perspectives. Throughout the work, we adopt the standard, flat ΛCDM cosmology [39] with rounded parameter values: a matter density

Theoretical Background and Data Analysis
In this section, we recall the basics of the fractional gravity framework, illustrate how this can be exploited to derive the pressure profile of the ICM in hydrostatic equilibrium, and describe our Bayesian analysis to constrain the fractional gravity parameters.

DM in Fractional Gravity
The density distribution of virialized halos for collisionless DM as extracted from N-body simulations in the standard ΛCDM model is routinely described via the Navarro-Frenk-White profile [40]: where r s is a scale radius and ρ s a characteristic density. The associated cumulative mass is given by with M s ≡ 4π ρ s r 3 s . In the standard (Newtonian) case, the potential Φ N (r) associated to a given density distribution ρ(r) is computed from the Poisson equation supplemented with appropriate boundary conditions (usually taken as a vanishing potential at infinity): where ∆ is the Laplacian operator; this is an inherently local equation, in that the potential at a point depends only on the value of the density there. For the spherically symmetric NFW profile, one easily finds that from the above expressions of the mass and potential, it is straightforward to verify that |dΦ N /dr| = G M(< r)/r 2 , as a direct consequence of Birkhoff's theorem.
In fractional gravity, the potential Φ F (r) is instead derived from the modified Poisson equation [31] ( where (−∆) s is the fractional Laplacian operator (see [1,31] for details), s ∈ [1, 3/2] is the fractional index (this range of values for s is required to avoid divergences; see Appendix A in [1]), and ℓ is a fractional length scale that must be introduced for dimensional reasons. At variance with the standard case, the fractional Laplacian is inherently nonlocal; the index s measures the strength of this 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 standard case).
In [1], we solved the fractional Poisson equation sourced by the NFW density distribution. For s ∈ [1, 3/2), the solution reads with Γ(s) = ∞ 0 dx x s−1 e −x being the Euler Gamma function and 2 F 1 (a, b, c; x) = ∑ ∞ k=0 (a) k (b) k x k /(c) k k! being the ordinary hypergeometric function in terms of the Pochammer symbols (q) k defined as (q) 0 = 1 and (q) k = q (q + 1) . . . (q + k − 1); plainly, Φ F (r) for s = 1 coincides with the usual expression Φ N (r) of Equation (4). For the limiting case s = 3/2, the computation requires some principal-value regularization and the solution reads with Li 2 (x) = ∑ ∞ k=1 x k /k 2 being the dilogarithm function. Being a nonlocal framework, in fractional gravity, the Birkhoff theorem does not hold, but one can insist in writing |dΦ F /dr| = G M F (< r)/r 2 in terms of an effective mass M F (< r), which plainly will be a function of the fractional gravity parameters s and ℓ. We illustrate the effective mass profile in Figure 1, suitably normalized so as to remove the dependence of dimensional quantities (including ℓ), for different values of the fractional index s. With s increasing from unity, the effective mass profile steepens: in the inner region, a uniform sphere behavior (corresponding to a cored density profile) tends to be enforced, while in the outskirts the effective profile resembles that of an isothermal sphere. Note that all the normalized mass profiles intersect at very close values of r/r s ; more in detail, the profile with a given s crosses the one with s = 1 at r/r s ≈ 1.58 for s = 1.1, at r/r s ≈ 1.49 for s = 1.3, and at r/r s ≈ 1.36 for s = 1.5; plainly, in log scale, all these points appear clustered around log r/r s ≈ 0.15 and are barely discernible by eye.
To have a quantitative grasp on the overall effect of fractional gravity, consider the s = 3/2 case where the effective mass can be computed in terms of a relatively simple analytical expression; it reads and it is easily found to behave as M F (< r) ∝ [1 − 3 log(r/r s )] r 3 for r ≪ r s and as M F (< r) ∝ r ln(r/r s ) for r ≫ r s ; besides minor logarithmic corrections, the overall behavior is very similar to that of a cored isothermal sphere.

Forward Modeling of the ICM Thermodynamics
Assuming hydrostatic equilibrium and spherical symmetry, the distribution of the ICM in the overall gravitational potential well is ruled by the equation where Φ = Φ DM + Φ gas is the total gravitational potential with main contributions from DM and gas, ρ gas is the gas mass density, and P gas is the gas pressure. One can conveniently write ρ gas = µm p n gas in terms of the mean molecular weight µ ≈ 0.6 and of the gas number density n gas , which is in turn easily related to the electron density by the expression n gas ≈ 1.8 n e , applying for a fully ionized plasma at high temperatures and a subsolar chemical composition typical of the ICM. The observed electron density profile n e (r) of individual clusters inferred from X-ray observations is often empirically rendered by the (simplified version of the) Vikhlinin profile [41] n e (r) where n 0 is the central density, r c and r t are a core and a transition radius (r c < r t ), α, β, and ϵ < 5 are three slopes characterizing the inner, intermediate, and outer radial behavior.
The gas mass can then be computed as M gas (< r) = 4π r 0 dr ′ r ′2 ρ gas (r ′ ) and the gas contribution to the hydrostatic balance is fully specified by |dΦ gas /dr| = G M gas (< r)/r 2 .
As to the DM contribution, we can exploit the results of the previous section and write |dΦ DM /dr| = G M F (< r)/r 2 in terms of the fractional gravity's effective mass M F (< r) illustrated in Figure 1, which depends on the parameters s and ℓ; in the standard case (corresponding to s = 1), this is just the DM mass profile of Equation (2). The mass profile is also a function of the NFW scale radius r s and mass M s ; for the present analysis, it is convenient to trade off these parameters for the mass M 500 and the concentration c 500 ≡ R 500 /r s at the reference radius R 500 . The conversion between these variables can be performed easily using the relations M 500 = 4π 500 ρ c (z) R 3 500 /3 and M 500 = M s /[ln(1 + c 500 ) − c 500 /(1 + c 500 )] stemming from the definition of R 500 and from the adopted NFW mass distribution.
Then, the solution to the hydrostatic equilibrium equation is given by where the zero pressure at infinity has been taken as a boundary condition. 3 Observationally, X-ray surface brightness and spectroscopic data can probe the electron density n e (r) and the temperature T gas (r) profiles, whence the pressure profile P gas (r) ∝ n gas (r) T gas can be derived, although sensitivity and background issues make such a determination robust only in the region out to R 500 . In the outskirts, SZ observations can complement X-ray data in probing the pressure profile, though with some caveats about conversion from line-of-sightintegrated to spherically averaged quantities. The rationale of the above forward modeling of the hydrostatic equilibrium is to test the fractional gravity parameters entering in the effective mass profile M F (< r) by simultaneously fitting the observed electron density profile via Equation (10) and the observed pressure profile via Equation (11).

Bayesian Data Analysis
We tested the fractional gravity framework by exploiting the X-COP sample [35][36][37][38] of 12 massive galaxy clusters. The clusters are in the redshift range 0.04 ≲ z ≲ 0.1 and feature typical sizes R 500 ∼ 1-1.5 Mpc and masses M 500 ∼ 10 14 -10 15 M ⊙ . The X-COP clusters were selected to allow a robust reconstruction of the electron density and gas pressure profiles out to R 200 via a joint analysis of high-quality X-ray data from XMM-Newton and of high signal-to-noise SZ observations from Planck. Another important property of the sample is that the hydrostatic equilibrium holds to a high accuracy, with at most mild levels on nonthermal pressure support in the outermost regions, as demonstrated by the analysis of [36,38]; this is particularly important, since nonthermal effects can appreciably affect the mass estimation in the outer regions [42,43] and potentially induce spurious effects in constraining modified gravity parameters [44]. All in all, X-COP is currently the largest cluster sample available so far for robust mass-modeling studies over an extended radial range, and as such it has been extensively exploited to probe modified-gravity scenarios [10][11][12][13].
To estimate the parameters θ F ≡ (s, ℓ, c 500 , M 500 ) describing the effective mass profile M F (< r), alongside with those θ n e ≡ (n 0 , α, β, ϵ, r c , r t ) describing the electron density profile n e (r), we adopted a Bayesian framework and built the joint log-likelihood log L(θ) = log L P gas (θ F , θ n e ) + log L n e (θ n e ) .
Each term in the log-likelihood reads log L(θ) = −χ 2 (θ)/2, where the chi-square χ 2 (θ) = ∑ i [M(θ, r i ) − D(r i )] 2 /σ 2 D (r i ) was obtained by comparing our empirical model expectations M(θ, r i ) to the data values D(r i ) with their uncertainties σ D (r i ), summing over the different radial coordinates r i of the data (approximately 65 points for n e and 20 points per P gas , with small variations around these numbers from cluster to cluster); note that for the pressure data from SZ observations, we took into account the full covariance matrix.
We adopted flat priors π(θ) on all the parameters; specifically, for those entering the effective mass profile in fractional gravity we took s ∈ [1, 3/2], log ℓ (Mpc) ∈ [−3, 3], log c 500 ∈ [−2, 2], log M 500 (M ⊙ ) ∈ [13,16]. We then sampled the parameter posterior distributions P (θ) ∝ L(θ) π(θ) via the MCMC Python package emcee [45], running it with 10 4 steps and 200 walkers for every individual cluster; each walker was initialized with a random position uniformly sampled from the (flat) priors. After checking the autocorrelation time, we removed the first 20% of the flattened chain to ensure burn-in; the typical acceptance fractions of the various runs were in the range 30-40%.

Results
In Figures 2 and 3, we illustrate the outcome of the fitting procedure on the 12 individual pressure and density profiles of the X-COP sample. In each panel, the best fit (solid lines) and the 2σ credible intervals sampled from the posterior (shaded areas) are shown. The reduced χ 2 r value of the joint fit to the pressure and density profiles is also reported in Figure 2. Overall, the fits in the fractional gravity framework are very good. In a few cases (such as A3266 and A2319), the reduced χ 2 r is somewhat large, but this should not raise any alarm, since the outcome is caused by some peculiar feature in the density profile reconstructed from X-ray data (oscillation in the data points at intermediate radii) or because of some outlier data in the pressure profile reconstructed from SZ (especially in the innermost or outermost radii); note that we retained all data points in our analysis, including them in the reduced χ 2 r computation. In Figure 4, we illustrate the MCMC posterior distributions for two representative clusters in the sample, namely A2255 and ZW1215; for clarity, we restricted the plot to the subspace of parameters entering the effective mass profile. Magenta/contour lines display the results in our fiducial setup, where no mass prior was imposed; the white cross marks the best-fit value of the parameters.
The corner plots illustrate a clear degeneracy between the fractional length-scale parameter ℓ and the DM mass M 500 . This is somewhat expected since the effective mass profile entering the hydrostatic equilibrium equation scales like M 500 ℓ 2−2s . Therefore, it is possible to obtain the same normalization of the pressure profile, at a given density profile, by changing M 500 and ℓ in the same direction. Since s does not deviate much from unity, the ℓ dependence is weak, implying that to compensate a rather small change in mass requires a substantial variation in ℓ; on the other hand, this is also at the origin of the rather loose constraints that can be derived on the parameter ℓ with the present cluster sample.
The situation is expected to improve if a mass prior from other probes such as weak lensing (WL) is introduced in the analysis. However, one must be careful and use WL mass estimates that are independent from assumptions on the shape of the lensing potential; this is because in fractional gravity, the lensing potential corresponding to a given mass distribution would be different from the standard case, thus causing an inconsistency. For five X-COP clusters (A85, A1795, A2029, A2142, and ZW1215), such nonparametric WL mass determinations are available in the literature [46]. 4 The outcome of exploiting the WL mass prior on the marginalized distributions of the parameters is illustrated by the cyan contours/lines in Figure 4. The DM mass posterior estimate of ZW1215 is made considerably more precise, and as a consequence of the above degeneracy, the estimate of the fractional length scale ℓ is also appreciably tightened. In any case, the posterior distributions on all the parameters for the analysis without and with the WL mass prior are consistent within 1σ.
We also tested the performance of fractional gravity by stacking the X-COP data of all the clusters in the sample. Specifically, we built stacked electron density and pressure profiles by normalizing the individual profiles of the 12 clusters at a reference radius R 500 , by co-adding them in radial bins of normalized radii r/R 500 , and by computing the corresponding mean and standard deviation. The outcome of this procedure is illustrated in Figure 5: the crosses mark the stacked profiles, and for reference, the gray lines show the individual ones. All in all, the fractional gravity frameworks fit the stacked profiles to a remarkable accuracy. Table 1. Marginalized posterior estimates (mean and 68% confidence limits are reported) for the parameters from the MCMC analysis of the X-COP sample in fractional gravity. Columns report the values of the fractional parameter s, the fractional length scale ℓ, the DM mass M 500 , the halo concentration c 500 , and the reduced χ 2 r of the joint fit to the density and pressure profiles. The different portions of the table refer to the fit to individual clusters with no mass prior, to individual clusters with weak lensing mass priors (marked WL, and only available for five clusters from [46]), and to the stacked sample. For reference, the last column reports the best-fit DM mass M 500 from [38].  Fits to the individual pressure profiles of the X-COP clusters in fractional gravity, according to the Bayesian analysis described in Section 2. Circles refer to X-ray and squares to SZ data. Solid lines illustrate the median, and the shaded areas show the 2σ credible interval from sampling the posterior distribution of the fits. The vertical dotted lines mark the reference radius R 500 . The reduced χ 2 r value of the joint fit to the pressure and density profiles is reported in each panel. Figure 6 summarizes the posterior distributions of the fractional index s, fractional length scale ℓ, concentration c 500 , and DM mass M 500 . Table 1 reports the marginalized posterior estimates (mean and 1σ credible intervals) of these parameters for all the individual X-COP clusters (including the WL mass prior when available), and for the stacked sample. On average, it is seen that the deviations of the fractional index s from unity are modest in clusters, and this originates rather loose constraints on the length scale ℓ. The inferred values of the DM mass M 500 and concentration c 500 are reasonable and consistent with that estimated by a variety of other methods in standard gravity [38]; we also checked that the same agreement applied for the gas fraction, as expected given the very good fits to the gas density profiles.

Cluster
In Figure 7, we checked the concentration vs. the DM mass relation for the X-COP sample in fractional gravity. To fairly compare with the relation expected from N-body simulations in the ΛCDM framework, we converted our fitting variables c 500 and M 500 at a reference radius R 500 to the corresponding values c 200 and M 200 at R 200 ; this is a trivial rescaling given the adopted NFW density profile. In Figure 7, we show as filled magenta circles the outcome for individual X-COP clusters and with a magenta cross that for the stacked sample. It is seen that the estimates of c 200 and M 200 in fractional gravity are fairly consistent in shape and scatter with the concentration vs. mass relation extracted from N-body simulations [47].  In passing, we note that the clusters A644, A1644, A2255, and A2319 have been shown not to favor an NFW mass profile, but rather a cored Burkert-like one (e.g., a Burkert, Hernquist, or pseudo-isothermal distribution) [36]. When forward modeling the pressure profiles in standard gravity with the NFW density distribution, this causes inconsistent results (especially in mass and concentration values) and/or a poor fit [10]. Contrariwise, such values and fits in fractional gravity stay reasonable and good, since the mass profile entering the hydrostatic equilibrium equation is not the true NFW mass, but the effective mass, which, as mentioned in Section 2, mirrors that of a cored profile. For these four clusters, we also checked that using a cored Burkert-like density distributions in place of the NFW one as an input in our fractional gravity framework did not substantially improve the fits to pressure profiles, and rather forced the fractional index to values s ≈ 1 compatible with pure Newtonian gravity. In fact, fractional gravity actually reconciles the NFW density distribution from simulations with the observed galactic dynamics, which are empirically described via cored, Burkert-like profiles. Moreover, A2319 have been shown to be characterized by an appreciable nonthermal support in the outskirts [48], which causes some difficulties in forward modeling and fitting the pressure profiles via the usual hydrostatic equilibrium equation in standard gravity. Instead, curiously, in fractional gravity, the fits stay good, suggesting that such a nonlocal framework may constitute an effective rendition for the effects of a nonthermal support on the pressure distribution.   In Figure 8, we explore the scaling of the fractional gravity parameters with the DM mass. For this purpose, we put together the analysis of the X-COP clusters from this work, and the constraints coming from the fitting of stacked galaxy rotation curves by [1]. These joint datasets covered six orders of magnitude in DM mass from M 200 ∼ 10 9 M ⊙ to 10 15 M ⊙ . As to the fractional index s, we confirmed the decreasing trend with the DM mass, passing from values around s ≈ 1.4 in dwarf galaxies, to s ∼ 1.2 − 1.3 in intermediate mass galaxies, to s ∼ 1.1 in massive galaxies and clusters. We described the s vs. M 200 relation by a linear fit (dashed line) with shape s = a + b (log M 200 (M ⊙ ) − 11) via an orthogonal distance regression (ODR) algorithm that took into account the error bars on both axis; we obtained the best-fit parameters a = 1.24 ± 0.02 and b = −0.057 ± 0.006 and a reduced χ 2 r ≈ 1.87; a nonlinear fit (solid line) s = (5/4) interpolating between asymptotic values s = 1 and 1.5 at small and large masses yielded the bestfit parameters c = −0.39 ± 0.06, d = 10.76 ± 0.25 and a reduced χ 2 r ≈ 1.34. As to the fractional length scale, there was an increasing trend with the DM mass, extending the finding by [1] at the cluster scales. We fit the ℓ vs. M 200 relation with a linear shape ℓ (Mpc) = a + b (log M 200 (M ⊙ ) − 11) via an ODR algorithm, to obtain the best-fit parameters a = −2.66 ± 0.09, b = 0.66 ± 0.06 and a reduced χ 2 r ≈ 1.09. This relation was somewhat steeper than the scaling with the DM mass of the NFW scale radius r s , in such a way that in dwarf galaxies ℓ/r s ≈ 0.25 but this ratio increased to around one at the cluster scales.  Scaling relations Figure 8. Scaling relations in fractional gravity from galaxies to clusters: fractional index s (left) and length-scale ℓ (right) vs. DM mass M 200 . Magenta circles refer to individual clusters and magenta crosses to the stacked cluster sample from this work, while cyan stars refer to the stacked rotation curve of galaxies from [1]. Dashed lines display an ODR linear fit to the overall data, while a solid line in the left panel shows a nonlinear ODR fit with limiting values 1 and 1.5 from large to small masses, and a dotted line in the right panel illustrates the scale radius r s of the NFW profile (see Section 3 for further details).

Summary and Outlook
Extending the analysis carried out by [1] on galactic scales, in this paper, we tested fractional gravity in galaxy clusters. Our aim was twofold: (i) to perform an independent sanity check that fractional gravity can accurately describe such large and massive structures; (ii) to derive a clear-cut trend for the strength of fractional gravity effects in systems with different DM masses.
To fulfill this program, we forward modeled the density and pressure distributions of the intracluster medium (ICM), working out the hydrostatic equilibrium equation in fractional gravity. Then, we performed a Bayesian analysis of the X-COP galaxy cluster sample to infer constraints on the fractional gravity parameters for individual clusters and also by stacking them.
We found that fractional gravity performed remarkably well in modeling the ICM profiles for the X-COP sample. We also checked that the relationship between the concentration of the DM profile and the DM mass still remained consistent with the expectations of N-body simulations in the ΛCDM framework. Finally, we confirmed the weakening of the fractional gravity effects toward more massive systems and derived the overall scaling of the fractional gravity parameters from dwarf galaxies to massive clusters, over six orders of magnitude in DM mass. Such an overall trend implies that fractional gravity can substantially alleviate the small-scale issues of the standard DM paradigm, while remaining successful on large cosmological scales.
In future work, we plan to investigate a theoretical explanation for the empirical scaling of fractional gravity parameters with the DM mass. Hints may come from the connection of these parameters with different MONDian and fractional gravity theories, as partly explored by [1,31,34]. In fact, it has been pointed out that all these frameworks are characterized by an index (in our case s) interpolating between the Newtonian and a MOND-like regime, and by a length scale ℓ ∼ √ G M/a 0 that dimensionally can be written in terms of a MOND-like characteristic acceleration scale a 0 and of the system's mass M (baryons in the basic MOND theory, total mass dominated by DM in our case). However, the empirical scaling between ℓ ∝ M ∼2/3 found here is barely consistent with this law within the uncertainties, and an ab initio explanation of the inverse dependence of s with mass M is difficult to be envisaged even in simple terms. This indicates that some crucial ingredient is missing to build a robust theoretical background behind DM in fractional gravity.
On the observational side, the present work clearly shows that whatever the ultimate origin of the fractional gravity behavior in the DM component is, most of its effects manifest in small DM masses; according to the canonical structure formation scenario, these objects must have formed at early cosmic times. Therefore, in the near future, we plan to look for signs of fractional gravity via kinematic (and possibly gravitational lensing) observations of low-mass galaxies at intermediate/high redshift, and of their relics in the local Universe; this could shed light on the mechanisms responsible for the origin and the emergence of fractional gravity across cosmic times.