UvA-DARE (Digital Academic Repository) Accurate modeling of the strong and weak lensing profiles for the galaxy clusters Abell 1689 and 1835

An accurate, spherically symmetric description of the mass distribution is presented for two quite virialized galaxy clusters, Abell 1689 and Abell 1835. A suitable regularization of the small eigenvalues of the covariance matrices is introduced. A stretched exponential proﬁle is assumed for the brightest cluster galaxy. A similar stretched exponential for the dark matter and halo galaxies combined, functions well, as do thermal fermions for the dark matter and a standard proﬁle for the halo galaxies. To discriminate between them, sensitive tests have been identiﬁed and applied. A deﬁnite verdict can follow from sharp data near the cluster centers and beyond 1 Mpc


Introduction
In the present era of precision cosmology, various cosmological parameters are known at the percent level, while serious tension remains, in particular, concerning the value of the Hubble parameter [1].As a next step, a similar precision is desired for galaxy clusters, to be called clusters from now on.Cluster masses can be several times 10 15 M and their size, several Mpc.A good understanding of suitable clusters may provide an additional grip on the nature of dark matter.
While interesting details can be derived from dynamical clusters such as the Bullet Cluster [2] and the Cosmic Train Wreck Cluster Abell 520 [3,4], their structure is complicated and their analysis subject to questions.We shall focus here on clusters that are reasonably spherically symmetric, so that spherical symmetry is a good, and often employed, approximation.
Apart from interest for its own right, study of clusters [5,6] provides information on the dark matter versus modified Newton force dispute.In particular, the MOND theory [7] has achieved success for galaxies [8].However, for fat clusters with M 200 ∼ 10 15 M , the gravitational acceleration starts out large in the center, and the Newton regime holds up to the MOND radius GM 200 /a 0 ≈ 1 Mpc, given that a 0 ≈ 1.2 × 10 −10 m/s 2 [9].In fact, self-gravitating isothermal spheres are unstable in Newtonian dynamics, hence they would expand to fill up their MOND radius, causing even highacceleration systems like galaxy clusters to be affected by MOND in the sense that their size should correspond a e-mail: t.m.nieuwenhuizen@uva.nl(corresponding author) to their MOND radius [10].Likewise, the observed velocity dispersion profile of Dragonfly 44 falsifies MOG at 5.5σ [11].In short, modifications of gravity, such as MOND, but also EG [12,13], MOG [14] and f (R) are under severe stress [15][16][17].Like the related f (T ) and f (R, T ) theories, they do not change matters appreciably inside the huge 1 Mpc domain, so that Zwicky's conundrum-there must be dark matter or something alike-remains unsolved [16].
In principle, dark matter and modified gravity may co-exist, and this combination may actually be required to fully explain the properties of galaxy clusters like El Gordo [18].In a MOND context, the best developed such proposal is the νHDM framework [19], in which galaxy clusters are explained together with the CMB anisotropies using 11 eV/c 2 sterile neutrinos with the same overall density as the CDM in ΛCDM.This framework might also account for the Hubble tension [20].
Another issue of importance is to establish whether dark matter is self-interacting.Analysis of clusters puts forward a cross section-to-mass ratio σ/m ∼ 1 cm 2 /gr [21,22], although the question is not settled [23].This large value excludes a lot of parameter space for various models of dark matter.Indeed, for MACHOs of Earth mass, the cross section would be comparable to the size of the Earth orbit (π AU 2 ), while in reality its cross section is πR 2  ⊕ with R ⊕ the Earth radius.Hence any type of 100% MACHO dark matter, even if consisting of axion stars or of (primordial) black holes, would be ruled out.Even dark matter particles heavier than, say, 0.4 GeV/c 2 , need mediators to establish the selfinteraction [24].
Clusters are mostly dynamical, meaning that they are an aggregate of subclusters, which still need giga years to get into an equilibrium to form a (meta)stable cluster.For such situations, at best a reasonable description can be achieved.Two well studied exceptions are the clusters Abell 1689 and Abell 1835.A1689 actually has one in-falling subcluster, far away from the center, with not very large mass, so that including or excluding the quadrant in which it lies, does not cause marked differences.A proper description of A1689 requires triaxiality [25], but actually, the mostly employed spherical approximation functions rather well.A1835 looks even more symmetrical visually, and the spherical approximation works well, but see [26] for triaxiality in its X-ray gas and lensing.
The setup of this article is as follows.In Sect. 2 we present lensing and gas data for A1689 and A1835.In Sect. 3 we discuss models for their mass distributions.In Sect. 4 the fits to A1689 are discussed and in Sect. 5 those to A1835.In Sect.6 we compare these fits and we close with a summary.

Data for Abell 1689 and Abell 1835 2.1 Strong lensing data for A1689
The galaxy cluster Abell 1689 lies at redshift z = 0.183 and acts as a strong lens, lensing many background galaxies into a number of up to 5 arclets, i.e., pieces of the ideal Einstein ring.From the observed strong lensing (SL) arclets, 2D mass maps have been generated using the program Lenstool [27], a strong lensing inversion algorithm.The code utilises the positions, magnitudes, shapes, multiplicity and spectroscopic redshifts for the multiply imaged background galaxies to derive the detailed mass distribution of the cluster.The overall mass distribution in cluster lenses is modeled in Lenstool as a superposition of smoother largescale potentials and small scale substructure that is associated with the locations of bright, cluster member galaxies.Both potentials are described using parametric mass models.The parameters are adjusted in a Bayesian way, i.e., their posterior probability is probed with a MCMC sampler.This process allows an easy and reliable estimate of the errors on derived quantities such as the amplification maps and the mass maps.
This inversion is an underdetermined problem, so that an ensemble N = 1001 of maps compatible with the data is produced.Integrated over the interior of circles around the cluster center, this yields data for M 2d (r), the mass inside a cylinder of projected radius r around the sightline to the cluster centre [17].For each map m, these M (m) 2d values are evaluated at radii r n = r 1 a n−1 with n = 1, . . ., 149 and a = 1.0388, such that (r 1 , r 149 ) = (3.15,879) kpc.Only N = 117 of these r n contain physical data, the other ones are omited.The ensemble averages M 2d (r n ) define the cylindrical mass density Σ n = M 2d (r n )/πr 2 n , while their covariances Γ mn also follow as averages over the maps [17].The Σ n data with error bars equal to (Γ nn ) 1/2 are pre-Fig.1 Strong lensing data for the cylindrical mass density Σ as function of the radius for the clusters A1689 (upper, black) and A1835 (lower, blue).Both clusters behave similarly around their centers; A1689 has more mass around 100 kpc; the clusters behave similarly beyond 500 kpc sented in Fig 1, for the present cluster A1689 and for the later discussed cluster A1835.For small and large r, they show quite similar behavior, while at intermediate r, A1689 is denser than A1835.

Regularization of the covariance matrix
We shall fit theoretical models for Σ(r) by minimizing In principle, one has C ij = Γ ij .However, the matrix Γ has a big spread of eigenvalues, from 4.4×10 −15 to 0.216 gr 2 /cm 4 .The near-degeneracies arise since the 2d maps are based on bins of which many are empty.The small eigenvalues are somewhat reflected in the small diagonal element around 120 kpc, see Fig. 1.But Eq. ( 1) with C = Γ is dominated by the very small eigenvalues, which are numerical artefacts.To regularize them, it is customary to employ a Tikhonov regularization counting for further scatter, by adding a constant to the diagonal elements of Γ [28][29][30], In [29], where the data are binned, we take σ SL = 0.19 gr/cm 2 and in [30] 0.16 gr/cm 2 ; the latter value is employed in Fig. 2. It is seen that the Tikhonovregularized C ii lie for large r much above the empirical values Γ ii , so that those data hardly play any role in the analysis.. To acknowledge the decay of Σ as function of r, we add, instead instead of (2), a constant times Σ 2 to the diagonal, so that we adopt instead the "poor man's regularization" it is seen to generalize Eq. ( 2) for cases where Σ changes appreciably, as happens in cluster lensing.This regularization weighs the contributions to χ 2 with basically equal weight for all r i .In Fig. 2 the lower (black) data present the Γ ii .The red points (upper on the right side), present Eq.(2), while the blue ones (middle on the right), presenting Eq. ( 3), do more justice to the Γ ii data.

The transversal shear
The transversal shear is defined as where Σ was introduced above, Σ c is a constant, and Σ is the line-of-sight density (projected density, 2d density) at transversal distance r of the center, Weak lensing (WL) data for Σ and for g t with Σ c = 0.6815 gr/cm 2 in A1689 have been reported by Umetsu et al. [31] and employed by us [30].They are represented in Figs. 5 and 8, respectively.Since here all bins were filled, their covariance matrices are well behaved.

Mass profile of the X-ray gas
The mass density of the X-ray gas follows from the observation of the electron density.Recent data at r < 900 kpc have been presented in [17], which fit well to a cored Sérsic profile [30] n Data for r > 1 Mpc, taken from Planck-ROSAT [32], fit to These behaviors combine into the global fit so that n e (0) = n 0 e .The best fit parameters are For a typical Z = 0.3 in units of Solar metallicity [33,34], the mass density of the X-ray gas reads where α ≈ 15/13 is the average number of nucleons per proton, m H is the mass of a neutral hydrogen atom, and n e (r) is the electron number density.

Generating Σ data by numerical differentiation
It is useful to derive data for Σ, which can be obtained from the SL data for Σ.We start from the relation which follows from the relation between the cylindrical mass density Σ and the line-of-sight mass density Σ, where M 2d is the mass in a cylinder of radius r around the cluster center.By taking dΣ and d log r from 58 pairs of adjacent data points, we work with the numerical derivative to be considered at the geometrical average position, Here Σ (1) n and r (1) n denote the original unbinned data (having bin size b = 1).The data for Σ follow from these as In these definitions, the superscript denotes that two Σ data points are employed for each Σ value.The covariance matrix of the Σ, to be denoted as Γ(Σ (2) ) follows by the rule ( 16) from the covariance matrix Γ(Σ) ≡ Γ.It has to be regularized as well.In analogy with (3) we take It turns out that several of the Σ (2) n are negative, an unphysical effect arising from noise in the data and/or lack of perfect sphericality.This can be overcome when first binning the data in bins of b = 3 points.For general b one has The binned position is located at After the binning, the quantities are defined as in (16), and located at the binned position With the same rules taken bilinearly, the correlator Γ(Σ (2b) ) follows the Γ(Σ), and its regularization and the common value α = √ 0.001 = 0.032.We bin the data in bins of b = 3; this produces 19 pairs of data points from which Σ is determined.Likewise, we produce the related 19 data points for the combinations Σ − Σ and 2Σ − Σ.
Their covariance matrices are constructed along similar lines, with a common value of α.Like the Σ, these observables are made up by combining two adjacent bins, with coefficients determined by ( 18) and (20).Their covariances, inherit small eigenvalues from Γ(Σ), so regularization is again needed.We take It is our philosophy to use all information available, and in particular, not to skip the data for, say, r < 300 kpc.While such a restricted data set would involve small scatter also in the differentiated data, its fit would be stronger driven by the regularization, see Figs. 2  and 3.The global outcome of such analyses is that the present fits still work, though less decisive, while other fits would be less ruled out, or even comparably acceptable.

Data for A1835
The galaxy cluster Abell 1835 lies at redshift 0.2532 and shares many characteristics with Abell 1689.In Fig. 1 it is seen that it has a similar mass in the center, less mass around 100 kpc, and quite similar mass beyond 500 kpc.
Strong lensing and X-ray data for A1835 were presented recently by us [17].Our routine generates mass maps as in A1689, at the radii r n = r 1 a n−1 for n = 1, • • • , 149.This involves the same ratio a but different r 1 , such that with (r 1 , r 149 ) = (4.027,1120) kpc.Again 117 points contain physical information.We produced 1001 mass maps M (m) 2d , which are averaged to yield data for Σ.
In A1835 the covariance matrix for Σ encounters the same problem as in A1689: there are very small eigenvalues (they range from 4.4 ×10 −15 to 0.22 gr 2 /cm 4 ) and the diagonal elements have a minimum at 187 kpc, see Fig. 3. Instead of the Tikhonov regularization (2), we again adopt the "poor man's regularization" (3), here with α = 0.05.The elements Γ ii and C ii are presented in Fig. 3.The values of C ii /Σ 2 i vary over a factor ii 2 4 4 Fig. 3 The empirical covariance elements Γii as function of the ri (lower data, black).The Tikhonov regularization (2) with σSL = 0.03 gr/cm 2 (red) exceeds the data strongly at large r.The "poor man's regularization" (3) with αSL = 0.01 (blue) does more justice to the large-r data 9.3.Making α still smaller would enhance this ratio and induce an overfitting of the data around 187 kpc.
As before, we also construct 19 data points for Σ, Σ − Σ and 2Σ − Σ from the Σ data in A1835.
The X-ray data for the electron density have been presented by us as well [29].We found that the following profile explains the data well, The best fit parameters are n 0 e = 0.0927 ± 0.0070 cm −3 , R 0 = 91 ± 13 kpc, R 1 = 31.8± 2.9 kpc, R 2 = 169 ± 15 kpc.
In these clusters, the gas mass density only becomes significant at large r, because the galaxies are dominant at small r.
3 Theoretical mass profiles

Generalities
A celebrated profile in astrophysics is the Sérsic profile for the line-of-sight (2d) luminosity of galaxies, which has the form of a stretched exponential, see Eq. ( 26) below.In cosmology, the most popular mass profile is the so-called NFW profile, (see Eq. ( 30) below), which has a cusp at the origin and falls of as a power law [35].
With further coauthors, the NFW authors observe in dark-matter-only simulations that the stretched exponential profile describes the 3d mass density better than the NFW profile [36].As exposed in their Fig. 6, the authors find good fits for families of dwarf galaxies, families of galaxies and families of galaxy clus-ters, with mass densities ranging from ×10 4 M /kpc 3 to 10 8 M /kpc 3 .The scales range from 0.2 kpc for dwarf galaxies, to several hundred kpc for galaxies, up to 1.5 Mpc for fat galaxy clusters, respectively.The mean Sérsic n values are 3.0 for dwarf-and galaxy-sized halos and 2.4 for cluster-sized halos, similar to the values that characterize luminous elliptical galaxies [36].We consider two mass profiles: the double stretched exponential profile (DSE) and thermal fermions.

The stretched exponential BCG mass profile
A stretched exponential profile has three parameters, the amplitude ρ 0 , the scale R, and the index n, It has total mass where Γ(n) is the standard generalization of (n − 1)!It corresponds to a central line-of-sight mass density

Double stretched exponential profile
The stretched exponential is an interesting candidate to model the combined mass density of the dark matter and the galaxies.For this aim, one assumes a sort of equilibration between them.To put it bluntly, for this profile one works with the adagio "where there are stars, there can not be dark matter".Since the central, brightest cluster galaxy is much heavier than the cluster halo extrapolated towards the center, we adopt an additional stretched exponential for it and arrive at the double stretched exponential profile (DSE).For the modeling of clusters, we thus assume a stretched exponentials for the total cluster (halo, h), and an additional one for the brightest cluster galaxy (BCG, b).Incorporating the gas, the total mass density reads In the central regions, ρ g is much smaller than the other two, and hence irrelevant.While ρ h decays exponentially fast to zero for r beyond R h , ρ g only decays as a power law, so it assures a power law decay of the total density.

NFW profiles for the halo
A popular profile is the NFW profile, and its generalization the gNFW profile (any n), This often employed profile has first been inferred from dark matter-only simulations, and it is often supposed to hold with the baryon density is included.This has the benefit of very few fit parameters (2 for NFW, 3 for gNFW).A drawback is then that it does not provide a handle on the matter density of the galaxies.
Putting things together, we have a stretched exponential for the BCG, a (g)NFW for the halo, on top of the gas density, viz.

Thermal fermionic dark matter
In 2009 we found the first indications that thermal fermions provide a good fit for lensing data of the cluster Abell 1689 [34].Our followup studies have confirmed this [17,29,30].Here we subject this profile to a more stringent test.
Consider g thermal fermion species of mass m at temperature T = mσ 2 , where σ is the velocity dispersion and μ chemical potential per unit mass, in the local gravitational potential ϕ(r).Their mass density reads The index ν expresses that a possible realization lies in sterile neutrinos, as suggested by [19,20] in the context of MOND [7], for which also arguments from the El Gordo cluster were given [37].The logarithmic integral is in general defined as with Γ(α) Euler's Gamma function and Li α (y) the standard logarithmic integral.For (α) > 0 the integral is well defined at all real x.The sum converges for x ≤ 0, so that Li α (x) ≈ e x for x −1.A better approximation is Exact to order e 6x is the Padé approximant For any α > 0 the coefficients are positive.In our case α = 3 2 , Eq. ( 36) takes at x = 0 the value 0.768095, close to the exact value (1 − 1 2 1/2 )ζ 3/2 = 0.765147.In this approach, the baryonic mass has to be specified.Various components contribute to the baryonic matter: the brightest cluster galaxy, the other ("halo") galaxies, globular clusters, cold gas clouds, the X-ray gas, etc.A fit to data for the X-ray gas has been discussed above.For the brightest cluster galaxy (BCG, "b") we adopt the previous stretched exponential form All other parts are lumped into the term "mass density of galaxies" (G).An adequate profile with total mass M G , inner scale R c and outer scale R g is [38] These components model the total mass density of galaxies ρ b (r) + ρ G (r) which has at r = 0 the property The gravitational potential ϕ, which enters ρ ν in Eq. ( 33), is solved self-consistently from the Poisson equation (40)

Lensing observables
We focus on the line-of-sight mass density (2d-density) and the cylindrical mass density In the fermionic application, the Poisson equation allows to express Σ(r) as [29] Σ and Σ(r) as the simpler expression [34] Σ We also consider the combinations Σ − Σ and 2Σ − Σ.If the mass density is cored at the origin, Σ − Σ will vanish there, so this combination tests the central behaviors.2Σ − Σ, on the other hand, tests the decay at large r.Indeed, consider an isothermal fall off ρ ≈ σ 2 /2πGr 2 , for which In our cases where always exceeds at intermediate r its large-r asymptote, the "excess" mass M 0 is positive.While Σ − Σ starts linearly from 0 at the origin, it will decay as σ 2 /2Gr + M 0 /r 2 for large r.On the contrary, the combination 2Σ−Σ starts at some finite Σ(0), while at large r the leading terms cancel, leaving a −M 0 /r 2 decay.Obviously, it changes sign at some finite r; hence 2Σ − Σ is a sensitive quantity for testing the large-r regime.It actually holds by definition that so its zero crossing occurs at the maximum of M 2d /r.This is reminiscent of the circular rotation speed of objects in the cluster, v rot (r) = GM 3d (r)/r.The circular speed can have a maximum and, with it, M 3d /r.Eq. ( 46) changes sign at an intermediate r.

Fits for A1689
We fit models for the mass distribution in A1689 to the strong lensing data for Σ and the weak lensing for Σ and g t data.We combine with the SL data for Σ, Σ − Σ and 2Σ − Σ, derived numerically from Σ, while neglecting their mutual correlations.Hereto one may imagine that each of them derives from averages over 250 of the in total 1001 M (m) 2d maps.Alternatively, one may view the derived values of Σ, Σ − Σ and 2Σ − Σ just as tools to optimize fit to the Σ data.
In A1689 the total χ 2 is taken as The first term involves the 117 SL data for Σ, with a weight factor 1/4 adopted to compensate for notbinning these data.χ 2 W L (Σ) involves the 14 WL data points for Σ and their covariance matrix, and χ 2 W L (g t ) involves the 13 WL data points for g t and their diagonal covariance matrix.The correlation matrix C(Σ) is regularized by adding α 2 Σ 2 on the diagonal of Γ(Σ).Likewise, for C(Σ − Σ) and C(2Σ − Σ) we add α 2 (Σ − Σ) 2 and α 2 (2Σ − Σ) 2 , respectively, to their diagonals.In A1689 we adopt the values We neglect the correlation between the various SL terms in (47), but keep the off-diagonal elements in each one.
We have attempted various further regularization schemes, without much improvement of the fits.
As typical when working with SL data that involve very small eigenvalues, the choice of our regularization parameter α in Eq. (3) needs some care.Extreme cases are to be avoided: a too large α effectively discards all information in the covariance matrix, while taking it too small gives too much weigth to the numerical artifacts in it.Hence it has been selected to get values χ 2 /ν of order 1 for the best fit.It then serves to establish the relative quality of fits, and varying it within a reasonable range does not alter the relative quality much.

Double stretched exponentials
The minimum of χ 2 is determined as function of the parameters.The errors in the parameters {p i } are set by linear regression.First, the Hessian H ij = 1 2 ∂ pi ∂ pj χ 2 is calculated by numerical differentiation and inverted.The diagonal elements provide the 1σ error bars δp i = (H −1 ) ii , while the off-diagonal elements represent parameter covariances.The best fit in A1689 reads 123 While the halo is quite well confined, the brightest cluster galaxy is less so; this is no big surprise, since the central data are scarce and relatively noisy.The M b value is to be compared with (13.0 ± 2.7)10 11 M from Loubser et al [39].The values of the separate χ 2 terms of Eq. ( 47) are listed in the second column of Table 1.The last line gives χ 2 /ν, with the number of free parameters ν = N − 6 for the DSE model.The number of data points corresponding to (47) is N = 152

NFW-type profiles in A1689
Fitting the 2-parameter NFW profile combined with the gas profile to secure a proper fall off at large r, to the SL data of Σ alone yields the good fit χ 2 /ν = 0.88.However, this deteriorates when other data are included, the reason being in particular the behavior at r ∼ 2 − 3 Mpc.With the weak lensing data for Σ and g t added, the fit is already bad, with χ 2 /ν ∼ 10.When the Σ − Σ and 2Σ − Σ are included, the situation worsens considerably; then the failure at large r already becomes relevant for r < 1 Mpc.To employ only 2 fit parameters is too poor, given the precise data.
Taking the sum of 2 NFW profiles does not work well either.Neither does the case regularly adopted in literature of one NFW profile at small r and another at larger r (which corresponds to taking their maximum).
One of our further attempts to improve the fit involves a gNFW profile for the halo and a stretched exponential for the BCG, to which the gas density is again added for the behavior at large r.The best case is found when the gNFW index is n = 0, that is to say: no cusp in the halo part, the BCG being fully accounted for by the stretched exponential.This cored case n = 0 goes against the philosophy of a cuspy NFW profile.
Since the remaining fit is by far not so good as in the DSE and thermal fermion model, we refrain from presenting further details.

Fits for Abell 1835
In this cluster we take half of the regularization parameters (48) in A1689, This sharpening is permissible, since the cluster looks more regular, and, perhaps, because there are no WL data upto 2-3 Mpc that would put further constraints.

Double stretched exponentials in A1835
We repeat the above procedure.

Thermal fermions in A1835
For the thermal fermion profile in this cluster it appears that the scale parameters R c and R g of the galaxy distribution are both about 100 kpc.Taking them equal, R g = R c , eliminates one free parameter, and giving as (58)

NFW-type profiles in A1835
The NFW situation is comparable to that in A1689, but here no WL data exist, hence no lensing data beyond 1.12 Mpc, which allows more flexibility.The Σ-only fit has again a very good χ 2 /ν = 0.72.Similarly to the case in A1689, this deteriorates when more data are included, the basic reason seeming to be that the NFW profile decays too slowly at large r.

Summary
We have considered precise strong lensing data for the clusters A1689 and A1835.For A1689 we include exist-ing weak lensing data.In both cases, the X-ray gas density is known from observations and fit to an analytical profile.
The strong lensing data have been gathered for cylindrical mass M 2d (r) or, equivalently, the cylindrical mass density Σ(r) = M 2d /πr 2 , where r is the projected distance to the cluster center.After binning, the data allow numerical differentiation, to yield data for the line-ofsight mass density Σ and, consequently, for the combinations Σ−Σ and 2Σ−Σ.They emphasize the behavior at the origin and at large r, respectively.
Fits to a double stretched exponential (DSE) profile and to thermal fermion profile are considered.It is observed that both profiles fit reasonably well, with the fermions fitting better in A1689 and the DSE better in A1835.Somewhat surprisingly, NFW-type profiles fit considerably less well at this level of accuracy, even though the NFW profile is expected in the ΛCDM framework [35]; this is compensated, however, by the good fit of the DSE profile, of which the halo part was put forward as a better profile by the same authors with their collaborators [36].Sharp data beyond 1 Mpc may help to discriminate more between the profiles, with the "winner" expected to be the double exponential profile.NFW and NFW-type profiles fit these precise data considerably less well.
The covariance matrix of the strong lensing data has very small eigenvalues, hence a cutoff is needed, as is well known.We propose a "poor man's" regularization, which is better suited in situations such as lensing, where the observable decays at large distance.The regularization parameter has been chosen to get a reasonable value for χ 2 /ν, but a theoretical criterion to fix it, such as by some maximal entropy condition, would be welcome.
Acknowledgements It is a pleasure to thank the anonymous referee for various suggestions, in particular to put the study better in the MOND-like framework.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Fig. 6
Fig. 6 Data of Σ − Σ in A1689 with the best double stretched exponential fit (red) and the best fermion fit (blue)

Fig. 7 Fig. 8
Fig. 7 Data of r(2Σ − Σ) in A1689 with the best DSE fit (red, lower) and the best fermion fit (blue, upper).Additional data beyond 1 Mpc may settle their difference statistically

Fig. 9 Fig. 10 Fig. 11
Fig. 9 Data of Σ in A1835 with the best Double Stretched Exponential (red) and best fermion fit (blue) fits

Fig. 12
Fig.12 Data of r(2Σ − Σ) in A1835 with the best Double Stretched Exponential fit (red) and the best fermion fit (blue).While the difference between the two models is statistically relevant, precise data beyond 1 Mpc could be decisive

Table 1
The separate χ 2 values in A1689, for the double stretched exponential (DSE) model, and for the thermal fermion dark matter model 1/4m, which we employ for comparison with earlier work.The first six compare well with earlier fits, while the last three, referring to the BCG, are new.The value for M b is adopted from Loubser et al[?].The χ 2 values for the various components are presented in the right column of Table1.The overall value χ 2 /ν = 1.32 represents a good fit.The enclosed masses at r 200 and r 500 in the fermion model are very close to the ones (50) of the DSE model,