# MG-MAMPOSSt, a code to test gravity at galaxy-cluster scales: a technical
introduction

Developed by L. Pizzuti  Osservatorio Astronomico della Regione Autonoma Valle
d’Aosta, Loc. Lignan 39, I-11020, Nus, Italy I. D. Saltas  CEICO, Institute of
Physics of the Czech Academy of Sciences, Na Slovance 2, 182 21 Praha 8,
Czechia  A. Biviano  INAF-Osservatorio Astronomico di Trieste, via Tiepolo 11,
34143 Trieste, Italy  IFPU - Institute for Fundamental Physics of the
Universe, via Beirut 2, 34014 Trieste, Italy  G. A. Mamon  Institut
d’Astrophysique de Paris (UMR 7095: CNRS & Sorbonne Université), 98 bis Bd
Arago, F-75014 Paris, France  L. Amendola  Institute of Theoretical Physics,
Philosophenweg 16, Heidelberg University, 69120, Heidelberg, Germany

###### Abstract

The MG-MAMPOSSt code is a license-free Fortran95 code to perform tests of
General Relativity (GR) through the analyses of kinematical data of galaxy
clusters based on the Jeans’ equation. The code is based on the MAMPOSSt
method, and extends the original code through new parametrisations of the
gravitational potential for general families of gravity theories beyond GR
aimed to explain dark energy. MG-MAMPOSSt is further supplemented with a new
capability to produce weak lensing forecasts for joint kinematic+lensing
analysis.

This document provides a technical description of the code’s new features,
functionality with respect to the original version, and instructions on its
installation and use. Finally, we explain how the code could be further
modified to include a wider family of gravity models and/or density profiles,
that could allow its application in broader theoretical frameworks as well as
other physical systems such as stellar clusters. A detailed forecast analysis
for the modified gravity models currently implemented in the code can be found
in the paper of Pizzuti et al., 2021.

##  1 Introduction

In view of the growing interest amongst the communities of cosmology and
astrophysics to test the nature of gravity and dark energy theories at large
scales, this manual is intended to provide a guide for the functionality and
the basic usage of MG-MAMPOSSt, a Fortran95 code capable of performing tests
of gravity models with the kinematics of member galaxies in clusters. The code
is based upon the original MAMPOSSt method developed by G. Mamon, A. Biviano
and G. Boué (Ref. [1], hereafter MAM13) 111A public version of MAMPOSSt by G.
Mamon can be found at https://gitlab.com/gmamon/MAMPOSSt..But whereas the
original MAMPOSSt code relies on the assumption of a standard Newtonian
gravitational potential, MG-MAMPOSSt implements general and viable models of
gravity beyond General Relativity (GR), with the aim of placing constraints on
their theory space at galaxy-cluster scales, as well as of investigating the
essential statistical degeneracy between model parameters. In addition, the
code is capable of producing complementary weak-lensing forecasts towards
joint kinematics+lensing analyses.222The original MAMPOSSt has been used in a
joint kinematics+lensing analysis by Ref. [2].

MAMPOSSt (Modelling Anisotropy and Mass Profile of Spherical Observed Systems)
determines mass profiles of galaxy clusters (or in general, spherical systems
in dynamical equilibrium) by analysing the internal kinematics of the cluster
members.333MAMPOSSt has also been used for elliptical galaxies traced by
globular clusters and dwarf spheroidals traced by their stars [3], and has
been recently extended into MAMPOSSt-PM to handle proper motions in star
clusters (Ref. 4, see Ref. 5). Given an input of projected positions and line-
of-sight (l.o.s) velocities of the member galaxies, and under the assumptions
of spherical symmetry and dynamical relaxation, the code solves the Jeans
equation to reconstruct the gravitational potential, the velocity anisotropy
profile, and (optionally) the projected number density profile.

MG-MAMPOSSt extends the method to gravity scenarios beyond GR, where the
gravitational potential is explicitly modified by the presence of an
additional scalar degree of freedom, resulting in effective mass profiles
different from GR. The resulting effective mass profile can be confronted
against real or synthetic data provided as input to the code. In addition, the
code requires a parametric modelling for the profiles of velocity anisotropy,
mass and number density, with several available choices which can be tuned as
input. The code’s main output is a tabulated likelihood as a function of all
the free model parameters of gravity and other input physics.

We emphasise that the goal of this document is to deliver the necessary
instructions for the operation of the new features in MG-MAMPOSSt. We will
defer from discussing the details of the original MAMPOSSt method on which the
code is built upon, and we refer to Ref. MAM13 for an exposition of MAMPOSSt.
What is more, for a detailed discussion of the application of MG-MAMPOSSt in
theories beyond GR with synthetic data we refer to Ref. [6].

The current licence-free version has been developed by L. Pizzuti based on the
original Fortran77 code of A. Biviano. I.D. Saltas, L. Amendola contributed to
the theoretical background and equations as implemented in the code, assisted
with the implementation of the statistical analysis in the code and advised in
various other parts. G. Mamon contributed to the development of the weak
lensing forecasts and to the refining of the user interfaces.

This document is organized as follows: In Section 2 we explain the basics
around running the code, and in Section 3 we provide a brief overview of the
method on which MG-MAMPOSSt and the original MAMPOSSt codes are based. In
Section 4 we describe the new modified gravity models implemented in the code.
In Sections 5 and 6 we list all the input parameters that can be defined in
MG-MAMPOSSt. In Section 7 we briefly discuss the output files of the code.
Finally, Section 8 provides an illustrative example of a complete run, and an
outlook is provided in Section 9.

##  2 Preliminaries

Running the code requires a Fortran95 compiler or higher. The code can be
downloaded from the following link

https://github.com/Pizzuti92/MG-MAMPOSSt.

The downloaded folder contains the main source code (gomamposstoptS.f) and all
the support files. The configuration and installation requires at least CMake
3.17.1. In the main MG-MAMPOSSt folder execute the following commands:

mkdir build
cd build/
cmake ..
cmake --build .
sudo cmake --install .
cd ..

To run and test the code, type:

gomamposst < gomamposst_x.inp

which produces the main outputs, or

./scripts/script_runmam.sh --options
which also generates additional plots if the MCMC mode is selected (see
below). Note that to run the above script, permissions should be changed to
make it executable. Otherwise, one can simply use the sh command. Installing
the code and running ./scripts/script_runmam.sh without changing the input
parameters and the operating options should produce the last example provided
in Section 8. The additional options are -d <directorypath>, which select the
directory where the file gomamposst_x.inp (see below) is located - default is
the current working directory - and -t which prints the execution time of the
code.

Input data should be stored in a text file named datphys.dat in the folder
data. The file is structured as a table where the number of rows coincides
with the number of data points. The first column is the projected radius in
units of kpc, the second and thirds columns represent the l.o.s. velocities
and the associated errors in units of km/s. Note that the first two lines are
considered as comment lines when MG-MAMPOSSt reads the data.

Input/output file names and locations are indicated in gomamposst_x.inp; in
particular, input files correspond to the first two lines. The role of
pars_test.txt (second line of gomamposst_x.inp) will be described in detail
within Section 5, while the output files are summarized in Section 7.

##  3 Generalities

MG-MAMPOSSt operates relying on the input of the projected phase space
$(R,v_{\text{z}})$ of the cluster member galaxies. Here, $R$ is the projected
distance from the cluster center at which a galaxy is seen by the observer,
and $v_{\text{z}}$ the velocity measured along the l.o.s. in the rest frame of
the cluster. Only those tracers lying in the projected range
$[R_{\text{low}},R_{\text{up}}]$ will be considered in the fitting procedure
of the model parameters (see Section 5). As mentioned before, the current
version of the code assumes parametric expressions for the various kinematical
quantities. Moreover, the 3-dimensional velocity distribution is taken to be a
Gaussian. The latter assumption has been well-tested through cosmological
simulations, as explained in the original MAMPOSSt paper (MAM13).

The output likelihood is computed for a given set of models and parameters
from the solution of the spherical Jeans’ equation (see e.g. Ref. [7]),

$\sigma^{2}_{r}(r)=\frac{1}{\nu(r)}\int_{r}^{\infty}{\exp\left[2\int_{r}^{s}{\frac{\beta(t)}{t}\text{d}t}\right]\nu(s)\frac{\text{d}\Phi}{\text{d}s}\text{d}s},$
(1)

projected in the phase space. The above equation captures the main necessary
input required for the MAMPOSSt method. In particular, the gradient of the
gravitational potential, which in turn contains information about the model of
gravity, as well as the velocity anisotropy profile ($\beta(r)$, eq. [2]) and
the projected number density of tracers ($\nu(r)$). The parametrisations for
all these quantities implemented within the code are defined in the following
sections.

The current version of the code can handle up to a seven-dimensional parameter
space, with two parameters defining the mass profile, one parameter for the
velocity anisotropy profile, one for the number density profile, and finally,
two parameters related to the modified gravity framework. Each parameter can
be either treated as free in the fitting procedure, or it can be assigned pre-
defined values (see also Section 5). The list of possible parameters and the
relevant equations for the modified gravity parametrizations implemented are
summarized in Table 1. MG-MAMPOSSt additionally returns a group of data files,
stored in the Output folder, which contain the binned and fitted projected
distribution of galaxies ($N(R)$), and the binned/fitted velocity dispersion
profiles (VDPs), as explained in Section 7. In Figure 1 we present a scheme
summarizing all the relevant aspects in order to perform a complete run of the
code.

Figure 1:  Sketch of the MG-MAMPOSSt algorithm overviewing the hierarchy and
connection between the input/output files in the code. Detailed information
about each file is described in the text.

A general analytic expression of the likelihood function can be found in
MAM13.

##  4 Implemented modelling

###  4.1 Mass and velocity anisotropy profiles

In the original code of MAM13 there are several possible choices for the
modelling of the dark matter mass profile in the $\Lambda$CDM scenario. MG-
MAMPOSSt adds new parametrizations to handle popular modifications of gravity.
At the moment, all the implemented non-standard profiles rely on the Navarro-
Frenk-White (NFW, Ref. [8]) mass density profile to model the matter density
distribution, which has been shown to provide a good description for simulated
and observed galaxy clusters, both in GR and in modified gravity (see e.g.
Refs. [9, 10, 11]). Nevertheless, other mass models are going to be included
in the code with upcoming versions. Since the mass profile, or equivalently,
the gravitational potentials, enters only in the expression of the radial
velocity dispersion eq. (1), the implementation of new models can be performed
directly by the user, modifying the subroutines where the above equation is
involved in the auxiliary file MAM.f. In particular, the functions sr2int(alr)
and fa(tlog) are the only parts of the MG-MAMPOSSt code where parametrizations
for the mass profile appear.

In the current version, ten choices of mass profiles are implemented: seven
(GR) standard models, and three modified NFW models for MG frameworks. Each
option in the code is identified by a specific acronym as described also in
Section 5. The implemented profiles are the following: NFW, Hernquist (Ref.
[12]), Pseudo Isothermal Elliptical Mass Distribution (PIEMD, Ref. [13]),
Burkert (Ref. [14]), Softened Isothermal Sphere, Einasto with exponent $m=5$
(Ref. [15]), "Linear Horndeski" NFW, Vainsthein screening (VS) NFW, Chameleon
Screening (CS) NFW. We have further implemented a generalized NFW profile
(gNFW hereafter) in GR, where the slope of the density is controlled by an
exponent $\gamma$, equal to 1 for the standard NFW of eq. (11), which is now
treated as a free parameter to be fitted together with the mass profile
parameters. In this last case, the exponent $\gamma$ is given by the parameter
$\mathcal{A}_{1}$. Further details are given in Section 5 where all options in
the code are over viewed.

Except for the modified gravity profiles and the gNFW, each model is
characterized by two parameters, namely a virial radius $r_{200}$, and a scale
radius $r_{\text{s}}$. Though we do not give a detailed description of those
profiles here, in the next subsection we will nevertheless present the key
expressions for the modified gravity models implemented and the necessary
parameters that need to be defined in the code.

In galaxy-cluster kinematics analyses, one of the major unknowns of the mass
profile reconstruction is the so-called velocity anisotropy profile

$\beta(r)=1-\frac{\sigma^{2}_{\theta}(r)+\sigma^{2}_{\phi}(r)}{2\sigma^{2}(r)},$
(2)

which accounts for the difference between the components of the galaxies’
velocity dispersion. As for the mass profile, in the current version of MG-
MAMPOSSt the velocity anisotropy is modeled as a parametric profile to be
fitted within the procedure. Eight models have been implemented so far,
identified by a specific label as for the mass modeling. Below we list the
expressions of those models:

  * 1

Constant anisotropy: $\beta(r)=\beta_{C}$,

  * 2

Mamon & Lokas profile (Ref. [7]):

$\beta_{\text{ML}}(r)=\frac{1}{2}\,\frac{r}{r+r_{\beta}}.$ (3)

  * 3

Osipkov-Merritt model (Refs. [16, 17]):

$\beta_{\text{OM}}(r)=\frac{r^{2}}{r^{2}+r^{2}_{\beta}}$ (4)

  * 4

Simplified Wojtak:

$\beta_{\text{SW}}(r)=\frac{r^{3/2}}{r^{3/2}+r_{\beta}^{3/2}}.$ (5)

  * 5

Tiret model (Ref. [18]):

$\beta_{\text{T}}(r)=\beta_{\infty}\,\frac{r}{r+r_{\beta}}.$ (6)

  * 6

Modified Tiret ("Opposite") model (see e.g. Ref. [19]):

$\beta_{\text{O}}(r)=\beta_{\infty}\,\frac{r-r_{\beta}}{r+r_{\beta}}.$ (7)

  * 7

Generalized Osipkov-Merritt model (Refs. [16, 17], see also Ref. [20]):

$\beta_{\text{gOM}}(r)=\beta_{0}+(\beta_{\infty}-\beta_{0})\frac{r^{2}}{r^{2}+r^{2}_{\beta}}$
(8)

  * 8

Generalized Tiret model (Ref. [18], see also Ref. [20]):

$\beta_{\text{gT}}(r)=\beta_{0}+(\beta_{\infty}-\beta_{0})\frac{r}{r+r_{\beta}}$
(9)

Note that models from one to four are characterized by only one free
parameter, while models five and six have two parameters: $\beta_{\infty}$,
the velocity anisotropy for $r\to\infty$, and $r_{\beta}$, the characteristic
scale radius (anisotropy radius). Finally, the last two models requires an
additional third normalization parameter $\beta_{0}$, which is the value of
anisotropy at $r=0$ However, in MG-MAMPOSSt the anisotropy radius for models
five to eight is assumed to be proportional to the scale radius of the mass
profile $r_{\beta}\equiv\alpha\,r_{\text{-2}}$, where $\alpha=0.5,1.51,1$ for
Hernquist model, Burkert model, and all the other models, respectively. This
leaves a maximum of two free parameters for the orbit anisotropy to be fitted.
Note that in the case of constant anisotropy, Tiret (normal, modified or
generalized) and generalized Osipkov-Merritt, the normalization parameters
used within the code are $\mathcal{A}_{C}=(1-\beta_{C})^{-1/2}$,
$\mathcal{A}_{\infty}=(1-\beta_{\infty})^{-1/2}$ and
$\mathcal{A}_{0}=(1-\beta_{0})^{-1/2}$ instead of $\beta_{C}$,
$\beta_{\infty}$ and $\beta_{0}$ , respectively.

###  4.2 Number density profiles

The computation of the MG-MAMPOSSt likelihood requires the knowledge of the
number density distribution of the tracers. Generally, this can be
reconstructed directly from the phase space before running the code, but one
can choose to fit the profile together with the anisotropy and the
gravitational potential. Currently, MG-MAMPOSSt implements three models for
$\nu(r)$, namely NFW model, eq. (11), Hernquist model and beta-model (Ref.
[21]), which can be selected from the input file pars_test.txt as explained in
Section 5. All the models are specified by a normalization factor and a scale
radius $r_{\nu}$. Moreover, the beta-model requires a (fixed) negative
exponent $n_{b}$ which is set from input. Note that in eq. (1), the
normalization factor of $\nu(r)$ cancels out. Thus, the only relevant
parameter of the profile is the scale radius $r_{\nu}$.

###  4.3 NFW models for Vainsthein and Chameleon gravity

The code is equipped with the two most popular and observationally viable
classes of dark energy models beyond GR based on a single, extra scalar field.
These correspond to the so–called chameleon models (“Models I") and Beyond
Horndeski/DHOST models (“Models II"). Common ground between the two families
of models is the presence of the extra dynamical scalar degree of freedom
($\phi$) which introduces a new gravitational force. However, the structure of
the fifth force in each of them exhibits different characteristics. For a
detailed exposition of these models and the associated equations we refer to
our main paper (Ref. [6]), as well as to the original papers where the models
were first introduced (Refs.[22, 23, 24]). Note that there is a third
possibility implemented in the current version of the code: the "Linear
Horndenski" models (“Models 0"), characterized by a Yukawa-type fifth force
without the inclusion of any screening mechanism. The model, first introduced
in the code by Ref. [25], is actually a sub-case of the chameleon class,
specified by two free parameters, namely the characteristic mass of the field
$\mathcal{A}_{1}=m$ and the strength of the coupling
$\mathcal{A}_{2}=\mathcal{Q}$. The difference between Models I and this linear
framework is in the full treatment of the Yukawa suppression, which is
neglected in the chameleon screening approximation. Models 0 will be removed
in the upcoming versions of the code as its phenomenology at cluster scales is
already encompassed by Models I.

As shown in Ref. [6], internal kinematics alone is generally not enough to
provide stringent bounds on the modified gravity parameters, due to the strong
degeneracy between model parameters. For this reason, MG-MAMPOSSt gives the
possibility to include a simulated lensing information to your kinematics
analysis in modified gravity, a feature which is particularly useful for
forecasting the constraining power of the method. The structure of the lensing
Likelihood depends on the chosen model of gravity and the details are given
below.

Models I: This is the family of chameleon models; they can be selected from
pars_test.txt as M(r)=mNFW_GC. Their key characteristic is the existence of a
scale $S$, the so–called screening radius, below (beyond) which the force is
suppressed (unsupressed). The Newtonian potential gradients in the Jeans
equation assumes the following form:

$\frac{\text{d}\Phi}{\text{d}r}=\frac{GM(r)}{r^{2}}+\frac{\mathcal{Q}}{M_{\text{P}}}\frac{\text{d}\phi}{\text{d}r},$
(10)

with $M(r)$ the total mass enclosed within a sphere of radius $r$,
$M_{\text{P}}=(8\pi G)^{-1/2}$ is the reduced Planck mass, $G$ is the
gravitational constant and $\mathcal{Q}$ a dimensionless coupling defining the
strength of the chameleon force. For the scalar field $\phi$, we implement an
analytic, approximated solution of eq. (10), which is explained in e.g. Ref.
[26], assuming the NFW model for the matter distribution. The NFW profile is
defined through

$\rho(r)=\frac{\rho_{\text{s}}}{r/r_{\text{s}}(1+r/r_{\text{s}})^{2}},$ (11)

where the characteristic density $\rho_{\text{s}}$ is a function of
$r_{\text{s}}$ and $r_{200}$. With this ansatz the field profile acquires the
form

$\phi(x)=\begin{cases}\phi_{\text{s}}[x(1+x)^{2}]^{1/(n+1)}\simeq 0&r\leq S\\\
-\displaystyle{\frac{\mathcal{Q}\rho_{\text{s}}r_{\text{s}}^{2}}{M_{\text{P}}}}\,\displaystyle{\frac{\ln(1+x)}{x}}-\displaystyle{\frac{C}{x}}+\phi_{\infty}&r>S.\end{cases}\
$ (12)

In the last expression, $x\equiv r/r_{\text{s}}$, $\phi_{\infty}$ is the
background value of the field and $C$ is an integration constant. Finally,
$\phi_{\text{s}}<<\phi_{\infty}$ is the field value in the core of the matter
distribution. Under the approximation that the interior solution is everywhere
close to zero, $C$ is fully determined by the free parameters of this gravity
models (i.e. the fifth-force coupling $\mathcal{Q}$ and the cosmological value
$\phi_{\infty}$) and the mass profile parameters:
$C=C(\mathcal{Q},\rho_{s},r_{s},\phi_{\infty})$. The _screening radius_ $S$
represents the transition scale between the regime where the force is screened
(i.e suppressed), and the one where the force is unscreened (i.e fully
operating). The screening radius has a similar functional dependence,
$S=S(\mathcal{Q},\rho_{\rm{s}},r_{\rm{s}},\phi_{\infty})$. We refer to Ref.
[6] for the explicit expressions. Note that, since the Yukawa suppression is
neglected in this approximation, in order to avoid an overestimation of the
fifth force the additional contribution to the gravitational potential is set
to zero if $S>1.5\,R_{\text{up}}$, where $R_{\text{up}}$ is the maximum radius
up to which galaxies are considered in the projected phase space of the MG-
MAMPOSSt fit, as explained in the next sections.

If the number of iteration for the coupling constant $n_{A2}$ (see Section 5)
is set to -1 for this class of models, the chameleon $f(\mathcal{R})$ sub-case
is selected, which correspond to a fixed $\mathcal{Q}=\sqrt{1/6}$. In this
framework, the only free parameter is the so called scalaron field
$f_{\mathcal{R}}=\partial f/\partial\mathcal{R}$, where $\mathcal{R}$ is the
Ricci scalar. The scalaron is connected to the chameleon field through:

$\phi_{\infty}=-\sqrt{\frac{3}{2}}\ln(1+f_{\mathcal{R}0})M_{\text{P}}c^{2}_{l}\,,$
(13)

where $c_{l}$ is the speed of light. In MG-MAMPOSSt the first modified gravity
parameter, generally dubbed as $\mathcal{A}_{1}$, corresponds to the
background field over the Planck mass $\phi_{\infty}/M_{\text{P}}$, given in
units of $10^{-5}\,c_{l}^{2}$, while the second modified gravity parameter is
$\mathcal{A}_{2}=\mathcal{Q}$ (see 5).

As for the lensing simulation, in chameleon gravity the fifth force doesn’t
affect photon propagation, except for a negligible rescaling of the
gravitational constant $G$. Thus, lensing analyses are sensitive only to the
GR parameters of the mass profile $r_{\text{s}},\,r_{200}$. As such, the
additional information in MG-MAMPOSSt is provided as a bivariate Gaussian
$P(r_{\text{s}},\,r_{200})$ whose parameters (namely, the means and the
covariance matrix) are given as input by the user in the Options.txt file (see
Section 6).

Models II: This class of models corresponds to the so–called Vainshtein
breaking theories. The corresponding modified gravitational potential can be
selected in MG-MAMPOSSt with M(r)=mNFW_BH. Here, the fifth force within the
spherical overdensity is only described by a single coupling parameter
$\mathcal{A}_{1}=Y_{1}$. The Newtonian potential is defined through (e.g.
[27])

$\frac{\text{d}\Phi}{\text{d}r}=\frac{GM(r)}{r^{2}}+\frac{Y_{1}G}{4}\,\frac{\text{d}^{2}M(r)}{\text{d}r^{2}}.$
(14)

The dimensionless coupling $Y_{1}$ is the additional free parameter that can
be fitted within the code. We note that, the second derivative
$d^{2}M(r)/dr^{2}$ in equation (14) is expressed in terms of the density
profile and its first derivative through
$\text{d}M(r)/\text{d}r=4\pi\rho(r)r^{2}$.

In these theories, the weak lensing potential, defined through
$\Phi_{\text{lens}}=\frac{1}{2}(\Phi+\Psi)$, is explicitly modified through
the appearance of a new term in the relativistic potential $\Psi$ as

$\displaystyle\frac{\text{d}\Psi}{\text{d}r}=\frac{GM(r)}{r^{2}}-\frac{5\,Y_{2}G}{4\,r}\,\frac{\text{d}M}{\text{d}r},$
(15)

with $Y_{2}\equiv\mathcal{A}_{2}$ a dimensional coupling. Assuming the NFW
mass model, the gradients of the potentials $\Phi,\,\Phi_{\text{lens}}$ read:

$\frac{\text{d}\Phi}{\text{d}r}=\frac{GM_{\text{NFW}}(r)}{r^{2}}\left\\{1+\frac{Y_{1}}{4}\,\frac{r^{2}(r_{\text{s}}-r)/(r_{\text{s}}+r)^{3}}{[\ln(1+r/r_{\text{s}})-r/(r_{\text{s}}+r)]}\right\\},$
(16)
$\frac{\text{d}\Phi_{\text{lens}}}{\text{d}r}=\frac{GM_{\text{NFW}}(r)}{r^{2}}\left\\{1+\frac{\left[Y_{1}(r_{\text{s}}-r)-5\,Y_{2}(r_{\text{s}}+r)\right]}{8\,[\ln(1+r/r_{\text{s}})-r/(r_{\text{s}}+r)]}\frac{r^{2}}{(r_{\text{s}}+r)^{3}}\right\\},$
(17)

where $M_{\text{NFW}}(r)$ is the standard NFW mass profile in GR, obtained by
integrating eq. (11).

As for Models I, in MG-MAMPOSSt we have implemented the possibility to
forecast the constraints obtainable on a given MG model when adding a lensing
information to the kinematics analysis only. In the case of Vainsthein
screening, the code is able to perform a weak lensing simulation which
produces a likelihood
$\mathcal{L}_{\text{lens}}(r_{\text{s}},r_{200},Y_{1},Y_{2})$ to be added to
the kinematic likelihood. The simulation builds a (GR) mock reduced tangential
shear profile assuming again a NFW model, and computes the log-likelihood as

$\ln\mathcal{L}_{\text{lens}}(\bm{\theta}_{\text{l}})=-\frac{1}{2}\sum_{i=1}^{N_{\text{b}}}\frac{\left[\langle
g_{\text{t}}(R_{i})\rangle-\langle
g_{\text{t,vs}}(R_{i}|\bm{\theta}_{\text{l}})\rangle\right]^{2}}{\sigma^{2}_{\text{l},i}}.$
(18)

In the above equation, $\langle g_{\text{t}}(R_{i})\rangle$ is the simulated
averaged reduced tangential shear profile at projected position $R_{i}$,
$\langle g_{\text{t,vs}}(R_{i})|\theta_{\text{l}})\rangle$ is the theoretical
profile computed for the Vainsthein model and the set of parameters
$\bm{\theta}_{\text{l}}=(r_{\text{s}},r_{200},Y_{1},Y_{2})$. The number of
bins is set to be $N_{\text{b}}=10$, in agreement with current lensing surveys
(e.g. Ref. [28]) in the projected radial range
$(0.12\,R_{\text{up}},2.9\,R_{\text{up}})$, where $R_{\text{up}}$ is the
maximum radius up to which tracers in the phase space are considered in the
MG-MAMPOSSt fit (see Section 5). As for the uncertainties,

$\sigma^{2}_{\text{l},i}=\sigma^{2}_{\text{e},i}+\sigma^{2}_{\text{lss}},$
(19)

where
$\sigma^{2}_{\text{e},i}=\sigma^{2}_{\text{g}}/[\pi(\alpha^{2}_{\text{up}}-\alpha^{2}_{\text{low}})n_{\text{g}}]$
is the noise due to the intrinsic ellipticity $\sigma^{2}_{\text{g}}$ of the
sources lying within an annulus between the angles $\alpha_{\text{low}}$ and
$\alpha_{\text{up}}$, and $\sigma^{2}_{\text{lss}}$ expresses the effect of
the uncorrelated projected large scale structure. Finally, $n_{\text{g}}$ is
is the average number of source galaxies per $\text{arcmin}^{2}$. The values
of $\sigma_{\text{g}},\,\sigma_{\text{lss}},\,n_{\text{g}}$ can be customized
in Options.txt

Family of parameters | Parameters | Relevant equations
---|---|---
Total matter profile (e.g. NFW) |  $r_{200}$, $r_{\rm{s}}$ |  (11)
Number density profile (e.g. NFW) | $r_{\nu}$ |  (11)
Velocity anisotropy profile (e.g. Tiret) |  $\beta_{\infty}$ or $r_{\beta}$, $\beta_{0}$ |  (6)
Model 1 (Chameleon) |  $Q$, $\phi_{\infty}$ |  (10), (12)
Model 2 (Vainshtein) |  $Y_{1}$, $Y_{2}$ |  (14), (15)
Table 1: The main parameters associated with the modelling of the dark matter
and velocity anisotropy profile, as well as the models of gravity implemented
in MG-MAMPOSSt. They can be either pre-set in the code or left as free
parameters to be fit when confronting against real or simulated data.

##  5 Overview of input parameters

The file input_pars/ pars_test.txt contains the input parameters for the MG-
MAMPOSSt procedure. All the parameters should be provided in a format "<label>
= <value>". The "<label>"s are mandatory while the "<value>"s, if not given,
are set by default, except for the case of the guess values of the model’s
free parameters. Inline comments can be inserted with a "!". Input parameters
are divided into four main groups:

1\. Number of iterations,

2\. Guess values,

3\. Model/setup options,

4\. Parameter limits .

1\. Number of iterations (integer): For the grid search mode, they represent
the number of points in each free parameter over which the likelihood is
computed. If set to $0,1$, the parameter is fixed to its guess value, except
for specific cases described below. The total number of likelihood evaluations
is given by

$N_{\text{tot}}=\prod_{i}n_{i},$ (20)

where $i$ runs over the number of parameters which are fitted by the code.
Currently, the maximum number of free parameters are seven.

When MG-MAMPOSSt is in MCMC mode, if number of iterations is different from
zero, then the corresponding parameters are optimized within the chain.
Otherwise the parameters are fixed to the guess value.

In some specific cases other options are available; we will briefly describe
them when listing all the input variables below.

2\. Parameter guess values (real): They serve as an initial guess for the MG-
MAMPOSSt fit both in grid and MCMC mode. If the corresponding number of
iterations is set to zero, or one when in the grid-search option, the
parameter guess is kept fixed within the code.

3\. Model options (real, integer, string): This family of values allow to
select cosmological environment such as the value of the Hubble parameter and
the average redshift of the cluster, the mass or number density profiles, the
velocity anisotropy profiles, as well as the optimization algorithm choices.

4\. Parameter limits (real): They are labelled as $X^{\text{low}}$,
$X^{\text{up}}$, where

$X\in\\{r_{200},r_{s},r_{\nu},\beta,\mathcal{A}_{1},\mathcal{A}_{2}\\},$ (21)

and are assumed as the limiting values of the parameter space. In the grid-
search mode, they are considered only if the option kpro $=1$ (see Options.txt
below). The grid exploration starts from the guess values and the likelihood
is computed over a logarithmically-spaced set of points where the width is set
such that the extremes corresponds to $X^{\text{low}}$ and $X^{\text{up}}$. If
kpro $=0$, the width is fixed to a default value chosen optimally by the code
for each parameter.

In the MCMC mode, the range within the limiting values should be assumed as a
flat prior, i.e.

$P(X)=\frac{1}{X^{\text{up}}-X^{\text{low}}}\begin{dcases}0&X<X^{\text{low}}\\\
1&X^{\text{low}}\leq X\leq X^{\text{up}}\\\ 0&X>X^{\text{up}}.\\\
\end{dcases}$ (22)

The pre-factor $1/(X^{\text{up}}-X^{\text{low}})$ ensures the normalization of
the probability distribution. One should be warned that, for Vainsthein
screening, some values of $Y_{1},\,Y_{2}$ can give rise to a negative
effective cluster mass at some radii, which is physically inconsistent. Thus,
in the MCMC run, in addition to the flat prior all the couples
$(Y_{1},\,Y_{2})$ related to negative masses are discarded. If kpro $\neq 1$
the parameter limits are chosen by default.

Below we list all the input parameters in the same order as they appear in
pars_test.txt:

  * -

nr200 (integer): Number of iteration for the virial radius $r_{200}$.

  * -

nrc (integer): Number of iteration for the tracers scale radius $r_{\nu}$. For
$\texttt{nrc}=-1$, it assumes that $r_{\nu}$ is equal to the scale radius of
the mass profile $r_{s}$ in the fit (option: “Light Follows Mass"). If
$\texttt{nrc}=-2$ the code first fits the number density profile alone and the
best fit found is then used as fixed value in the MG-MAMPOSSt procedure.

  * -

nrs (integer): Number of iteration for the mass profile scale radius $r_{s}$.
For $\texttt{nrs}=-1$ assumes that $r_{s}$ is equal to the scale radius of the
number density profile $r_{\nu}$ in the fit (option: “Mass Follows Light"). If
$\texttt{nrs}=-2$ the mass scale radius is computed by using the theoretical
relation of Ref. [29] ($\Lambda$CDM option).

  * -

nbeta (integer): Number of iteration for the first anisotropy parameter. If
$\texttt{nbeta}=-1$ the anisotropy profile is forced to be a Mamon & Lokas
profile (eq. 3) with $\beta\equiv r_{\beta}=r_{s}$. If $\texttt{nbeta}=-2$ the
Hansen-Moore model of Ref. [30] is assumed ($\beta(r)$ related to
$\rho_{m}(r)$).

  * -

nA1 (integer): Number of iteration for the first MG parameter
$\mathcal{A}_{1}$.

  * -

nA2 (integer): Number of iteration for the second MG parameter. For the case
of general chameleon gravity, $\mathcal{A}_{2}$ corresponds to the coupling
constant $\mathcal{Q}$. If $n_{A2}=-1$ it forces the case of $f(\mathcal{R})$
gravity ($\mathcal{Q}=1/\sqrt{6}$).

  * -

nbeta2 (integer): Number of iteration for the second anisotropy parameter,
used only if beta(r)=gT,gOM

  * -

r200g (real): guess starting value of the characteristic "virial" radius for
the mass profile $r_{200}$, measured in units of $[\text{Mpc}]$.

  * -

rcg(real): guess starting vaule for the scale radius of the number density
profile $r_{\nu}$, units of $[\text{Mpc}]$.

  * -

rsg(real): guess starting value for the scale radius of the mass profile
$r_{s}$, in units of $[\text{Mpc}]$.

  * -

betag (real): starting guess for the first velocity anisotropy profile
parameter $\beta$ . If the selected model is "constant", "Tiret", "Modified
Tiret", "generalized Tiret" or "Generalized Osipkov-Merritt" (see below) the
parameter is dimensionless, and corresponds to $\mathcal{A}_{\infty/0}$.
Otherwise it will be interpreted in units of $[\text{Mpc}]$.

  * -

A1g (real): Initial guess for the first modified gravity parameter
$\mathcal{A}_{1}$ (see Section 4). If the mass model is mNFW_GC, the parameter
is given in units of $10^{-5}c_{l}^{2}$, where $c_{l}$ is the speed of light
(e.g. a value of $10$ means $10^{-4}c_{l}^{2}$). For the case of gNFW model,
$\mathcal{A}_{1}=\gamma$ is the free exponent which characterizes the profile.

  * -

A2g (real): Initial guess for the second modified gravity parameter
$\mathcal{A}_{2}$ (see Section 4).

  * -

beta2g (real): starting guess for the second velocity anisotropy profile
parameter $\beta_{0}$. Used only for the generalized Tiret and the generalized
Osipkov-Merritt models. It is always dimensionless.

  * -

H0 (real): The value of the Hubble parameter evaluated at redshift $z=0$,
measured in units of $[\text{km}\text{\, s}^{-1}\,\text{Mpc}^{-1}]$. Default
is 70.

  * -

za (real): Average redshift of the cluster’s center of mass. Default is 0.0.

  * -

Olam (real): Value of the $\Omega_{\Lambda}$ density parameter today. Default
is 0.7.

  * -

Omegam (real): Value of the $\Omega_{\text{m}}$ density parameter today.
Default is 0.3.

  * -

Rlow (real): Inner projected radius $R_{\text{low}}$, defining the minimum
projected radius from the cluster center (given in $[\text{Mpc}]$) at which
galaxies should be considered in the MG-MAMPOSSt fit. Default is 0.05.

  * -

Rup (real): Outer projected radius $R_{\text{up}}$, defining the maximum
projected radius from the cluster center (given in $[\text{Mpc}]$) at which
galaxies should be considered in the MG-MAMPOSSt fit. Default is the
(mandatory) value of r200g.

  * -

N(R) (string): Number density profile model. It selects the model for the
tracer’s projected number density profile. In the current version three
possible choices are allowed, indicated by the following labels:
pNFW for projected NFW,
pHer for projected Hernquist,
beta for beta-model.
Default is pNFW.

  * -

al (real): Exponent for beta-model. The beta-model choice requires an
additional input exponent which should be negative. The input parameter is
ignored otherwise. Default is 0.0.

  * -

M(r) (string): Mass profile/gravitational potential model. It selects the
allowed parametrisations of the mass profile to be used in the Jeans’
equation. The labels corresponding to the implemented models are:
NFW for NFW,
Her for Hernquist,
PIEMD for PIEMD,
Bur for Burkert,
SoftIS for Softened Isothermal Sphere,
Eis5 for Einansto with $m=5$,
gNFW for generalized NFW,
mNFW_LH for modified NFW in Linear Horndeski gravity,
mNFW_BH for modified NFW in Vainsthein screening (Beyond Horndeski),
mNFW_GC for modified NFW in general Chameleon screening gravity.
Details on each modified gravity model are given in Section 4. Default is NFW.

  * -

Beta(r) (integer): Anisotropy profile model. It selects the velocity
anisotropy $\beta(r)$ in the Jeans’ equation. The allowed profiles (see
Section 4) are: C for constant anisotropy,
ML for Mamon & Lokas,
OM for Osipkov-Merritt,
WJ for simplified Wojtak,
T for Tiret,
O for modified Tiret (Opposite),
gT for generalized Tiret,
gOM for generalized Osipkov-Merritt.
Default is C.

  * -

rcut (real): Truncation radius. Value of the truncation radius needed in the
pseudo-isothermal elliptical mass distribution (PIEMD). It is ignored for
other mass profile models.

  * -

FASTMODE (integer): If equal to 1, the likelihood is estimated by using a grid
of values (default 60 points) in the phase space $(R,v_{z})$ of the galaxies
and then bispline-interpolating over the data points.

  * -

OPT (integer): Optimization choice. If required, an optimization algorithm is
launched before the grid/MCMC likelihood computation to find the minimum of
$-\ln\mathcal{L}$. Eventually, the resulting best fit parameters are used as
guess values for the subsequent parameter space exploration. Currently, three
choices are available: BOBYQA ($\texttt{OPT}=0$), NEWUOA ($\texttt{OPT}=1$) or
POWELL ($\texttt{OPT}=2$)444We point out that the POWELL algorithm does not
work efficiently when considering the gT and gOM models for the velocity
anisotropy, due to round-off errors in the integration routines, while the
other two alghoritms do.. For details, see MAM13. If $\texttt{OPT}=-1$ the
optimization is skipped.

  * -

screen (integer): Screening mode (available only for M(r)=mNFW_LH). In linear
Horndeski, one can choose to adopt the $f(\mathcal{R})$ sub-case,
$\texttt{screen}=0$, where $\mathcal{A}_{2}$ is fixed to $\sqrt{1/6}$. In this
framework, there is the possibility to include a model-dependent screened
$f(\mathcal{R})$ model with Hu&Sawicki functional form,implemented by assuming
a simple analytical approximation. The transition between the screened and
linear regime can be instantaneous ($\texttt{screen}=1$), or smoothed with an
additional parameter controlling the sharpness ( $\texttt{screen}=2$). For
$\texttt{screen}=-1$, the general linear Horndeski with two free parameters is
selected.

  * -

r2low: Lower limit for $r_{200}$ $[\text{Mpc}]$.

  * -

r2up: Upper limit for $r_{200}$ $[\text{Mpc}]$.

  * -

rclow: Lower limit for $r_{\nu}$ $[\text{Mpc}]$.

  * -

rcup: Upper limit for $r_{\nu}$ $[\text{Mpc}]$.

  * -

rslow Lower limit for $r_{\text{s}}$ $[\text{Mpc}]$.

  * -

rsup: Upper limit for $r_{\text{s}}$ $[\text{Mpc}]$.

  * -

blow: Lower limit for the first anisotropy parameter $\beta$.

  * -

bup: Upper limit for the first anisotropy parameter $\beta$.

  * -

A1low: Lower limit for the first modified gravity parameter.

  * -

A1up: Upper limit for the first modified gravity parameter.

  * -

A2low: Lower limit for the second modified gravity parameter.

  * -

A2low: Upper limit for the second modified gravity parameter.

  * -

b2low: Lower limit for the second anisotropy parameter $\beta_{0}$.

  * -

b2up: Upper limit for the second anisotropy parameter $\beta_{0}$.

Note that the parameter limits in the last block have as default values
$10^{-}4$ and $10^{2}$ for all the lower bounds and upper bounds respectively.

##  6 Operating Options

The file Options.txt contains various options and switches for the new
features in MG-MAMPOSSt. These are mostly related to the numerical analysis
and evaluation of the posterior likelihood. Notice that, the input parameters
can be binary integers (with values $0$ or $1$), integers*4 or reals*8. As for
the case of pars_test.txt, all the parameters must be given in a format
"<label> = <value>". The "<label>"s are mandatory while the "<value>"s, if not
given, are set by default.

  * •

nmcmc (binary). Select between grid-search mode ($=0$), and MCMC sampling
($=1$). The default value is $=1$.

  * •

Nsample (integer). Number of points in the MCMC run. Default is $400000$.

  * •

nlens (binary). Lensing information: If equal to 1, it adds to the kinematics
likelihood, a probability distribution simulating additional information such
as provided by a lensing mass profile reconstruction. For each set of values
of the parameters, the joint (log) likelihood is then computed as

$\ln\mathcal{L}_{\text{joint}}=\ln\mathcal{L}_{\text{dyn}}+\ln\mathcal{L}_{\text{lens}},$
(23)

where $\ln\mathcal{L}_{\text{dyn}}$ and $\ln\mathcal{L}_{\text{lens}}$ are the
log-likelihoods of the MG-MAMPOSSt procedure and the simulated lensing
distribution respectively. For linear Horndeski and Chameleon screening, where
photon propagation is not affected by the new degrees of freedom, the lensing
likelihood has the form of a bivariate Gaussian distribution
$\mathcal{L}_{\text{lens}}=P_{\text{lens}}(r_{\text{s}},r_{200})$, specified
by central values, standard deviations and correlation (see below). In
Vainsthein screening, where lensing mass profile is explicitly modified by the
MG parameters, the likelihood is computed by simulating a full tangential
shear profile, as explained in Ref. [6]. The default value is $=0$.

  * •

r200t (real). "True" value of the cluster’s virial radius (in unit of Mpc)
around which the lensing distribution is centered. This are mandatory. In case
it is not provided, the code switches to $\texttt{nlens}=0$.

  * •

rst (real). "True" value of the cluster’s scale radius (in unit of Mpc) around
which the lensing distribution is centered. This are mandatory. In case it is
not provided, the code switches to $\texttt{nlens}=0$.

  * •

delta1 (real). For Vainshtein screening, it represents the intrinsic
ellipticity of galaxies in the (weak) lensing simulations. For Chameleon
screening, it is the relative uncertainty on the virial radius in the Gaussian
distribution $\sigma_{r_{200}}/r_{200}$. Default for Vainsthein screening is
$=0.3$, default for Chameleon screening is $=0.1$.

  * •

delta2 (real). For Vainshtein screening, it represents the large scale
structure contribution to the errors in the (weak) lensing simulations. For
Chameleon screening, it is the relative uncertainty on the scale radius in the
Gaussian distribution $\sigma_{r_{\text{s}}}/r_{\text{s}}$. Default for
Vainsthein screening is $=0.005$, default for Chameleon screening is $=0.3$.

  * •

delta3 (integer/real). For Vainshtein screening, it represents the number of
galaxies per arcminute square in the (weak) lensing simulations. For Chameleon
screening, it is the correlation $\rho$ in the Gaussian distribution. Default
for Vainsthein screening is $=30$, default for Chameleon screening is $=0.5$.

  * •

kpro (binary). If it is equal to 1, the parameter space exploration is made
over a given interval $[X^{\text{low}},X^{\text{up}}]$, where $X$ indicates a
generic free parameter, as mentioned in Section 5. Default is 1.

  * •

Nclust (integer). MG-MAMPOSSt allows for efficient statistical forecast
analyses of the constraints on the implemented MG models. In particular, it is
possible to input Nclust realizations of phase spaces at the same time to
compute directly the joint likelihood for a given set of parameters, obtained
from the combination of the likelihood from each single data-set. These data
files should be located in "/data" folder and named as datphys_<i>.dat, where
<i> labels the number of the file in ascending order, starting from 2 (e.g.
datphys_2.dat, datphys_3.dat, …). The file format is the same as the main
input datphys.dat (see Section 3). Default is $=1$. Note that in order to
obtain meaningful results using this option, all the data files should be a
realization of the same cluster (i.e. characterized by the same true values of
all parameters).

  * •

nsingle (binary). If it is equal to 1, the Nclust-clusters likelihood is
computed by simply multiplying by Nclust the likelihood from a single data-
set. Useful to forecast a fast estimation of the limiting behaviour of the
constraints when increasing the number of clusters. Default is $=0$.

Additional Options

  * •

istop (binary). When equal to 1, the program stops after the preliminary
optimization algorithm if the relative difference between the best fit found
and the guess values is larger than
$\vec{\epsilon}=(\epsilon_{r_{200}},\epsilon_{r_{\nu}},\epsilon_{r_{\text{s}}},\epsilon_{\beta},\epsilon_{\mathcal{A}_{1}},\epsilon_{\mathcal{A}_{2}},\epsilon_{\beta_{0}})$,
or if the relative difference of the logarithm of the likelihood, computed at
the best fit and at the guess values is larger than
$\Delta_{\ln(\mathcal{L})}$ (see below). Default value is 0.

  * •

teps (array, real). Threshold values for $\vec{\epsilon}$. Default values are:

$\epsilon_{i}=0.1.$

  * •

delik (real). Threshold value for the relative difference of the logarithmic
likelihood

$\Delta_{\ln(\mathcal{L})}=\frac{\ln(\mathcal{L})_{\text{bf}}-\ln(\mathcal{L})_{\text{guess}}}{\ln(\mathcal{L})_{\text{guess}}},$

where $\ln(\mathcal{L})_{\text{bf}}$ is the best fit likelihood according to
MG-MAMPOSSt and $\ln(\mathcal{L})_{\text{guess}}$ is the likelihood computed
at the parameters’ guess values. Default is $=0.2$.

  * •

nskip (binary). For the grid case, if equal to 1, the exploration starts from
the free parameters guess values even if the preliminary optimization is not
skipped. This could be useful in modified gravity frameworks where the MG-
MAMPOSSt likelihood presents more than one maximum, due to statistical
degeneracy between model parameters. One can be interested in knowing the
position of the global peak but sampling the likelihood only around a specific
local maximum. Default is $=0$

  * •

nres (binary). For the grid exploration without fixed bounds (kpro$=0$), it
selects a larger (0) or a smaller-size (1) grid step. Note that the steps are
different for each parameter and adjusted depending on the number of available
tracers in the fit, unless specified (see below). Default is $=0$

  * •

nequ (binary). For the grid exploration without fixed bounds (kpro$=0$), if
equal to 1, it removes the re-scaling of the grid steps by the number of
available tracers. Default is $=0$

  * •

nsame (binary). If equal to 1, likelihoods are computed on specific values of
the parameters, given from an external input file. The name of the file should
be MaxLik_input.dat, structured as a table with seven columns (one for each
parameter), following the same order as in the output likelihood MaxLik.dat,
i.e.
$r_{200},\,r_{\nu},\,r_{\text{s}},\,\beta,\,\mathcal{A}_{1},\,\mathcal{A}_{2},\,\beta_{0}$
(see Section 7).
This option only works if the MCMC mode (nmcmc$=1$) is selected. Default is
$=0$

##  7 Output

The MG-MAMPOSSt code produces several output files which are automatically
saved in the Output folder. Below we list the files in the same order as they
appear in the gomamposst_x.inp script, summarizing their content and focusing
in particular on the output likelihood/posterior555Note that we refer to
“likelihood” when no fixed bounds are assumed (kpro$=0$), i.e. no prior
information, and to “posterior” otherwise..

  * rnvn.dat. Projected distances $R_{i}$ (in Mpc), velocities $v_{z,i}$ (in km/s) of the tracers used in the fit and the associated probability $p(R_{i},v_{z,i})$ to find a galaxy in the phase space $(R_{i},v_{z,i})$ given the theoretical model. Note that the (log) likelihood is given by the sum of $\ln p$ over the tracers.

  * NRpar.dat. Single-row data file where the choice of the number density profile, the best fit value of the scale radius $r_{\nu}$ and the value of the negative exponent for the beta-model ($k_{\text{n}}=3$) are stored.

  * NRbin.dat. Binned projected number density profile $N_{\text{bin}}(R)$ where the projected positions are in units of Mpc. The third column lists the associated errors.

  * NRfit.dat. Fitted projected number density profile $N_{\text{fit}}(R)$ where the projected positions are in units of Mpc.

  * MaxLik.dat. File containing the tabulated likelihood/posterior computed by the MG-MAMPOSSt procedure. This is the main output of the code. The table is structured as follows: columns from 1 to 7 store the values of the parameters, ordered as $r_{200},\,r_{\nu},$ $r_{\text{s}},\,\beta,\,\mathcal{A}_{1},\,\mathcal{A}_{2},\,\beta_{0}$. Column 8 is $-\ln{\mathcal{L}}$, while the integer number in the last column is related to the choice of the anisotropy profile, according to the following: 0=’C’, 1=’ML’, 2=’OM’, 21=’gOM’, 3=’WJ’, 4=’T’, 41=’gT’, 5=’O’. The two rows at the end of the file highlight the minimum found by the exploration and the minimum found by the preliminary optimization respectively, which may not coincide for distributions with more than one local maximum.

  * svbin.dat. Binned velocity dispersion profile (VDP). Each row shows the following quantities: number of galaxies in the i-th bin, $R_{i}$ (Mpc), VDP$(R_{i})$ (km/s), VDP’s upper error (km/s), VDP’s lower error (km/s)

  * svfit.dat. Fitted VDP obtained by using the best fit values found by the MG-MAMPOSSt procedure. First column: projected position (in Mpc). Second column: VDP (km/s).

In addition, running ./script_runmam.sh produces plots of the marginalized
distributions for the free parameters if the MCMC exploration is selected.
Plots are generated with the plot.py Python3 script which makes use of the
GetDist package 666https://getdist.readthedocs.io/en/latest/index.html, a
license-free Python package developed by A. Lewis (Ref. [31]). For some mass
models, the script produces also a plot of the best fit mass profile with the
associated one sigma and two sigma bands. The models for which this feature is
implemented are NFW, Burkert, Hernquist, generalized NFW, and the modified
gravity NFW for Vainsthein and Chameleon screening.

##  8 A simple walk-through tutorial

In this section we provide a tutorial on how to use the main functionality of
the code to analyse kinematics and simulated weak-lensing data.

In particular, we show an example of a complete run for the Vainsthein
screening case (mNFW_BH) using a single data-set shipped together with the
code base. The data file is a mock phase space drawn from an isolated,
spherically symmetric distribution of collisionless particles whose velocity
dispersion is given by the solution of the Jeans’ equation for a certain mass
profile and velocity anisotropy model (see Refs. [32, 6]). In particular, the
data-set consists in $\sim 2700$ particles obtained from a NFW profile in GR
(i.e. the fiducial value of the MG parameters is 0) with
$r_{200}=1.4\,\text{Mpc}$, $r_{\text{s}}=0.33\,\text{Mpc}$ and $r_{\nu}\equiv
r_{\text{s}}$. The phase space is obtained requiring that $\sim 600$ tracers
are included in a sphere of radius $r_{200}$. As for the velocity anisotropy
profile, we adopt the Tiret model of eq. (6), with $\beta_{\infty}=0.5$,
$\beta_{0}=0.0$ and $r_{\beta}=r_{\text{s}}$.

Let us first consider the kinematical data only, without the additional
lensing simulation. For comparison, we perform a grid and an MCMC parameter
space exploration. In Table 2 we indicate the computational time of the two
modes on a common laptop as a function of the number of free parameters.

As a first step, we open the gomamposst_x.inp script and adjust the name and
the location of the data-set file. In our case we type:
data/datphys_test.dat. Then, in input_pars/pars_test.txt, write all the
correct settings for the case of interest. For simplicity, we assume as guess
values the true values of the cluster’s parameters. The second block of
parameters in the input file should appear as:

⬇

************ free parameters initial guess values *********************

r200g = 1.41 ! mass profile r200 initial guess (Mpc)

rcg = 0.33 ! N(R) scale radius initial guess (Mpc)

rsg = 0.47 ! mass profile scale radius initial guess (Mpc)

betag = 1.41 ! Anisotropy initial guess,

! beta_C, a_ML, a_OM, a_W, beta_inf

A1g = 0.00 ! first MG parameter initial guess

A2g = 0.00 ! second MG parameter initial guess

beta2g = 1.0 ! Second Anisotropy parameter initial guess beta0 for

! gOM and gT models

************************************************************************

As for the cosmology, the synthetic halo is generated at average redshift 0
with $H=70\,\text{km/s\,Mpc}^{-1}$, $\Omega_{\Lambda}=0.7$ and
$\Omega_{m}=0.3$. For the upper and lower limits of the fit, we consider
$R_{\text{up}}=r_{200}$ and $R_{\text{low}}=0.05\,\text{Mpc}$. For a detailed
discussion about the physical meaning of this limits, see e.g. Refs. [1, 19,
25].

We consider the modified NFW model for Vainshtein screening gravity of eq.
(16), mNFW_BH, with a Tiret anisotropy profile ’T’ and a NFW for the tracers
number density distribution pNFW. We adopt the fast-mode and the optimization
algorithm NEWUOA, a derivative-free method developed by Ref. [33], OPT=1, to
find a preliminary best-fit. The remaining parameters are dummy in this case.
Thus, the third block of parameters is given by

⬇

************ Model options *******************************************

H0 = 70.0 ! Hubble constant at z=0

za = 0.0 ! average redshift of the cluster (needed to evaluate Hz)

Olam = 0.7 ! Omega lambda density parameter

Omegam = 0.3 ! Omega matter density parameter

Rlow = 0.05 ! Inner projected radius for sample selection (Mpc)

Rup = 1.41 ! Outer projected radius for sample selection (Mpc)

N(R) = pNFW

! model for number density profile N(R)

! Here we list the labels for N(R) and

! then the corresponding profile:

! pNFW=projected NFW / pHer=projected Hernquist /

! beta=beta-model

al = -0.00

! N(R) negative exponent (only if knfit=3)

M(r) = mNFW_BH

! mass profile model.

! Here we list the labels for M(r) and

! then the corresponding profile:

! NFW=NFW/ Her=Hernquist/ PIEMD=PIEMD/ Bur=Burkert/

! SoftIS=SoftIS/! Eis5=Einasto_m=5/

! mNFW_LH = mod_NFW linear Horndeski/

! mNFW_BH = mod_NFW beyond Horndeski/

! mNFW_GC= mod_NFW general chameleon/

! gNFW = generalized NFW

Beta(r) = T

! Anisotropy model.

! Here we list the labels for beta(r) and

! then the corresponding profile:

! C=Constant/ ML=Mamon&Lokas/

! OM=Osipkov-Merritt/ WJ=simplified Wojtak/

! T=simplified Tiret/ gT= generalized Tiret

! O= (Opposite) modified Tiret/

! gOM= generalized Osipkov-Merritt

rcut= 1.0

! PIEMD model rcut in Mpc (only used if M(r)=PIEMD)

FASTMODE= 1

! Run MG-MAMPOSSt in fast mode interpolating over data points

OPT= 1

! optimization algorithm: 0/1/2=bobyqa/newuao/powell

! -1 skip optimization

screen= 0

! for linear f(R) kmp.eq.7, one can decide to set

! an instantaneous transition between

! screeening and linear regime, by using the

! analytical approximation of Lombriser+12.

!-1/0/1/2=noscreen (general Hordenski)/noscreen f(R)/

! screen(instantaneous transition)

! /screen (arctan transition)

************************************************************************

Finally, one may select the parameter limits (priors) for the exploration in
the case of kpro$=1$. For this example run we choose the intervals as (see
Section 5 for the ordering details):

⬇

********* parameter limits *********************************************

r2low= 0.4 !r200 lower bound

r2up= 5.0 !r200 upper bound

rclow= 0.05 !rc lower bound

rcup= 3.9 !rc upper bound

rslow= 0.04 !rs lower bound

rsup= 3.9 !rs upper bound

blow= 0.5 !beta lower bound

bup= 7.1 !beta upper bound

A1low = -0.67 !first MG parameter lower bound

A1up = 7.2 !first MG parameter upper bound

A2low= -12.0 !second MG parameter lower bound

A2up= 12.0 !second MG parameter upper bound

b2low = 0.5 !second beta parameter lower bound

b2up = 7.1 !second beta parameter upper bound

************************************************************************

Before going to the Options.txt file, choose the free parameters in the first
group of quantities inpars_test.txt. We recall that, for the MCMC mode, a
value different from zero means that the corresponding parameter is optimized
within the chain, while in the grid mode this value indicate the number of
grid points where the likelihood is evaluated. Below the example setup for
five free parameters
$(r_{200},\,r_{\text{s}},\,\beta,\,\mathcal{A}_{1},\,\mathcal{A}_{2})$:

⬇

************** number of steps in free parameters ********************

nr200 = 40 ! number of steps for r200 fit

nrc = 0 ! number of steps for rc fit, scale radius of N(R)

! [if = 0 takes guess value]

! [if = -1 forces LfM, c_nu=c]

! [if = -2 fit N(R) outside MAMPOSSt]

nrs = 40 ! number of steps for rs fit, scale radius of M(r)

! [if = 0 takes guess value]

! [if = -1 forces MfL, c=c_nu]

! [if = -2 forces LCDM, c=c(M)]

nbeta = 40 ! number of steps for anisotropy parameter

! [if = -1 forces a_ML=r_s]

! [if = -2 forces Hansen+Moore]

nA1 = 40 ! number of steps in the first MG parameter

nA2 = 40 ! number of steps in the second MG parameter

! If equal to -1 force the case of chameleon f(R) gravity

nbeta2 = 0 ! number of steps for the second anisotropy parameter

************************************************************************

In Options.txt we can now select the grid or MCMC search switch; let’s first
consider the grid exploration, nmcmc $=0$. We furthermore set kpro$=1$ and
nlens=0 (i.e. we do not include the lensing simulation for now). We do not
modify the other parameters which are irrelevant in this first case.
Important: we stress again that while the <value> can be omitted to ensure the
default choice, all the <labels> described in Section 6 are mandatory, do not
remove them from the file.

The code can be now executed form a linux terminal by means of the command
lines mentioned in Section 2. Since the grid mode has been selected, no plots
are shown at the end of the run. After $\sim 0.1$ seconds ($\sim 1$ second for
normal mode) from execution, the terminal should print the following lines:

⬇

optimiz. results: r200,rc,rs,beta,A1,A2,cbe0,-logL

1.406 0.330 0.369 1.569 0.083 0.000 1.000 5908.378

likelihood of the initial value and Delta chi square:

r200,rc,rs,beta,A1,A2, -logL, Delta chi square:

1.410 0.330 0.470 1.410 0.00 0.000 1.000 5909.734 2.711

Then, the parameter exploration starts around the minimum found and the values
of the posterior/likelihood are written on MaxLik.dat for each set of
parameters777In Vainsthein screening, kinematics is not sensitive to $Y_{2}$,
which means that the likelihood does not depend on this parameter.. Here is an
example of the first lines of the tabulated likelihood for the case of five
free parameters:

0.42844 0.33000 0.04000 0.50000 -0.67000 -12.00100 1.00000 9331.69590 4
0.42844 0.33000 0.04580 0.50000 -0.67000 -12.00100 1.00000 9298.36094 4
0.42844 0.33000 0.05245 0.50000 -0.67000 -12.00100 1.00000 9265.74595 4
0.42844 0.33000 0.06005 0.50000 -0.67000 -12.00100 1.00000 9234.35828 4
$\dots$

The exploration is made by means of nested cycles over the parameters of
interest. The cycles are ordered in the following way, from the innermost to
the outermost: $r_{\text{s}}\to
r_{\nu}\to\beta_{0}\to\beta\to\mathcal{A}_{1}\to\mathcal{A}_{2}\to r_{200}$.
In Table 2 we can see how the computational time becomes extremely large when
increasing the number of free parameters, even when the fast mode is selected.
Moreover, the structure of the table makes the output file enormous and very
difficult to handle. Thus, the MCMC mode is suggested when the number of free
parameters is $\geq 4$. Note that in general the execution depends also on the
number of data points where the likelihood is computed. However, since for
real clusters the number of tracers can be around $\sim 10^{3}$ in very few
optimistic cases, this factor has a mild influence on the total time.

Computational time needed
---
| Number of model parameters to fit
Mode | 1 | 2 | 3 | 4 | 5
grid | 0.3 s | 10 s | 605 s |  $5\times 10^{4}$ s |  $\sim 2$ weeks
grid fast | 0.3 s | 2 s | 106 s |  $6\times 10^{3}$ s |  $\sim 2.5$ days
MCMC |  $1.6\times 10^{4}$ s |  $1.6\times 10^{4}$ s |  $1.7\times 10^{4}$ s |  $1.7\times 10^{4}$ s |  $1.7\times 10^{4}$ s
MCMC fast |  $1.3\times 10^{3}$ s |  $1.4\times 10^{3}$ s |  $1.4\times 10^{3}$ s |  $1.4\times 10^{3}$ s |  $1.4\times 10^{3}$ s
Table 2:  The required computational time in order to perform a complete run
of MG-MAMPOSSt for the various modes as a function of the free parameters to
be fitted. The grid search assumes 40 points for each parameter (i.e. $40^{n}$
total points with $n$ the number of free parameters), while the MCMC has been
made sampling $10^{5}$ points every time. The run has been performed on an
i7-8565U processor laptop (8th Gen) at $1.80$ GHz, with 16 GB RAM.

To select the MCMC mode, set nmcmc $=1$ in Options.txt. As for the number of
samples in the run, we choose $10^{5}$ points by typing Nsample $=100000$ in
the same file. Leave all the other options unchanged, save the file and run
again the executable ./script_runmam.sh. Now the tabulated likelihood looks
like as:

1.43094 0.33000 0.32708 1.52905 0.07640 -0.04110 1.00000 5894.40454 4
1.38378 0.33000 0.31959 1.67060 0.24586 -0.03153 1.00000 5894.64370 4
1.34847 0.33000 0.32015 1.60624 0.14262 -0.06625 1.00000 5895.73375 4
$\dots$

To include the weak lensing simulation in Vainsthein screening, open again
Options.txt and choose nlens $=1$, together with the related parameters. When
MG-MAMPOSSt is run, for each vector of parameters the probability distribution
is then computed with eq. (23). In this example, we consider r200t $=1.4$ Mpc,
rst $=0.33$ Mpc, togehter with delta1 $=0.25$, delta2 $=0.005$ and delta3
$=30$ for the intrinsic ellipticity, the large-scale structure noise and the
number of galaxies per $\text{arcmin}^{2}$, respectively. An important note is
in order here: the weak lensing mock analysis is performed by assuming a value
$z=0.44$ for the average cluster redshift. We further note that the background
cosmology is almost irrelevant for our simulation, as in the kinematics
analysis the phase spaces are always given in the cluster’s rest frame.

Figure 2: Marginalized distributions obtained from the output posterior of the
MG-MAMPOSSt procedure in Vainsthein screening. The distributions refer to a
run with five free parameters and kinematics+lensing analysis. The darker and
lighter shaded areas represent the $1\sigma$ and $2\sigma$ regions,
respectively

In the case of an MCMC exploration, at the end of the MG-MAMPOSSt routine, the
marginalized distributions are plotted, as shown in Figure 2 for the case of a
joint lensing+kinematic analysis in Vainsthein screening and with five free
parameters.

##  9 Summary

We provided a detailed reference guide for the key features and usage of the
MG-MAMPOSSt code. The code is based on the original MAMPOSSt method and code,
and presents with important new features aiming at testing large classes of
gravity models with the internal kinematics of galaxy clusters. What is more,
MG-MAMPOSSt is now capable to complement kinematical analyses with additional,
simulated lensing information towards consistency tests of gravity based on
the reconstructed kinematical and lensing mass. In addition, the code is well-
suited for analyses based on either synthetic or real data, and provides fast
optimisation techniques for parameter exploration.

Although MG-MAMPOSSt is equipped with the currently most general theories for
dark energy based on new scalar fields, there is certainly a big space for
improvement. For example, it would be interesting to extend the theories
currently implemented to new ones or add efficient parametrisations for the
gravitational field, or other number/matter density profiles. Also, a very
interesting, yet challenging task would be its extension for the case of other
physical systems such as stellar clusters and dwarf galaxies. In this regard,
we hope that the code will provide a useful tool for future analyses of
gravity within the communities of cosmology and astrophysics, and see further
extensions and applications.

## Acknowledgements

LP is partially supported by a 2019 "Research and Education" grant from
Fondazione CRT. The OAVdA is managed by the Fondazione Cleḿent Fillietroz-
ONLUS, which is supported by the Regional Government of the Aosta Valley, the
Town Municipality of Nus and the "United́es Communes valdotaines Mont-Eḿilius.
I. D. Saltas is supported by the Grant Agency of the Czech Republic (GAČR),
under the grant number 21-16583M. The authors further acknowledge S. Sartor
for useful comments and suggestions, as well as all the developers of the free
Fortran routines used in MG-MAMPOSSt. Credits are given in the header of each
routine.

## References

  * [1] G. A. Mamon, A. Biviano, and G. Boué, “MAMPOSSt: Modelling Anisotropy and Mass Profiles of Observed Spherical Systems - I. Gaussian 3D velocities,” MNRAS, vol. 429, pp. 3079–3098, Mar. 2013.
  * [2] T. Verdugo, M. Limousin, V. Motta, G. A. Mamon, G. Foëx, F. Gastaldello, E. Jullo, A. Biviano, K. Rojas, R. P. Muñoz, R. Cabanac, J. Magaña, J. G. Fernández-Trincado, L. Adame, and M. A. De Leo, “Combining strong lensing and dynamics in galaxy clusters: integrating mamposst within lenstool. i. application on sl2s j02140-0535,” A&A, vol. 595, p. A30, Oct. 2016.
  * [3] G. A. Mamon, J. Chevalier, A. J. Romanowsky, and R. Wojtak, “Dark matter inner slope and concentration in galaxies: from the fornax dwarf to m87,” in Galaxy Masses as Constraints of Formation Models (M. Cappellari and S. Courteau, eds.), vol. 311, pp. 16–19, Apr. 2015.
  * [4] G. A. Mamon and E. Vitral, “Mamposst-pm: a bayesian code for the joint analysis of line-of-sight and plane-of-sky kinematics of spherical systems,” MNRAS, 2022.  in preparation.
  * [5] E. Vitral and G. A. Mamon, “Does ngc 6397 contain an intermediate-mass black hole or a more diffuse inner subcluster?,” A&A, vol. 646, p. A63, Feb. 2021.
  * [6] L. Pizzuti, I. D. Saltas, and L. Amendola, “mg-mamposst: a code to test modifications of gravity with internal kinematics and lensing analyses of galaxy clusters,” MNRAS, vol. 506, pp. 595–612, 06 2021.
  * [7] G. A. Mamon and E. L. Łokas, “Dark matter in elliptical galaxies - II. Estimating the mass within the virial radius,” MNRAS, vol. 363, pp. 705–722, Nov. 2005.
  * [8] J. F. Navarro, C. S. Frenk, and S. D. M. White, “A Universal Density Profile from Hierarchical Clustering,” ApJ, vol. 490, pp. 493–+, Dec. 1997.
  * [9] K. Umetsu, “Cluster–galaxy weak lensing,” A&AR, vol. 28, Nov 2020\.
  * [10] S. Peirani, Y. Dubois, M. Volonteri, J. Devriendt, K. Bundy, J. Silk, C. Pichon, S. Kaviraj, R. Gavazzi, and M. Habouzit, “Density profile of dark matter haloes and galaxies in the HORIZON-AGN simulation: the impact of AGN feedback,” MNRAS, vol. 472, pp. 2153–2169, Dec. 2017.
  * [11] H. Wilcox, R. C. Nichol, G.-B. Zhao, D. Bacon, K. Koyama, and A. K. Romer, “Simulation tests of galaxy cluster constraints on chameleon gravity,” MNRAS, vol. 462, no. 1, pp. 715–725, 2016.
  * [12] L. Hernquist, “An analytical model for spherical galaxies and bulges,” ApJ, vol. 356, pp. 359–364, June 1990.
  * [13] A. Kassiola and I. Kovner, “Elliptic Mass Distributions versus Elliptic Potentials in Gravitational Lenses,” ApJ, vol. 417, p. 450, Nov. 1993\.
  * [14] A. Burkert, “The Structure of Dark Matter Halos in Dwarf Galaxies,” ApJ, vol. 447, p. L25, July 1995.
  * [15] J. Einasto, “On the Construction of a Composite Model for the Galaxy and on the Determination of the System of Galactic Parameters,” Trudy Astrofizicheskogo Instituta Alma-Ata, vol. 5, pp. 87–100, 1965.
  * [16] L. P. Osipkov, “Spherical systems of gravitating bodies with an ellipsoidal velocity distribution,” Soviet Astronomy Letters, vol. 5, pp. 42–44, Jan. 1979.
  * [17] D. Merritt, “Distribution functions for spherical galaxies,” MNRAS, vol. 214, pp. 25P–28P, June 1985.
  * [18] O. Tiret, F. Combes, G. W. Angus, B. Famaey, and H. S. Zhao, “Velocity dispersion around ellipticals in MOND,” A&A, vol. 476, pp. L1–L4, Dec. 2007.
  * [19] A. Biviano, P. Rosati, I. Balestra, A. Mercurio, M. Girardi, M. Nonino, C. Grillo, M. Scodeggio, D. Lemze, D. Kelson, and others, “CLASH-VLT: The mass, velocity-anisotropy, and pseudo-phase-space density profiles of the z = 0.44 galaxy cluster MACS J1206.2-0847,” A&A, vol. 558, p. A1, Oct. 2013.
  * [20] G. A. Mamon, A. Cava, A. Biviano, A. Moretti, B. Poggianti, and D. Bettoni, “Structural and dynamical modeling of WINGS clusters. II. The orbital anisotropies of elliptical, spiral, and lenticular galaxies,” A&A, vol. 631, p. A131, Nov. 2019.
  * [21] A. Cavaliere and R. Fusco-Femiano, “The Distribution of Hot Gas in Clusters of Galaxies,” A&A, vol. 70, p. 677, Nov. 1978.
  * [22] T. Kobayashi, Y. Watanabe, and D. Yamauchi, “Breaking of Vainshtein screening in scalar-tensor theories beyond Horndeski,” Phys. Rev. D, vol. D91, no. 6, p. 064013, 2015.
  * [23] M. Crisostomi and K. Koyama, “Vainshtein mechanism after GW170817,” Phys. Rev. D, vol. D97, no. 2, p. 021301, 2018.
  * [24] A. Dima and F. Vernizzi, “Vainshtein Screening in Scalar-Tensor Theories before and after GW170817: Constraints on Theories beyond Horndeski,” Phys. Rev. D, vol. D97, no. 10, p. 101302, 2018.
  * [25] L. Pizzuti et al., “CLASH-VLT: constraints on $f(R)$ gravity models with galaxy clusters using lensing and kinematic analyses,” JCAP, vol. 1707, no. 07, p. 023, 2017.
  * [26] A. Terukina, L. Lombriser, K. Yamamoto, D. Bacon, K. Koyama, and R. C. Nichol, “Testing chameleon gravity with the Coma cluster,” JCAP, vol. 1404, p. 013, 2014.
  * [27] I. D. Saltas and I. Lopes, “Obtaining Precision Constraints on Modified Gravity with Helioseismology,” Phys. Rev. Lett., vol. 123, no. 9, p. 091103, 2019\.
  * [28] K. Umetsu, A. Zitrin, D. Gruen, J. Merten, M. Donahue, and M. Postman, “CLASH: Joint Analysis of Strong-lensing, Weak-lensing Shear, and Magnification Data for 20 Galaxy Clusters,” ApJ, vol. 821, p. 116, Apr. 2016.
  * [29] A. V. Macciò, A. A. Dutton, and F. C. van den Bosch, “Concentration, spin and shape of dark matter haloes as a function of the cosmological model: WMAP1, WMAP3 and WMAP5 results,” MNRAS, vol. 391, pp. 1940–1954, Dec. 2008.
  * [30] S. H. Hansen and B. Moore, “A universal density slope Velocity anisotropy relation for relaxed structures,” New Astron., vol. 11, pp. 333–338, Mar. 2006.
  * [31] A. Lewis, “GetDist: a Python package for analysing Monte Carlo samples,” 2019\.
  * [32] L. Pizzuti, I. D. Saltas, S. Casas, L. Amendola, and A. Biviano, “Future constraints on the gravitational slip with the mass profiles of galaxy clusters,” MNRAS, vol. 486, no. 1, pp. 596–607, 2019.
  * [33] A. Zaslavski and M. Powell, The NEWUOA software for unconstrained optimization without derivatives, pp. 255–297.  01 2006.