Limit on graviton mass from galaxy cluster Abell 1689

To date, the only limit on graviton mass using galaxy clusters was obtained by Goldhaber and Nieto in 1974, using the fact that the orbits of galaxy clusters are bound and closed, and extend up to 580 kpc. From positing that only a Newtonian potential gives rise to such stable bound orbits, a limit on the graviton mass $m_g<10^{-29}$ eV was obtained (PRD 9,1119, 1974). Recently, it has been shown that one can obtain closed bound orbits for Yukawa potential (arXiv:1705.02444), thus invalidating the main \emph{ansatz} used in Goldhaber and Nieto to obtain the graviton mass bound. In order to obtain a revised estimate using galaxy clusters, we use dynamical mass models of the Abell 1689 (A1689) galaxy cluster to check their compatibility with a Yukawa gravitational potential. We assume mass models for the gas, dark matter, and galaxies for A1689 from arXiv:1703.10219 and arXiv:1610.01543, who used this cluster to test various alternate gravity theories, which dispense with the need for dark matter. We quantify the deviations in the acceleration profile using these mass models assuming a Yukawa potential and that obtained assuming a Newtonian potential by calculating the $\chi^2$ residuals between the two profiles. Our estimated bound on the graviton mass ($m_g$) is thereby given by, $m_g<1.37 \times 10^{-29}$ eV or in terms of the graviton Compton wavelength of, $\lambda_g>9.1 \times 10^{19}$ km at 90\% confidence level.


I. INTRODUCTION
A century after its inception, General relativity (GR) passes all observational tests at solar system and binary pulsar length scales with flying colors [5][6][7]. The recent direct detection of gravitational waves has confirmed the validity of general relativity in the dynamical strong-field regime [8]. Despite this, a whole slew of modified theories of gravity have been explored, ever since the equations of GR were first written more than a century ago. Most of the recent resurgence in studying and proposing the plethora of modified theories of gravity has been driven by the need to address problems in Cosmology such as Dark matter, Dark Energy, Inflation, and Baryogenesis [9][10][11][12][13][14][15][16][17][18], which cannot be explained using GR and the Standard model of particle physics. Independent of cosmological problems, a number of alternatives have also been extensively proposed to resolve conceptual problems in classical GR at the interface of fundamental Physics, such as Big Bang singularity [19], arrow of time [20,21], or the quantization of gravity [22][23][24]. An updated summary of almost all the modified theories of gravity can be found in monographs such as [25] and also in recent reviews [26][27][28] and references therein.
One such modification to GR, namely massive gravity in which a graviton is endowed with non-zero mass dates back to more than 70 years. The first ever theory of massive gravity in the perturbative limit was proposed by Pauli and Fierz [29]. However, this theory does not reproduce the GR result, in the limit when the graviton mass goes to zero, usually referred to in the literature as the vDVZ discontinuity [30,31]. However, Vainshtein showed that this discontinuity is due to how the * E-mail: shntn05@gmail.com gravitational degrees of freedom are treated during the linearization procedure and can be fixed in a non-linear version of massive gravity [32]. Bouleware and Deser then showed that this non-linear version has a ghost [33]. Therefore, the field of massive gravity theories lay dormant because of these conceptual problems. However in the last decade, the Bouleware-Deser ghost problem has been solved, leading to a resurgence of interest in these massive gravity theories [34][35][36][37]. These massive gravity models can address multiple problems in cosmology such as dark energy [38], dark matter [39,40], inflation [41] and also in fundamental physics related to quantization of gravity [42].
One generic feature of massive gravity models is that the gravitational potential has a Yukawa behavior in the linear weak field limit, typically parameterized as [1,[43][44][45]: where λ g is the Compton wavelength of the graviton and is given by λ g ≡ h mgc , where m g is the graviton mass. Basically there are three model-independent methods, which have been used to obtain graviton mass bounds [45]. The first method involves looking for a weakening of the gravitational force due to a Yukawalike potential. The second type of constraint comes from looking for fifth force interactions, which arise in massive gravity models. The third type of limit comes from the propagation of gravitational waves, either due to modified dispersion relations or from difference in arrival times between gravitational waves and other astrophysical messengers (photons, neutrinos). In the gravitational wave literature, the limits from the first two types of measurements are referred to as "static" bounds, whereas the limits from gravitational wave observations are referred to as "dynamical" bounds. In addition to these three traditional methods, one can also obtain bounds on graviton mass by studying its implications for cosmology, such as large scale structure and late time evolution. But these are strongly model-dependent. In addition, one can also get constraints on Yukawa gravity from the weak-field limit of certain modified gravity theories such as f (R) gravity or Moffat's Scalar-Vector-Tensor gravity [46,47]. A comprehensive summary of all the observational/experimental bounds on the mass of the graviton as well as future prospects can be found in Ref. [45] and a tabular summary can be found in Table 1 of the same paper. We now briefly recap the limits from these different types of methods.
Two years ago there was a watershed event in the history of physics, due to the direct detection of gravitational waves from the two LIGO detectors in Hanford and Livingston [48]. These observations enabled us to obtain the most stringent dynamical bounds on the graviton mass, by looking for dispersion in the observed signal as it propagated towards the detectors. The first detection from GW150914 [48] provided a limit of λ g > 10 13 km, or m g < 10 −22 eV [49], based on looking for a modified dispersion relation for a non-zero graviton mass. Subsequently, a more stringent bound of m g < 7.7 × 10 −23 eV has been obtained using GW170104 [50]. We note however that Deser [51] has pointed out that no strong field generation of radiation in massive gravity models can reproduce the observed ringdown patterns observed in LIGO. For a gravitational wave source in our galaxy, one could also constrain this mass using the line of sight Shapiro delay from the source of the gravitational wave [52]. Most recently, the direct detection of photons in coincidence with the gravitational waves from a binary neutron star merger have enabled us to constrain the mass of the graviton to m g 10 −22 eV [53]. In the future, eLISA could obtain bounds of m g < 10 −26 eV [43], and from the detection of inflationary gravitational waves from stage IV CMB experiments, one could get a bound of m g < 3 × 10 −29 eV [54].
Many massive gravity models give rise to a fifth force. However, these results are theory dependent and in particular depend on how the non-linear Vainshtein mechanism operates in these theories [45]. However, one common feature in these models is the existence of a Galileonlike scalar. As of now, the bounds on graviton mass in this category have been obtained from the decoupling limit of DGP [55] and dRGT [56] theories. Data from lunar laser ranging experiments give a mass bound of m g < 10 −30 eV within the context of the decoupling limit of dRGT theory [45]. From the corresponding decoupling limit of DGP, future surveys on galaxy-galaxy lensing could set a bound upto m g < 10 −33 eV [57].
We finally recap the limits on graviton mass by looking for Yukawa-type fall off of the gravitational force. The first such bound was obtained by Hare [58], by assuming that the gravitational force from the center of the galaxy is diminished by factor of less than 1 e . From this argument, a mass bound of m g < 6.7 × 10 −28 eV was obtained [58]. A similar reasoning was then extended by Goldhaber and Nieto to extragalactic observations of galaxy clusters [1].
The current best limit (from all the three types of methods) on the mass of a graviton comes from the measurements of weak lensing cosmic shear [59], obtained by comparing the variance of the modified shear convergence power spectrum in massive gravity models to the observed data [60]. 1 By imposing the condition that the observed deviations from the ΛCDM power spectrum are less than 1σ, a limit of m g < 6 × 10 −32 eV was obtained [60]. One assumption however made in obtaining this limit is that the graviton mass has no effect on the cosmological expansion, growth of structure and also the CDM transfer function. Furthermore, there is also a degeneracy in the modified power spectrum between a non-zero graviton mass and other cosmological parameters. To evade this degeneracy, the other parameters were determined using the values of the power spectrum at smaller values of the radius, for which the effect of a non-zero graviton mass is assumed to be negligible. These fitted parameters were then used for the limit on graviton mass using the measurements for larger values of the radius [60].
The constraint on graviton mass in Ref. [1] using galaxy clusters was obtained by assuming that the orbits of galaxies in clusters are bound as well as closed and using the fact that the maximum separation between galaxies from the Holmberg galaxy cluster catalog is about 580 kpc [61]. The limit on graviton mass was obtained by positing e −1 ≤ exp(−µ g r), (where µ g is the reciprocal of the reduced Compton wavelength) and assuming r = 580 kpc. This condition implies that there are at most O(1) departures from Newtonian gravity at the edge of the galaxy cluster. The estimated limit on graviton mass thus obtained was µ g < 5.6 × 10 −25 cm −1 or m g < 1.1 × 10 −29 eV or λ g > 10 20 km. We note that this is a very rough estimate. This limit does not use any dynamical mass information for the galaxy cluster or any ansatz for the potentials of the different cluster components (gas, galaxies, dark matter). Also, no confidence interval was given for this upper limit.
Furthermore, very recently it has been shown that Newtonian gravity is not the only central force that gives rise to bound orbits and one can also get bound orbits for potentials which do not satisfy Bertrand's theorem [2,62]. In particular, Mukherjee et al [2] have shown that one can get single-particle bound orbits in a Yukawa potential for certain values. Thus, the main edifice upon which this bound of m g < 10 −29 eV has been obtained [1] is no longer valid.
Galaxy clusters are the most massive gravitationally bound objects in the universe and provide an excellent laboratory for studying a diverse range of topics from galaxy evolution to cosmology (For reviews, see [63][64][65][66]). In the past two decades a large number of galaxy cluster surveys in the optical [67][68][69][70], microwave [71][72][73], and Xray [74][75][76][77][78][79] have come online. These surveys have enabled the discovery of a large number of galaxy clusters up to very high redshifts, allowing us to probe a wide range of questions in astrophysics and cosmology. Multiple observables from these surveys such as galaxy cluster counts, gas mass fraction, and dynamics of galaxies within clusters have been widely used to design tests and constrain a large class of modified theories of gravity, which dispense with dark energy [66,[80][81][82][83][84][85][86][87][88]. Galaxy clusters have also been used to constrain modified gravity theories, which dispense with dark matter, e.g., various incarnations of MOND-like theories, Verlinde's entropic gravity, Moffat's MOG theory, nonlocal gravity [89][90][91][92][93][94][95][96][97]. However, despite the wealth of exquisite multiwavelength galaxy cluster observations, no improvement to the initial estimate on graviton mass using galaxy clusters has been obtained after Goldhaber and Nieto's 1974 paper. 2 The only related result is by De Martino and Laurentis, who obtained a constraint on a variant of the Yukawa gravitational potential considered here (from the post-Newtonian limit of f (R) gravity), using the thermal SZE profile of the Coma cluster from the 2013 Planck observations [46,47]. In principle however, the analysis in this work could be extended to obtain a limit on the graviton mass.
Therefore, to the best of our knowledge, we are not aware of any direct constraint on graviton mass using galaxy clusters from completed stage II or ongoing stage III dark energy experiments, or any forecast on the estimated sensitivity to graviton mass from upcoming stage IV experiments such as LSST [98], Euclid [99], WFIRST [100], etc. This is despite the fact that one of the key science driver for these upcoming surveys is to test modified gravity theories [101].
Therefore, to rectify the situation and to see how sensitive current galaxy cluster data is to graviton mass compared to very rough estimates from four decades ago, we do a first end-to-end study to obtain a limit on graviton mass using the Abell 1689 galaxy cluster, for which exquisite multi-wavelength data is available, allowing the reconstruction of detailed mass models for this cluster in literature.
This manuscript is organized as follows. We discuss the dynamical modeling and mass estimates in Section II. Our analysis and results can be found in Section III. We examine the robustness of our limit to different mass models in Section IV. We conclude in Section V.

II. DYNAMICAL MODELING OF A1689
We use the galaxy cluster Abell 1689 (hereafter A1689) for our analysis. A1689 is one of the largest and most massive galaxy cluster located at a redshift of 0.18. In the past decade, it has been subjected to intensive dynamical modeling within the context of the ΛCDM cosmological paradigm, using multi-wavelength observations from weak and strong lensing, SZE and X-Ray observations [102][103][104][105] (and references therein). These observations have enabled us to obtain estimates separately for the dark matter, gas, and galaxy components for this cluster. Most recently, this cluster was extensively studied to see if its available data is compatible with MONDlike theories, which provide a solution for the dark matter problem from a modification of Newtonian gravity and without the need for dark matter [3,4,106]. We use the same modeling from these papers to test to what extent the data is viable with a Yukawa potential in order to constrain the graviton mass. The first step in this procedure involves estimating the total mass of the different components of the galaxy cluster, viz. its dark matter content, galaxies, and gas in the intra-cluster medium. We follow the same procedure as in Ref. [4].
The total dark matter mass can be obtained by assuming the density obeys the Navarro-Frenk-White profile [107] and is given by [4] : where r s and ρ s represent the dark matter halo scale radius and scale density respectively. They are usually obtained from the relation between the NFW concentration parameter (c 200 ) and the total mass at a radius 200 times the critical universe density (M 200 ). To calculate the total dark matter mass, we use the NFW concentration parameters for this cluster measured by Umetsu et al [104], viz. c 200 = 10.1 ± 0.82, M 200 = (1.32 ± 0.09) × 10 15 M h −1 , obtained using a combination of weak and strong lensing observations. Masses obtained from weak or strong lensing do not depend on the dynamical state of the cluster and hence do not rely on assumptions of hydrostatic equilibrium. We however note that spherical symmetry has been assumed in the dynamical mass modeling, whereas there is observational evidence that this cluster has triaxial symmetry [104].
Although this cluster has been modeled using ellipsoidal halo [104], for this work spherical symmetry has been assumed throughout. The central galaxy (often called BCG, which is an acronym for Brightest Cluster Galaxy) mass distribution is modeled by positing a density distribution of the form [3,102]: where M cg and R cg represent the BCG mass and core radius respectively; R co represents the core size. The values for M cg , R cg , and R co that we use for our analysis can be found in Refs. [3,4,106], which we use for our analysis. The gas mass is obtained using a cored Sersic profile and given by [3,108]: where m p is the proton mass, n e0 is the central electron density; R g represents the radial extent of the gas; while k g and n g are dimensionless parameters which control the shape of the gas profile. The values for all these parameters can be found in Refs. [3,4]. The total baryonic mass M bar upto a given radius R can be found by integrating the galaxy and gas density profiles from Eq. 3 and Eq. 4 : M bar = R 0 4π[ρ gal + ρ gas ]r 2 dr. We note that these mass estimates have been made by positing spherical symmetry.
Once, we have calculated the mass of the different components, the total acceleration assuming only Newtonian Gravity (a newt ) from the center of the galaxy cluster is given by

III. RESULTS
In order to test the viability of Yukawa gravity, the gravitational acceleration can be obtained from the derivative of the Yukawa potential (Eq. 1) and is given by [43]: In the limit that m g → 0, Eq. 6 will asymptote to Eq. 5. In order to test for the validity of this modified acceleration law, we assume that the total mass is the same as that in Newtonian gravity and we only look for deviations compared to ordinary gravity as a function of distance from the cluster center. This is similar to the approaches used to constrain modified theories of gravity, which dispense with the dark matter paradigm [3,4,106].
For this analysis, the NFW density profile, the BCG density profile, and gas density estimates from literature, used to calculate the total mass have been obtained after positing a Newtonian potential. Strictly speaking, the total estimated mass would be larger within the context of Yukawa gravity, because of the weakness of gravity in a Yukawa potential compared to the corresponding Newtonian one. However, a self-consistent constraint on the graviton mass is out of the scope of the current paper and at this time would require significant additional work from the community, since one would need to determine three unknown density profiles (the dark matter, gas, galaxy) in addition to the graviton mass from the observational data. For the dark matter component, one would need to do a suite of N -body simulations in Yukawa gravity as a function of graviton mass and then obtain a parametric estimate as a function of graviton mass. Similarly, the baryonic mass would need to be determined assuming hydrostatic equilibrium in such a modified potential. Therefore, because of the large number of additional free parameters, doing the whole problem self-consistently in order to obtain more robust bounds on graviton mass would pose formidable challenges.
However, if the graviton mass is very small, the deviations in the total mass estimates should be small compared to those obtained using Newtonian gravity and the errors in our estimate of graviton mass should be negligible. Furthermore, since we shall only be interested in deviations in the acceleration profile (compared to a Newtonian potential) for the same mass, the total mass would only be a normalization constant and would not make a difference to the final limit. Therefore, similar to what is usually done in constraining alternate gravity theories, which dispense with the dark matter paradigm (see e.g. [4] and references therein), we assume that the total density profile is the same as in Newtonian gravity and then look for deviations in the acceleration profile as a function of distance from the center of the cluster to constrain departures from a standard Newtonian acceleration profile. Due to the above assumption, the limit is of course conservative. In Sect. IV, we shall examine how the limit on graviton mass changes when varying the mass models for the cluster used here.
We also point out that if we posit that dark matter is made up of massive gravitons (see e.g. [39]), then only the dark matter potential would be modified while the other terms in the potential would be unchanged and our limits would be different. Here, we assume that all the distinct mass components (gas, galaxy, dark matter) uniformly obey the Yukawa potential and dark matter is some hypothetical elementary particle, with the same gravitational laws as the baryonic components.
To quantify the deviations between Newtonian and Yukawa gravity as a function of distance from the center of the cluster, we construct a χ 2 functional given by, where a newt and a yuk are given by Eqs. 5 and 6; σ a is the uncertainty in the estimated acceleration. To get the 90% c.l. upper limit on the mass of the graviton, we find the threshold value of m g for which ∆χ 2 > 2.71 [109], where ∆χ 2 = χ 2 − χ 2 min . We note that χ 2 min =0 corresponds to a zero graviton mass. Therefore, ∆χ 2 is identically equal to χ 2 from Eq. 7. Since, the mass of the graviton cannot be negative, m g = 0 is a physical boundary. Therefore, in such cases ∆χ 2 values for a given confidence interval could in principle get modified compared to the values in Ref. [109]. To obtain the modified ∆χ 2 , we use the procedure recommended by the Particle Data Group [110,111], which has previously been used for neutrino oscillation analysis [112]. The effect of the physical boundary is determined by the difference between the minimum value of χ 2 and the value of χ 2 at the boundary of the physical region. If the minimum value of χ 2 occurs at the physical boundary (which is true in our case), then ∆χ 2 intervals for a given confidence interval are the same as without a physical boundary [111]. Therefore, to obtain the 90% confidence level upper limit we choose ∆χ 2 = 2.71, which is the same as the value without a physical boundary. Alternately, the modified ∆χ 2 threshold can also be obtained using the Feldman-Cousins method [113], which requires extensive Monte-Carlo simulations. However, they have been shown to not differ too much compared to the method use here [111].
We calculated the ∆χ 2 for 24 points between roughly 1 and 3000 kpc, for which errors in acceleration have been estimated from existing observations [3], for which spherical symmetry has been assumed. The first 12 points were located at radii between 3 and 271 kpc, for which the errors in acceleration have been estimated from the line of sight mass density [108], obtained using strong lensing observations [102]. The remaining 12 data points were distributed between 125 kpc and 3 Mpc and the errors in acceleration were estimated from the weak lensing shear profiles [104]. We note that the errors in acceleration data do not include any errors in determination of the radii. χ 2 was then estimated from Eq. 7 for these 24 data points by calculating a newt and a yuk at these radii and using the errors in acceleration estimated in Ref. [3]. This plot is shown as a function of graviton mass in Fig. 1. The 90% c.l. upper limit on the mass of a graviton obtained from ∆χ 2 = 2.71, is given by m g < 1.37 × 10 −29 eV, corresponding to a Compton wavelength of λ g > 9.1 × 10 19 km. For this value of mass, we also show the fractional deviation between the ordinary Newtonian acceleration and that assuming the Yukawa potential in Fig. 2. We can see that for this graviton mass, the differences are less than 1% upto 200 kpc and about 10% at about 1 Mpc.

IV. EFFECT OF SYSTEMATIC ERRORS
We now examine the sensitivity of our limit on graviton mass to different mass models for the three components of A1689, compared to the results in the previous section. A tabular summary of all these upper limits on varying the mass models can be found in Table I. We start with the dark matter part. A large number of groups have obtained different NFW parameters for this cluster using weak and strong lensing data. (See Table 9 10 - 30 10 -29 10 -28 10 -27 m g (eV)    0.83 × 10 15 M , we get m g < 1.43 × 10 −29 eV. Even though, most of the mass in this galaxy cluster is made up of dark matter, we don't expect any major changes with different BCG or gas mass profiles. Nevertheless, to check this, instead of Eq. 3, we used the Hernquist profile [115] (similar to Ref. [94]) for the BCG mass with the same parameters as in Ref. [94]. With this new galaxy mass profile, the new graviton mass limit is the same as before.
Therefore, the change in the limit on the graviton mass by varying our ansatz for the mass models is less than 15% and does not change the ballpark estimate on the limit on graviton mass of m g 10 −29 eV.

V. CONCLUSIONS
In 1974, a limit on graviton mass of m g < 1.1 × 10 −29 eV was obtained from galaxy clusters, using the fact that the orbits of galaxy clusters are bound up to 580 kpc [1] and such closed bound orbits can only exist within Newtonian gravity. However, recently it has been shown that one can get closed bound orbits for a Yukawa potential [2]. Therefore, the main premise used to obtain the mass bound limit from galaxy clusters in Ref. [1] can no longer be justified and this result should no longer be quoted in the literature.
Subsequently, even though a huge amount of work has been done in testing a plethora of modified gravity theories with galaxy clusters using optical, X-ray, and SZE data, we are not aware of any other work on estimating a bound on graviton mass from clusters, despite a wealth of new precise observational data in the past decade, courtesy a whole slew of multi-wavelength surveys.
We obtain a limit on graviton mass from A1689 using an independent method compared to Ref. [1]. We use recent dynamical mass models of the different components of galaxy cluster A1689, obtained using X-ray, weak and strong lensing data [3,4,106] to obtain a limit on the graviton mass. For this purpose, we assume that the potential due to the gas, galaxy, and dark matter all follow a Yukawa behavior, due to non-zero graviton mass. We then look for deviations from the estimated acceleration data (assuming validity of Newtonian gravity) and a Yukawa potential, and find the critical graviton mass for which the ∆χ 2 difference between the two potentials crosses 2.71. This gives us a 90% c.l. upper bound on the graviton mass of m g < 1.37 × 10 −29 eV or on the Compton wavelength λ g > 9.1×10 19 km. We also checked how the limit varies with different mass models for the dark matter and BCG potential. We find that the maximum variation in the limit on graviton mass is about 15% and thus does not change the ballpark estimate of our limit, which is O(10 −29 ) eV.
We should point out that the fact that our upper limit is approximately of the same order of magnitude as that obtained by Goldhaber and Nieto [1] is only a coincidence. The maximum size they assumed for the galaxy cluster orbits is about 580 kpc, as this was the size of the largest known clusters in 1974. Using this estimate for the size, they obtained an upper limit of O(10 −29 ) eV, which is of the same order of magnitude as ours. In principle, one could trivially apply the same method [1] to some of the galaxy superclusters currently known. For example, the recently discovered Saraswati supercluster [117] (whose size is at least 200 Mpc) would yield a more stringent upper limit on the graviton mass of about 3 × 10 −32 eV. However, as mentioned earlier the underlying assumptions behind this argument used to obtain the limit are incorrect.
We however note that to obtain our limit, mass estimates for the different components (dark matter, gas, galaxy) have been obtained assuming Newtonian gravity, since otherwise the whole problem of simultaneously determining the mass of the three unknown components in addition to the graviton mass becomes currently intractable, given the large number of free parameters. However, in this case since the limit has been obtained by using deviations from Newtonian acceleration profile as a function of the distance from the galaxy cluster center and the total mass of each component would mainly act as a normalization constant and not make a big difference to our final limit.
Given the large number of upcoming Stage IV dark energy experiments such as LSST [98], Euclid [99], WFIRST [100] etc, it would be interesting to estimate the expected improvement in the limit on graviton mass compared to the result obtained in this work. We now carry out an order of magnitude estimate of the same.
Some of our errors in acceleration data (at radii less than about 150 kpc) come from strong lensing measurements [102], which use Hubble Space Telescope data. Therefore, we do not expect significant improvements in the strong lensing based error estimates. The acceleration errors at higher radii are obtained from weak lensing measurements using Suprime-Cam data [104]. For Euclid and other stage IV experiments, the multiplicative bias from the shear must be less than 0.1% [99,118].
If we evaluate the acceleration errors from these predicted shear errors for distances greater than 150 kpc and combine it with current errors from strong lensing for smaller radii, we expect a 90% confidence upper limit of m g < 2.75 × 10 −30 eV. This is still not as sensitive as the current best limit on graviton mass [60]. However, we caution that this is only a ballpark estimate of expected improvement sensitivity. More detailed forecasting studies need to be done by the relevant working groups from the various stage IV dark energy experiments. One key missing ingredient needed for that purpose is the generation of N -body simulations in Yukawa gravity and the calculation of the corresponding halo mass functions.