11institutetext: Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
22institutetext: Aix Marseille Univ, Université de Toulon, CNRS, CPT,
Marseille, France

# Towards an optimal marked correlation function analysis for the detection of
modified gravity

M. Kärcher 1122<EMAIL_ADDRESS>J. Bel 22 S. de la Torre 11

Modified gravity (MG) theories have emerged as a promising alternative to
explain the late-time acceleration of the Universe. However, the detection of
MG in observations of the large-scale structure remains challenging due to the
screening mechanisms that obscure any deviations from General Relativity (GR)
in high-density regions. The marked two-point correlation function, which is
particularly sensitive to environment, offers a promising approach to enhance
the discriminating power in clustering analysis and potentially detect MG
signals. This work investigates novel marks based on large-scale environment
estimates but also that exploit the anti-correlation between objects in low-
and high-density regions. This is the first time that the propagation of
discreteness effects in marked correlation functions is investigated in depth,
as in contrast to standard correlation functions, the density-dependent marked
correlation function as estimated from catalogues is affected in a non-trivial
way by shot noise. We assess the performance of various marks to distinguish
GR from MG. This is achieved through the use of the ELEPHANT suite of
simulations, which comprise five realisations of GR and two different MG
theories: $f(R)$ and nDGP. In addition, discreteness effects are thoroughly
studied using the high-density Covmos catalogues. We establish a robust method
to correct for shot-noise effects that can be used in practical analyses. This
methods allows the recovery of the true signal, with an accuracy below $5\%$
over the scales of $5\,h^{-1}\,{\rm Mpc}$ up to $150\,h^{-1}\,{\rm Mpc}$. We
find that such correction is absolutely crucial to measure the amplitude of
the marked correlation function in an unbiased manner. Furthermore, we
demonstrate that marks that anti-correlate objects in low- and high-density
regions are among the most effective in distinguishing between MG and GR, and
uniquely provide visible deviations on large scales, up to about
$80\,h^{-1}\,{\rm Mpc}$. We report differences in the marked correlation
function between $f(R)$ with $|f_{R0}|=10^{-6}$ and GR simulations of the
order of 3-5$\sigma$ in real space. The redshift-space monopole exhibits
similar features and performances as the real-space marked correlation
function. The combination of the proposed $\tanh$-mark with shot-noise
correction paves the way towards an optimal approach for the detection of MG
in current and future galaxy spectroscopic surveys.

###### Key Words.:

large-scale structure of Universe

##  1 Introduction

The seminal works of Riess et al. (1998) and Perlmutter et al. (1999) revived
the cosmological constant $\Lambda$ as a form of dark energy to explain the
late-time accelerated expansion of the Universe. Together with cold dark
matter (CDM), this settled the $\Lambda$CDM model as the current concordance
model of cosmology. Upon closer examination, however, the $\Lambda$CDM model
is found to exhibit certain inherent problems. On the theoretical side, the
fine-tuning problem of $\Lambda$, as extensively studied in Martin (2012),
represents a significant challenge. On the observational side, most recent
cosmological results show a growing tension between early- and late-time
measurements of the Hubble constant $H_{0}$, respectively extracted from the
cosmic microwave background (CMB) anisotropies (Planck Collaboration et al.,
2020) and local distance ladders (Riess et al., 2022). Another source of
contention, although not as significant, comes from an apparent mismatch in
the measured variance of matter fluctuations $\sigma_{8}$, between early-
(Planck Collaboration et al., 2020) and late-time large-scale structure
measurements (Tröster et al., 2020).

In order to address the aforementioned issues, numerous attempts have been
made at the theoretical level. Of particular interest to circumvent the
introduction of a cosmological constant or dark energy component in the first
place, are modified gravity (MG) theories (see Clifton et al., 2012). One of
the most popular modifications is the theory of inflation by Guth (1981)
accommodated with scalar-tensor models to resolve the classical flatness and
horizon problem in $\Lambda$CDM. MG theories are commonly defined and compared
through their respective action. Possibly the most straightforward extension
to the Einstein-Hilbert action in GR consists in replacing the Ricci scalar by
a free function of it, dubbed $f(R)$ theories (see De Felice & Tsujikawa,
2010, for a review). The latter have been subjected to comprehensive analysis,
resulting in tight constraints on viable $f(R)$ functions to ensure the
realisation of accelerated expansion without the necessity of a cosmological
constant, while simultaneously satisfying solar system GR tests (Cognola et
al., 2008). Of particular importance is the $f(R)$ model proposed by Hu &
Sawicki (2007), which simultaneously realises accelerated expansion and evades
solar system tests through the use of a so-called screening mechanism

For MG theories to be a viable replacement or extension to GR they have to
fulfil very stringent tests coming from solar system observations (see
Bertotti et al., 2003; Williams et al., 2004, 2012). On larger scales, the
situation is more complex but there is an increasing effort to tighten
constraints on MG theories by using CMB data (Planck Collaboration et al.,
2016) in combination with large-scale structure and supernovae observables
(Lombriser et al., 2009; Battye et al., 2018). In order for a modification to
standard gravity to shroud its effects on small scales to recover GR,
screening mechanisms are invoked. In general terms, the screening mechanism
describes a suppression of any fifth force to a negligible level such that
gravity follows GR in certain environments. Screening can happen in different
ways such as the chameleon screening in scalar-tensor-theories of gravity (see
Khoury & Weltman, 2004a, b). It also emerges in some $f(R)$ theories due to
the equivalence between scalar-tensor and $f(R)$ theories (Sotiriou & Faraoni,
2010). The DGP gravity model, originally developed by Dvali et al. (2000),
exhibits the screening mechanism first introduced by Vainshtein (1972). A
third popular screening mechanism by Damour & Polyakov (1994) is present in
the symmetron model (Hinterbichler & Khoury, 2010). For a detailed description
of the field of screening mechanisms we refer the reader to Brax et al.
(2022). Intuitively, screening mechanisms in the cosmological context,
particularly the chameleon one, can be understood as a density dependency,
where modification to GR should appear only in low-density regions compared to
the mean density of the Universe. In high-density regions, inside galaxies or
stellar systems for instance, any modification should be negligible. This
imprints a fundamental environmental dependency on the clustering of matter
predicted in those theories.

Since modifications to GR are expected to be small, the observational
detection of MG on cosmological scales poses a major challenge. Guzzo et al.
(2008) advocated the use of the growth rate of structure $f$ measured from
redshift-space distortions (RSD) in the galaxy clustering pattern as an
indicator of the validity of GR in the large-scale structure. Since then, $f$
has become a quantity of major interest and has been measured in large galaxy
redshift surveys (e.g. Blake et al., 2011; Beutler et al., 2012; de la Torre
et al., 2013; Bautista et al., 2021). It is now a standard probe that will be
measured by ongoing surveys, in particular the dark energy spectroscopic
instrument (DESI) (DESI Collaboration et al., 2016) and _Euclid_ mission
(Euclid Collaboration et al., 2024) with an exquisite precision. It is worth
mentioning the $E_{g}$ statistic developed by Zhang et al. (2007), a mixture
of galaxy clustering and weak lensing measurements to probe the properties of
the underlying gravity theory and that has been measured (e.g. Reyes et al.,
2010; de la Torre et al., 2017; Jullo et al., 2019; Blake et al., 2020). Other
quantities that can in principle be measured from observations are the
gravitational slip parameter $\eta$ and the growth index $\gamma$ (see Ishak,
2019, for a review). At the present time, any of the aforementioned
observables has enabled the detection of a deviation from the standard
gravitational field.

In order to improve on existing approaches and to exploit the additional
environmental dependency of MG in clustering analyses, White (2016) proposed
the marked correlation function as a tool to increase the difference in the
clustering signal between MG and GR. In that case, the marked correlation
function is a weighted correlation function normalised to the unweighted
correlation function, and where object weights or marks, are a function of the
local density. The latter is estimated from the density field inferred by dark
matter or galaxies. With this methodology, Hernández-Aguayo et al. (2018),
Armijo et al. (2018), and Alam et al. (2021) investigated marked correlation
functions in N-body simulations of MG. In addition to examining different mark
functions based on density, they also considered marks based on the local
gravitational potential or the host halo mass of the galaxy. They observed
significant differences between MG and GR for marks based on density on small
scales, below about $20\leavevmode\nobreak\ \,h^{-1}\,{\rm Mpc}$. White &
Padmanabhan (2009) also showed the potential of marked correlation functions
to break the degeneracy between of HOD and cosmological parameters. Similar
approaches using weighted statistics or transformation of the density field
have further been proposed. Llinares & McCullagh (2017) used logarithmic
transformations of the density field and computed power spectra of the
transformed field in N-Body simulations to improve on the detection of MG.
Boosting the constraining power in cosmological parameter inference using
power spectra has been shown by using Fisher forecasts by Valogiannis & Bean
(2018), where they compare the Fisher boost using the field transformation of
Llinares & McCullagh (2017), the clipping strategy to mask out high-density
regions (Simpson et al., 2011, 2013) and the mark proposed in White (2016).
Another application of clipping has been done by Lombriser et al. (2015) to
the power spectrum in order to better detect $f(R)$ theories with chameleon
screening. Recently, the use of marked power spectra has been extended to
constrain massive neutrinos (Massara et al., 2021) and tighten constraints on
cosmological parameters (Yang et al., 2020; Xiao et al., 2022).

While a lot of effort on marked statistics for MG has been carried out on
simulation, there have been several applications to observational data.
Satpathy et al. (2019), for the first time, measured marked correlation
functions from observations in the context of MG. They used the proposed
original mark introduced by White (2016) and investigated the monopole and
quadrupole of the marked correlation function measured over the LOWZ sample of
the Sloan Digital Sky Survey (SDSS) DR12 dataset (Alam et al., 2015). They
could not detect significant differences between MG and GR and they attributed
this to modelling uncertainties of the two-point correlation function (2PCF)
on scales of $6\,h^{-1}\,{\rm Mpc}<s<69\,h^{-1}\,{\rm Mpc}$. Armijo et al.
(2024b) applied the strategy introduced in Armijo et al. (2024a) to LOWZ and
CMASS catalogues of SDSS thereby incorporating uncertainties of the HOD on the
projected weighted clustering. They compare predictions from GR and $f(R)$ but
find no significant differences, both fit the LOWZ data and are within the
uncertainties for the investigated scales between $\leavevmode\nobreak\
0.5$$\,h^{-1}\,{\rm Mpc}$ and 40$\,h^{-1}\,{\rm Mpc}$. For the CMASS catalogue
the predictions for both GR and $f(R)$ models fail to properly follow the data
in the first place.

A number of the issues encountered in the literature regarding the use of
marked correlation functions to distinguish MG from GR can be identified as
arising from two main sources. The first is the choice of the mark function,
which, in the majority of cases, results in significant differences on small
scales only. On those scales, a thorough theoretical modelling is difficult as
a proper inclusion of non-linear effects of redshift-space distortions is
needed as well. The second issue is the propagation of discreteness effects in
the mark estimation, i.e. computing the local density from a finite point set,
into the measurement of the marked correlation function. To the best of our
knowledge, this has not been done so far and can lead to biased measurements
if not accounted for. The present work therefore aims at identifying an
optimal mark function that is able to significantly discriminate GR from MG on
larger scales where theoretical modelling is more tractable. For this, we
develop new ways to include environmental information into weighted statistics
as well as investigating new algebraic functions of the density contrast to be
used as a mark. Furthermore, we investigate the discreteness effects and
devise a new methodology to correct marked correlation function measurements
for the bias induced by estimating density-dependent marks on discrete point
sets. We demonstrate that by applying this methodology we are able to robustly
measure the amplitude of marked correlation function and mitigate possible
artefacts in the subsequent analysis of MG signatures.

This article is structured as follows. Section 2 describes the $f(R)$ and nDGP
gravity models that are later investigated and tested. Section 3 introduces
the basics of weighted two-point statistics and marked correlation function.
Section 4 presents the MG simulations used in this work and measurements of
unweighted statistics, which serve as a reference for comparison with the
marked correlation function. Section 5 presents new marks to be used in the
analysis of MG. This is followed in Section 6 by the study of the effects of
shot noise in weighted two-point statistics. Section 7 shows the main results
of this article, which are obtained by applying the previously-defined
methodology to MG simulations. Section 8 comprises a discussion on the optimal
methodology for marked correlation function and conclusions are provided in
Section 9.

##  2 Modified Gravity

We provide in this section a brief review of the theory behind the two classes
of MG models that are used later in this work. In particular, we report the
respective actions alongside with the equation of motion for the additional
scalar degree of freedom, which elucidates the different screening mechanisms
incorporated in those gravity theories.

###  2.1 f(R) Gravity

A general extension to the Einstein-Hilbert action in GR is accomplished by
adding a general function of the Ricci scalar $f(R)$, which then takes the
form

$S=\int\text{d}^{4}x\sqrt{-g}\left\\{\frac{R+f(R)}{16\pi
G}+\mathcal{L}_{m}\right\\},$ (1)

when including a matter Lagrangian $\mathcal{L}_{m}$. This leads to the field
equations

$G_{\alpha\beta}+f_{R}R_{\alpha\beta}-\left(\frac{f}{2}-\Box
f_{R}\right)g_{\alpha\beta}-\nabla_{\alpha}\nabla_{\beta}f_{R}=8\pi
GT_{\alpha\beta},$ (2)

where $\Box$ denotes the d’Alembertian operator and Greek indices are running
from 1 to 4. From these field equations an equation of motion for the scalaron
field $f_{R}=\partial f/\partial_{R}$ can be deduced by taking the trace. The
Ricci scalar is given by

$R=12H^{2}+6HH^{\prime},$ (3)

where a prime denotes a differentiation with respect to the natural logarithm
of the scale factor. In a $\Lambda$CDM universe, today’s Ricci scalar is

$R_{0}=12H_{0}^{2}-9H_{0}^{2}\Omega_{\textrm{m}}^{0}.$ (4)

Although the $f(R)$ function being completely general, there are several
constraints concerning its derivatives with respect to $R$ to obtain a theory
that is free from ghost instabilities (see Tsujikawa, 2010, for a derivation
of those stability conditions). Furthermore, specific functions can be chosen
depending on the context of the theory. Here we focus on a cosmological model
with a late-time accelerated expansion for which the Hu-Sawicki theory (Hu &
Sawicki, 2007) is the most promising. The $f(R)$ function in this model takes
the form

$f(R)=-m^{2}\frac{c_{1}(R/m^{2})^{n}}{c_{2}(R/m^{2})^{n}+1},$ (5)

with $m=8\pi G\rho_{0}/3$, and $c_{1}$, $c_{2}$, $n$ being constants. In the
simulations presented in the next sections, a value of $n=1$ was used. To
produce a background expansion as dictated by $\Lambda$CDM, the ratio between
$c_{1}$ and $c_{2}$ has to be chosen such that

$\frac{c_{1}}{c_{2}}=6\frac{\Omega_{\Lambda}^{0}}{\Omega_{\textrm{m}}^{0}}.$
(6)

From this follows a Lagrangian of the form $\mathcal{L}=R/16\pi G-\Lambda$ for
the gravitational sector in the $R>>m^{2}$ limit where
$f(R)\approx-m^{2}c_{1}/c_{2}$, which corresponds to the well-known Einstein-
Hilbert action with cosmological constant. Furthermore, by expanding the
$f(R)$ function in the aforementioned limit but keeping the next-to-leading
order term we arrive at

$f(R)=-\frac{c_{1}}{c_{2}}m^{2}\left(1-\frac{m^{2}}{Rc_{2}}\right)=-m^{2}6\frac{\Omega_{\Lambda}}{\Omega_{\textrm{m}}}-f_{R0}\frac{R_{0}^{2}}{R},$
(7)

where in the second equality we used the expression of the scalaron field

$f_{R}=-\frac{m^{4}}{R^{2}}\frac{c_{1}}{c_{2}^{2}}$ (8)

evaluated for the background Ricci scalar value today ($f_{R0}$). We replaced
$c_{1}/c_{2}$ with the previous expression to obtain a $\Lambda$CDM
background. In this approximation, and by fixing $n=1$, the $f(R)$ function
depends solely on the cosmological parameters and $f_{R0}$, the latter
encoding the strength of the modification to GR.

Having an $f(R)$ modification in the Lagrangian will introduce additional
force terms into the Poisson equation in the quasi-static and weak-field
limit, as can be derived from perturbed field equations (Bose et al., 2015)

$\nabla^{2}\Phi=4a^{2}\pi G(\rho-\bar{\rho})-\frac{1}{2}\nabla^{2}f_{R}$ (9)

and

$\nabla^{2}f_{R}=-\frac{a^{2}}{3}(R-\bar{R})-\frac{8\pi
G}{3}a^{2}(\rho-\bar{\rho}),$ (10)

where $\bar{\rho}$ and $\bar{R}$ are the matter density and Ricci scalar at
the background level. These additional terms should be suppressed in the
vicinity of massive objects, otherwise solar system tests might have detected
the fifth force. When $f(R)$ gravity is rewritten as a scalar-tensor gravity,
the potential of the scalar field receives a contribution from the matter
density (Khoury & Weltman, 2004a) as

$V_{\textrm{eff}}(\varphi)\equiv V(\varphi)+\rho e^{\varphi\beta/M_{pl}}.$
(11)

This leads in turn to a modified equation of motion for the scalar field
$\varphi$ that includes density-dependent potential. In this context a thin-
shell condition can be derived, stating that the difference between the scalar
field far away from the source $\varphi_{\infty}$ and inside the object
$\varphi_{c}$ should be small compared to the gravitational potential on the
surface of the object (Khoury & Weltman, 2004a). Exterior solutions for
$\varphi$ around compact objects satisfying the thin-shell condition will
reach the solution $\varphi_{\infty}$ at larger distances, thereby suppressing
the effect of the scalar field close to the object.

###  2.2 nDGP Gravity

The modification to standard gravity devised by Dvali, Gabadadze and Porrati
(Dvali et al., 2000), hereafter DGP gravity, is of a radically different kind
compared to $f(R)$ gravity. The setup is a 4D brane embedded in a 5D bulk and
the modification to gravity comes from the fith dimensional contribution. The
action is given by (Clifton et al., 2012)

$\begin{split}S=&M_{5}^{3}\int\text{d}^{5}x\sqrt{-g_{5}}\,R_{5}\\\
&+\int\text{d}^{4}x\sqrt{-g_{4}}\left\\{-2M_{5}^{3}K+\frac{M_{4}^{2}}{2}R_{4}-\sigma+\mathcal{L}_{m}\right\\},\end{split}$
(12)

where $g_{5}$ and $g_{4}$ are the 5D and 4D metric, respectively. The matter
Lagrangian $\mathcal{L}_{m}$ does live on the 4D brane as well as the brane
tension $\sigma$, which can act as a cosmological constant. Furthermore, there
is both a 5D Ricci scalar $R_{5}$ and its 4D counterpart $R_{4}$, and the
brane has an extrinsic curvature term $K$. Generally, both the brane and bulk
have their individual mass scales $M_{4}$ and $M_{5}$ and they give rise to a
specific cross-over scale $r_{c}$ defined as

$r_{c}=\frac{M_{4}^{2}}{2M_{5}^{3}},$ (13)

which regulates the contribution of 4D with respect to 5D gravity.

The modified Poisson equation for the gravitational potential and the equation
for the additional scalar degree of freedom $\varphi$ (also called brane-
bending mode as it describes the displacement of the brane) lead to the fifth
force. They are given in the quasi-static approximation by (see Koyama &
Silva, 2007)

$\begin{split}&\nabla^{2}\Phi=4\pi
Ga^{2}(\rho-\bar{\rho})+\frac{1}{2}\nabla^{2}\varphi\\\
&\nabla^{2}\varphi+\frac{r_{c}}{3\beta
a^{2}}\left((\nabla^{2}\varphi)^{2}-(\nabla_{i}\nabla_{j}\varphi)^{2}\right)=\frac{8\pi
Ga^{2}}{3\beta}(\rho-\bar{\rho}),\end{split}$ (14)

where $\beta$ is

$\beta(t)=1\pm 2Hr_{c}\left(1+\frac{\dot{H}}{3H^{2}}\right).$ (15)

The dot refers to a derivative with respect to metric time $t$. One important
feature of the DGP model is the existence of a normal branch and of a self-
accelerating branch, indicated respectively by the $+$ and $-$ signs in the
equation for $\beta$. While the self-accelerating branch appears appealing for
cosmology at first sight, as it can generate accelerated expansion without
cosmological constant (the limit of vanishing brane tension), it contains
unphysical ghost instabilities (Clifton et al., 2012). Hence the model used in
the simulation analysed in this work implements the normal branch, which does
need a non-vanishing brane tension to produce accelerated expansion. It is
interesting to study normal branch DGP models as it exhibits the Vainshtein
screening mechanism (see Schmidt, 2009; Barreira et al., 2015). To illustrate
that mechanism, the equation for the scalar field has to be studied around a
mass source. Far away from the source, only the linear term
$\nabla^{2}\varphi$ will dominate and this will contribute substantially to
the usual gravitational force as it will also scale $\propto 1/r$. However,
non-linear terms start to dominate once we are closer to the source than to
the Vainshtein radius $r_{V}$, defined by

$r_{V}\approx(r_{s}r_{c}^{2})^{1/3},$ (16)

with $r_{s}$ being the Schwarzschild radius of the source. At some point, non-
linear terms will dominate and the resulting force will scale as $\sqrt{r}$
and hence will be suppressed with respect to the gravitational force. A
derivation of the solution for $\varphi$ can be found in Koyama & Silva (2007)
for the general case, which includes linear and non-linear terms, and where
the same scaling are recovered in the respective regimes.

At fixed Schwarzschild radius, the cross-over scale determines the Vainshtein
radius, so by running simulations with different $r_{c}$ one will obtain
different strengths of the Vainshtein screening. Therefore, varying $r_{c}$
allows the tuning of the amount of deviation to GR that is required.

##  3 Weighted statistics and estimators

In Tab. 1 we summarise the notation used throughout this work to ease
distinguishing between the different discrete and continuous quantities.

$\langle\rangle$ | ensemble average (moment)
---|---
$\langle\rangle_{c}$ | ensemble cumulant
$\delta_{M}(\mathbf{x})$ | continuous weighted density contrast
$\delta(\mathbf{x})$ | continuous density contrast
$\delta_{f}(\mathbf{x})$ | discrete density contrast
$\delta_{R}(\mathbf{x})$ | continuous smoothed density contrast
$\delta_{Rf}(\mathbf{x})$ | discrete smoothed density contrast
$m(\mathbf{x})$ | mark field
$W(\mathbf{r})$ | weighted correlation function
$W_{f}(\mathbf{r})$ | estimated weighted correlation function
$\mathcal{M}(\mathbf{r})$ | marked correlation function
$\mathcal{M}_{f}(\mathbf{r})$ | estimated marked correlation function
$F(\mathbf{r})$ | smoothing kernel
$\bar{m}$ | mean mark
$\bar{m}_{f}$ | mean mark taken over a point set
$\bar{n}$ | mean number of points per volume
$\bar{N}$ | mean number of points per grid cell
$a$ | size of one grid cell
Table 1: Summary of notation used in this article.

###  3.1 Unweighted Statistics

The density contrast $\delta(\mathbf{x})$, which encodes the relative change
of the density field $\rho(\mathbf{x})$, is defined as

$\delta(\mathbf{x})=\frac{\rho(\mathbf{x})-\bar{\rho}}{\bar{\rho}},$ (17)

where $\bar{\rho}$ is the mean density. In order to study the matter
clustering in the cosmological context, one of the most common summary
statistic to characterise the density field, is the two-point correlation
function $\xi(\mathbf{x},\mathbf{y})$ or its Fourier counterpart the power
spectrum. The 2PCF is the cumulant
$\langle\delta(\mathbf{x})\delta(\mathbf{y})\rangle_{c}$ of the density
contrast at positions $\mathbf{x}$ and $\mathbf{y}$. For two-point
correlations, the cumulant and standard ensemble average are the same
quantity. They remain the same up to three-point correlations but start to
differ from four-point correlations onwards. Due to the assumed statistical
invariance by translation, the correlation function does only depend on the
separation vector $\mathbf{r}=\mathbf{x}-\mathbf{y}$. By inserting the
definition of the density contrast we have that

$\xi(\mathbf{r})=\frac{\langle\rho(\mathbf{x}+\mathbf{r})\rho(\mathbf{x})\rangle_{c}}{\bar{\rho}^{2}}.$
(18)

From the last equation one can see that the 2PCF is zero if the field is
totally uncorrelated at two different positions.

In order to estimate the 2PCF, we can deploy the commonly-used Landy-Szalay
pair-counting estimator proposed by Landy & Szalay (1993) to minimise the
variance, and which takes the form

$\xi(\mathbf{r})=\frac{DD(\mathbf{r})-2DR(\mathbf{r})+RR(\mathbf{r})}{RR(\mathbf{r})}.$
(19)

The terms $DD(\bf{r})$ and $RR(\bf{r})$ are the normalised pair counts
measured in the data sample and a random sample following the geometry of the
data sample, respectively. In addition, a cross term with pairs consisting of
one point in the data sample and the other in the random sample is given by
$DR(\bf{r})$. In this work, we only compute two-point correlation functions in
periodic boxes without selection function. In this case, the term $DR$
converges to the term $RR$ in the limit of many realisations of random
catalogues, and we can use the natural estimator given by Peebles & Hauser
(1974)

$\xi(r)=\frac{DD(r)-RR(r)}{RR(r)}.$ (20)

The distribution of pairs in real space is isotropic, and together with
periodic boundary conditions, lead the correlation function to only depend on
the modulus $r$ of the pair separation vector.

In redshift space, it is useful to compute the anisotropic correlation
function $\xi(s,\mu)$, binned in the norm of the pair separation vector
$\mathbf{s}$ and the cosine angle between the line of sight (LOS) and the pair
separation vector $\mu$. The 2PCF estimator for a periodic box is hence

$\xi(s,\mu)=\frac{DD(s,\mu)-RR(s,\mu)}{RR(s,\mu)},$ (21)

and normalised $RR$ counts are given by

$\begin{split}RR([[s,s+\Delta
s],[\mu+\Delta\mu]])&=\frac{1}{L^{3}}\frac{4}{3}\pi\Delta\mu\\{(s+\Delta
s)^{3}-s^{3}\\},\end{split}$ (22)

which can be derived by calculating the volume covered by the respective bins
in $s$ and $\mu$ relative to the total volume of the bin. For real space
measurements, the $RR$ counts can be evaluated analytically in a similar
fashion as in Eq. (22).

The 2PCF in redshift space, $\xi(s,\mu)$, can be decomposed into multipole
moments, which is a basis encoding the different angle dependencies of the
full 2PCF. Usually the decomposition is done into the first three non-
vanishing multipole moments, being the monopole, quadrupole and hexadecpole.
In the following, we will focus on the first two since the hexadecapole can be
quite noisy for small point sets. The multipole moment correlation functions
are obtained by decomposing the $\xi(s,\mu)$ in the basis of Legendre
polynomials as

$\xi_{\ell}(s)=\frac{(2\ell+1)}{2}\int_{-1}^{1}\text{d}\mu\,\xi(s,\mu)P_{\ell}(\mu),$
(23)

yielding for the monopole and quadrupole to

$\begin{split}\xi_{0}(s)&=\frac{1}{2}\int_{-1}^{1}\text{d}\mu\,\xi(s,\mu)\\\
\xi_{2}(s)&=\frac{5}{2}\int_{-1}^{1}\text{d}\mu\,\xi(s,\mu)\frac{1}{2}(3\mu^{2}-1).\end{split}$
(24)

In practice, these integrals are discretised and we measure $\xi(s,\mu)$ in
$100$ bins from $\mu=0$ to $\mu=1$ using the symmetry under interchange of
galaxies for a given pair, which is fulfilled in our periodic box simulations.
The discretised correlation function is then integrated by approximating the
integral as a Riemann sum.

###  3.2 Weighted statistics

Let us now define the weighted density contrast

$\delta_{M}(\mathbf{x})=\frac{\rho_{M}(\mathbf{x})-\overline{\rho_{M}}}{\overline{\rho_{M}}},$
(25)

where the weighted density field is given by
$\rho_{M}(\mathbf{x})=m(\mathbf{x})\rho(\mathbf{x})$,
$\overline{\rho_{M}}=\langle\rho_{M}(\mathbf{x})\rangle$, and $m(\bf{x})$ is
the mark field. The weighted correlation function is the ensemble average of
the weighted density contrast correlation,

$W(\mathbf{r})=\langle\delta_{M}(\mathbf{x})\delta_{M}(\mathbf{x}+\mathbf{r})\rangle,$
(26)

which when substituted with the definition of the density contrast, takes the
form

$\begin{split}1+W(\mathbf{r})&=\frac{\bar{\rho}^{2}}{\overline{\rho_{M}}^{2}}\langle
m(\mathbf{x})(1+\delta(\mathbf{x}))m(\mathbf{x}+\mathbf{r})(1+\delta(\mathbf{x}+\mathbf{r}))\rangle.\end{split}$
(27)

The mark field $m(\mathbf{x})$ can be continuous in space, or discrete and
defined on the point set (galaxy or halo catalogue). Each object in the
catalogue can be assigned a mark from the mark field, e.g. the $i$-th object
has a mark $m_{i}$. The normalised weighted pair counts are obtained as

$WW(r)=\frac{\sum_{i\neq j}m_{i}m_{j}}{(\sum m_{i})^{2}-\sum m_{i}^{2}},$ (28)

where the sum is computed over all pairs with a separation inside the bin
centred on $r$. The marked correlation function is then defined as (Beisbart &
Kerscher, 2000; Sheth, 2005)

$\mathcal{M}(r)\equiv\frac{1+W(r)}{1+\xi(r)}.$ (29)

It converges to $\mathcal{M}(r)=1$ on large scales as $W(r)$ and $\xi(r)$
approach zero.

In order to estimate the weighted correlation function from a catalogue, the
natural estimator can be generalised to include weighted $DD(r)$ so that we
can simply replace $DD(r)$ with $WW(r)$ counts, arriving at

$W(r)=\frac{WW(r)-RR(r)}{RR(r)}.$ (30)

Inserting this into the definition of the marked correlation function
$\mathcal{M}(r)$ we have

$\mathcal{M}(r)=\frac{1+\frac{WW(r)-RR(r)}{RR(r)}}{1+\frac{DD(r)-RR(r)}{RR(r)}}=\frac{WW(r)}{DD(r)}.$
(31)

If the LS estimator is employed instead, one has to compute $WR(r)$ and
$DR(r)$ terms in addition.

Computing the multipoles of the weighted correlation function is analogous to
the unweighted case. However, the multipoles of the marked correlation
function can be defined in two ways. The most intuitive definition is obtained
by decomposing the marked correlation function $\mathcal{M}(s,\mu)$ in the
basis of Legendre polynomials, yielding

$\mathcal{M}_{\ell}(s)=\frac{(2\ell+1)}{2}\int_{-1}^{1}\text{d}\mu\,\mathcal{M}(s,\mu)P_{\ell}(\mu).$
(32)

The second approach uses the following definition

$\mathcal{M}_{\ell}(s)=\frac{1+W_{\ell}(s)}{1+\xi_{\ell}(s)},$ (33)

which is motivated by the fact that the denominator is the actual multipole of
the unweighted 2PCF. This is not the case in the first definition in Eq. (32).
The second definition has been used for instance by White (2016) and Satpathy
et al. (2019). Throughout this work we will use the form of Eq. (32).

##  4 Simulations

###  4.1 Characteristics

In order to investigate different marked correlation functions and assess
their discriminating power regarding GR and MG, we use the Extended LEnsing
PHysics using ANalytic ray Tracing (ELEPHANT) simulation suite, thoroughly
discussed in Sec. II B. of Alam et al. (2021). We only provide a brief
description of it in the following. This simulation suite consists of 5
realisations of GR with $\Lambda$CDM cosmology, $f(R)$ gravity with three
different values of $|f_{R0}|=[10^{-6},10^{-5},10^{-4}]$, and nDGP gravity
with $H_{0}r_{c}=[5.0,1.0]$. Henceforth, we will refer to the different
simulations as GR, F6, F5, F4, N5, and N1, respectively. The background
cosmology is summarised in Tab. 2 and resembles the best-fitting cosmology
obtained from 9-year WMAP CMB analysis presented in Hinshaw et al. (2013).

| ELEPHANT | DEMNUni/Covmos
---|---|---
$\Omega_{b}$ | 0.046 | 0.05
$\Omega_{cdm}$ | 0.235 | 0.27
$\Omega_{m}$ | 0.281 | 0.32
$\Omega_{\Lambda}$ | 0.719 | 0.68
$h$ | 0.697 | 0.67
$n_{s}$ | 0.971 | 0.96
$\sigma_{8}$ | 0.82 | -
$A_{s}$ | - | $2.1265\times 10^{-9}$
Table 2: Reference cosmology of the ELEPHANT (first column) and DEMNUni
(second column) simulations. The PDF of the dark matter particles in the
DEMNUni simulation have been used to define the target PDF used to produce
Covmos realisations.

Key simulation parameters are summarised in Tab. 3. The dark matter halos have
been identified with the ROCKSTAR algorithm (Behroozi et al., 2013) and have
been subsequently populated with galaxies using the 5-parameter halo
occupation distribution (HOD) model of Zheng et al. (2007). For each
realisation, redshift-space coordinates have been calculated by fixing the LOS
to one of the three simulation box axes, individually, and by ’observing’ the
box from a distance equal to 100 times the box side length.

N-body code | ECOSMOG
---|---
Box side length $L$ | 1024$\,h^{-1}\,{\rm Mpc}$
Particle number $N_{p}$ | $1024^{3}$
Number of grid cells $N_{grid}$ | $1024^{3}$
Initial redshift $z_{ini}$ | 49.0
Initial conditions | Zel’dovich approx. (MPGrafic)
Mass particle $M_{p}$ | $7.79854\times 10^{10}\textrm{M}_{\odot}$
Final redshift z | 0.506
Table 3: Characteristics of the ELEPHANT simulation suite.

One crucial property of this suite of simulations, which makes it particularly
suitable for our studies, is the matching of the projected 2PCF of galaxies
$w_{p}(r_{p})$ predicted by GR in the MG simulations. The latter was done by
tuning the HOD parameters of the MG simulations. For the GR simulation, the
best-fit HOD parameters were taken from Manera et al. (2013).

We use a second set of simulations to assess discreteness effects in the
estimation of the mark and how they propagate into the marked correlation
function. For this, we make use of the Covmos realisations from Baratta et al.
(2023). These are not full N-body simulations, rather they reproduce dark-
matter particle one- and two-point statistics following the technique
described in Baratta et al. (2020). This procedure consists of applying a
local transformation to a Gaussian density field such that it follows a target
probability distribution function (PDF) and power spectrum. The point set is
then obtained by a local Poisson sampling on the linearly interpolated density
values. For the set of Covmos realisations, the target PDF and power spectrum
were set by the DEMNUni N-body simulation (Castorina et al., 2015) statistics.
The DEMNUni simulation assumes a $\Lambda$CDM cosmology with parameters
presented in Tab. 2. The Covmos catalogues contain about $20\times 10^{6}$
points in a box of 1 $h^{-1}$Gpc of side, resulting in a number density of
about $0.02\,h^{3}\leavevmode\nobreak\ \mathrm{Mpc}^{-3}$. Such an high
density allows treating those catalogues as being almost free from shot noise.

###  4.2 Two-point correlation function

Although the galaxy projected correlation function are matched in the ELEPHANT
suite, it is instructive to assess residual deviations in other statistics,
particularly for the interpretation of differences arising in the analysis of
marked correlation functions.

We measured both the real- and redshift-space correlation functions in 30
linear bins in $r$ and $s$, respectively, ranging from
$10^{-3}\leavevmode\nobreak\ \,h^{-1}\,{\rm Mpc}$ to $150\leavevmode\nobreak\
\,h^{-1}\,{\rm Mpc}$. For the redshift-space measurements, we used the
ELEPHANT catalogues with the LOS fixed to the $x$-direction. All correlation
function measurements in this work have been performed using the publicly
available package Corrfunc (Sinha & Garrison, 2019, 2020). In the upper panel
of Figure 1, we show the standard correlation function in real space for the
different gravity simulations. The measurements appear to be within the
respective uncertainties over all scales, albeit on very small scales, a more
careful assessment of possible deviations is advised as the error bars are
very small. In Figure 2, the monopole (right) and quadrupole (left) of the
anisotropic 2PCF in redshift space are presented in the upper panels.
Similarly as for the real space correlation function, the multipoles are
within the respective uncertainties on large scales, although the N1
measurement appears to deviate from the others in the quadrupole. On smaller
scales, discrepancies seem to appear as uncertainties are getting very small
and a visual inspection is not sufficient to quantify those differences.

In order to properly assess the difference between MG and GR in weighted or
unweighted correlation functions, we define the difference between MG and GR
as the mean

$\overline{\Delta\mathcal{M}}(r)=\frac{1}{5}\sum_{i=1}^{5}\\{\mathcal{M}_{i,\text{MG}}(r)-\mathcal{M}_{i,\text{GR}}(r)\\},$
(34)

where $i$ ranges over the number of realisations. However, the mean
differences alone does not tell about the significance as the data might
fluctuate much more than differences. We therefore divide the mean difference
by the standard deviation as

$\sigma_{\textrm{avg}}(r)=\sqrt{\frac{1}{N(N-1)}\sum_{i=1}^{5}\\{\Delta_{i}\mathcal{M}(r)-\overline{\Delta\mathcal{M}}(r)\\}^{2}}.$
(35)

The factor of $1/(N-1)$ is necessary in order to compute an unbiased standard
deviation, since we only have 5 realisations at hand. Furthermore, the
additional factor of $1/N$ comes from fact that we want the error on the mean
and not of a single measurement. In a similar manner, we compute the standard
deviation of a single marked correlation function as

$\sigma_{\textrm{s}}(r)=\sqrt{\frac{1}{N-1}\sum_{i=1}^{5}\\{\mathcal{M}_{i}(r)-\overline{\mathcal{M}}(r)\\}^{2}}.$
(36)

In the end, the ratio of interest is

$\textrm{SNR}(r)=\frac{\overline{\Delta\mathcal{M}}(r)}{\sigma_{\textrm{avg}}(r)},$
(37)

giving directly the difference in terms of standard deviations. If the
absolute value of this SNR is larger than 3 then we would advocate a
significant deviation between MG and GR.

Another quantity of interest that we use throughout this work is the ratio
between the error on a single measurement of the marked correlation function
and the noise $\sigma_{\textrm{avg}}$, as used in the signal-to-noise ratio.
We will refer to this ratio as

$\alpha(r)=\frac{\sigma_{\textrm{s}}(r)}{\sigma_{\textrm{avg}}(r)}$ (38)

and will include it the figures as shaded regions. The error of a single
measurement $\sigma_{s}(r)$ is hereby taken for the GR case. This ratio gives
an indication on the statistical significance of a difference, if we would
have only one simulation/measurement at hand. To assess this, we have to
compare $\textrm{SNR}(r)$ with $\alpha(r)$, and if
$\textrm{SNR}(r)>3\alpha(r)$ then we can claim a $3\sigma$ difference to be
detectable with a single measurement. Of course, care must be taken if the
error of a single measurement is significantly different between GR and MG,
since $\alpha(r)$ will differ depending on what simulations are used to
estimate the error, therefore possibly affecting conclusions.

In the lower panel of Fig. 1, we display the SNR$(r)$ as introduced above. The
differences rarely cross the limit of $3\sigma$ except for the very lowest
scales, below 20$\,h^{-1}\,{\rm Mpc}$, or for F4 at intermediate scales where
deviations can reach up to 6$\sigma$. However these large deviations happen
only for single scales and there is no general trend. This suggests that the
crossing of the 3$\sigma$ border might be caused by sample variance, as the
statistical power is limited with only 5 realisations. In addition, when
considering the error of a single realisation volume, as displayed by the blue
shaded region, we can see that the deviations for F4 are within the
uncertainty.

In the lower panels of Fig. 2 we show the SNR$(s)$ for the multipoles in
redshift space. Except for the smallest scales, the simulations F5, F6 and N5
show generally no significant differences to GR. In the case of the monopoles,
the SNR is close to 0 for almost all scales. For F4 and N1 significant
differences are present although mainly on scales below
$20\leavevmode\nobreak\ \,h^{-1}\,{\rm Mpc}$. On larger scales both
simulations show SNR varying around 3$\sigma$. These differences are mostly
within the uncertainty, if the error of a single volume is considered. This
confirms that by considering the standard correlation function only, in real
or redshift space, we cannot really distinguish between GR and those MG
models.

Figure 1: Difference in the measured standard correlation function $\xi(r)$
between GR and MG in real space. In the upper panel the correlation functions
themselves are plotted where different colours indicate the underlying gravity
theory. The curves show the average over 5 realisations and the errorbars
correspond to the mean standard deviation over these realisations. The lower
panel quantifies possible differences in terms of the SNR as introduced in
Section 4. Black dashed lines indicate a SNR of $\pm 3$. The shaded region
refers to the error of a single measurement divided by the mean error of the
difference as described in Eq. (38).

Figure 2: Differences in the standard measured correlation function multipoles
in redshift space between GR and MG. The upper panels present the mean
correlation function multipoles taken over 5 realisations with the monopole on
the left side and the quadrupole on the right side. The errorbars correspond
to the mean standard deviation over 5 realisations. The lower panels show the
respective SNR with 3$\sigma$ indicated by the black dashed lines. The colour
coding refers to different gravity simulations and the corresponding shaded
regions refer to the error of a single measurement divided by the mean error
of the difference (see Eq. (38)).

##  5 Marks for modified gravity

There is a vast space of possible marks that can be used and the specific
choice strongly depends on the context in which the marked correlation
function is studied. The most popular mark function $M[\rho(\bf{x})]$ in the
literature within the context of detecting MG was introduced by White (2016)
and takes the form

$m(\mathbf{x})=M_{W}[\rho(\mathbf{x})]\equiv\left(\frac{1+\rho_{*}}{\rho_{*}+\rho(\mathbf{x})}\right)^{p},$
(39)

where $\rho_{*}$ and $p$ are free parameters used to control the mark
upweighting of low- versus high-density regions. We will refer to this mark in
the following as the White mark and be indicated via the $W$ in the subscript.
The White mark can be seen as a _local_ transformation of the density field.
Choices for the free parameters range from $(\rho_{*},p)=(10.0,7.0)$ (Aviles
et al., 2020), $(\rho_{*},p)=(4.0,10.0)$ Alam et al. (2021); Valogiannis &
Bean (2018), to $(\rho_{*},p)=(10^{-6},1.0)$ (Hernández-Aguayo et al., 2018).
Upweighting galaxies in high-density regions using $(\rho_{*},p)=(1.0,-1.0)$
has also been explored in Alam et al. (2021). Other values have also been
investigated by Satpathy et al. (2019) and Massara et al. (2021) for the
marked power spectrum. This underlines the wide range of possible mark
functions and configurations to be used and the amount of freedom this can
introduce in the analysis.

Marks based on the local density require an estimation of the latter from a
finite point set in the first place and there exist several different
approaches to do so. While we will use an estimation based on mass assignment
schemes (MAS), adaptive approaches such as Delaunay (Schaap & van de Weygaert,
2000) or Voronoi tessellations as used in void finders (e.g. Neyrinck, 2008)
could also be used. With a MAS applied to a discrete density field (subscript
$f$ for finite)

$\rho_{f}(\mathbf{x})=m\sum_{i=0}^{N-1}\delta_{D}(\mathbf{x}-\mathbf{x}_{i}),$
(40)

where the $i$-th point is located at position $\mathbf{x}_{i}$, the estimated
density field on the grid takes the form (Sefusatti et al., 2016)

$\delta_{Rf}(\mathbf{x})=\frac{1}{\bar{N}}\sum_{i=0}^{N-1}F\left(\frac{\mathbf{x}-\mathbf{x}_{i}}{a}\right)-1,$
(41)

with $\bar{N}$ is the density of points per grid cell, $N$ is the number of
points, $a$ is the size of one grid-cell, and $F(\mathbf{x})$ the MAS kernel.
The coordinate $\mathbf{x}$ is only evaluated at grid points but in principle
can be placed anywhere. For this derivation, we assumed all points to have the
same mass $m$ and we made use of the simplified notation where
$F(\mathbf{x})\equiv F(x_{1})F(x_{2})F(x_{3})$. The density field obtained in
this way is related to the true density field by a convolution with the MAS
kernel. In this work, we mainly use a piece-wise cubic spline (PCS) for the
MAS but higher- and lower-order kernels are employed for specific tests. The
explicit form of the used kernels up to septic order can be found in the
appendix of Chaniotis & Poulikakos (2004)111It has to be noted that there are
two minor typos. Their octic spline is actually the septic spline and in its
expression the term $s^{7}/20$ should be replaced by $s^{7}/720$..

###  5.1 Beyond local density

A way to include information beyond the local density field is by using the
large-scale environment that can be divided into clusters, filaments, walls
and voids. Generally, there are different ways to define these structures from
a galaxy catalogue ranging from the sophisticated approach by Sousbie (2011)
based on topological considerations to the work of Falck et al. (2012) using
phase-space information. One of the most straightforward approaches utilises
the T-web formalism (Forero-Romero et al., 2009) based on the Hessian of the
gravitational potential. For a thorough comparison of the above mentioned
cosmic web classifications and many more we refer the reader to Libeskind et
al. (2018). In this analysis we will deploy the T-web classification that uses
the relation of the eigenvalues $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$
of the tidal tensor to the density evolution as given in Cautun et al. (2014)

$\rho(\mathbf{x})=\frac{\bar{\rho}}{(1-D(t)\lambda_{1})(1-D(t)\lambda_{2})(1-D(t)\lambda_{3})},$
(42)

where $D(t)$ is the growing solution to the growth factor. This expression can
be derived from Lagrangian perturbation theory to linear order (Zel’dovich,
1970). The dimensionality of the structure then depends on the number of
eigenvalues with positive sign. Three positive eigenvalues corresponds to a
cluster as it encodes a collapse among all three spatial directions. Two or
one positive eigenvalues result in a filament or wall, respectively. If all
eigenvalue are negative then $\rho(\mathbf{x})$ will never diverge and we can
interpret this as a void. A pitfall of this classification appears if some of
the eigenvalues are very small but positive as the corresponding structure
might not collapse in an Hubble time. To circumvent this issue while not
having to rely on thresholds for the eigenvalues we use the scheme as proposed
in Cautun et al. (2013). They give an environmental signature $\mathcal{S}$
for ordered eigenvalues $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ in their
Eq. 6 and 7 as signatures for clusters, filaments, walls and voids,
respectively. We adapted this scheme to be used on the eigenvalues of the
Tidal tensor instead of the Hessian of the density contrast as proposed in
Cautun et al. (2013).

In order to obtain the tidal tensor in a simulation we follow the grid-based
approach as used for the density field. With a density field on a grid at
hand, the gravitational potential or tidal tensor can be straightforwardly
deduced by a series of fast Fourier transforms (FT). For this we use the
Poisson equation to relate the density field to the gravitational potential

$\nabla^{2}\Phi(\mathbf{x})=4\pi
G\delta(\mathbf{x})\quad\overset{\textrm{FT}}{\Leftrightarrow}\quad
k^{2}\Phi(\mathbf{k})=-4\pi G\delta(\mathbf{k}).$ (43)

We absorb the constant $4\pi G$ into the definition of the gravitational
potential. There exists a singularity when the wavevector is equal to zero,
which we evade by simply setting the zeroth mode of $\Phi(\mathbf{k})$ to
zero, as we expect the gravitational potential sourced by the density contrast
to have a zero mean. The components of the tidal tensor $T_{ij}$ can then be
derived by taking successive derivatives in the respective directions as

$T_{ij}(\mathbf{x})=\partial_{i}\partial_{j}\Phi(\mathbf{x})\quad\overset{\textrm{FT}}{\Leftrightarrow}\quad
T_{ij}(\mathbf{k})=-k_{i}k_{j}\Phi(\mathbf{k}).$ (44)

The off-diagonal terms suffer from a break in the Fourier symmetry pairs when
evaluated on a finite grid. This leads to a non-vanishing imaginary part once
the tidal tensor in configuration space is obtained by inverse FT. We
circumvent this issue by setting the imaginary part to zero via applying a
filter that sets the symmetry breaking modes at the Nyquist frequency to zero.
In our implementation we will evaluate the environmental signatures on the
grid over which the eigenvalues of the tidal tensor have been computed. For
each grid cell the largest signatures defines the corresponding environment
and if all signatures are zero then the environment is set to be a void.

Once each galaxy has been classified, this information can be used to enhance
effects of MG in clustering measurements. The simplest use of the
environmental classification is to divide the catalogue into sub-catalogues
consisting of galaxies located in voids, walls, filaments and clusters,
respectively and compute auto-correlation functions. The difference in
clustering amplitude in the different environments is expected to be stronger
in MG, particularly in the correlation function of void galaxies. In the work
of Bonnaire et al. (2022), they computed the power spectra of density fields
that have been obtained by splitting the original density field into
respective contributions from galaxies in voids, walls, filaments and
clusters. They applied this approach to the dark matter particles of the
Quijote simulation (Villaescusa-Navarro et al., 2020) thereby having a much
larger set of points. A drawback in the analysis of our galaxy mock
catalogues, but also when using limited survey samples, is the loss of
information due to discarding many galaxies leading to an increase in
uncertainty of the measurements. Computing environmental correlation functions
is the same as computing weighted correlation functions but with a mark set to
one for all galaxies living in the respective environment and zero otherwise.
We refer the interested reader to Section A for some notes on the structure of
weighted correlation functions for those kind of weights. In the following we
will denote such marked correlation functions, consisting of the environmental
weighted correlation function divided by the total unweighted correlation
function as in Eq. (29), via their respective environment, e.g. _void marked
correlation function_.

Conversely, we can use the full catalogue of objects and put larger weights to
progressively more unscreened galaxies, as done by the following mark field

$m(\mathbf{x})=\begin{cases}4\quad\text{if void}\\\ 3\quad\text{if wall}\\\
2\quad\text{if filament}\\\ 1\quad\text{if
cluster}\end{cases}\quad\text{Void${}_{\rm LEM}$},$ (45)

where LEM stands for _linear environment mark_. This approach is similar to
the density split technique as used in Paillas et al. (2021), where cross-
correlation functions of galaxies living in differently dense regions have
been investigated. Our proposed mark can also be devided into a specific
combination of auto and cross-correlations. This mark could in principle be
extended to a $\textrm{Wall}_{\textrm{LEM}}$ mark, where wall galaxies get
assigned a weight of 4 and voids galaxies a weight of 3, as well as similarly
peaked functions for filaments and clusters. However, as we expect MG to be
the strongest in low-density regions we restrict ourselves to upweighing void
or wall galaxies only.

Yet another idea of using the environmental classification of galaxies as a
mark would be to further increase the anti-correlation present in low density
regions. This can be accommodated with the following mark

$m(\mathbf{x})=\begin{cases}-1\quad\text{if void}\\\
1\quad\text{else}\end{cases}\quad\text{Void${}_{\rm AC}$},$ (46)

where we abbreviated _anti-correlation_ with AC. In principle, there is no
difference if we switch signs of this mark because from Eq. (28) it is clear
that any overall factor of the marks would be cancelled by the division of the
normalisation. This mark leaves galaxy pairs that are in voids unweighted as
well as galaxy pairs not in voids. However, if one galaxy is in a void and the
other is not, the weight will be -1 thereby creating an anti-correlation.

Marks based on the tidal tensor components may appear promising to go beyond
the local density. An interesting quantity first introduced by Heavens &
Peacock (1988) and then used by Alam et al. (2019), is the tidal torque,
defined as

$t(\mathbf{x})=\frac{1}{2}\\{(\lambda_{3}-\lambda_{2})^{2}+(\lambda_{3}-\lambda_{1})^{2}+(\lambda_{2}-\lambda_{1})^{2}\\},$
(47)

with $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are the eigenvalues of the
tidal tensor. The larger the difference between the eigenvalues the more
anisotropic is the structure. Hence we expect the tidal torque to be large for
filaments and walls and small for clusters or voids. Another field depending
directly on tidal tensor components is the tidal field, also known as the
second Galileon $\mathcal{G}_{2}$ (Nicolis et al., 2009), and which was used
extensively for the emergence of non-local bias between galaxies and dark
matter (Chan et al., 2012). The tidal field is defined as
$\mathcal{G}_{2}=(\partial_{i}\partial_{j}\Phi)^{2}-(\partial^{i}\partial_{i}\Phi)^{2}$,
where we can identify the components of the tidal tensor as introduced in Eq.
(44). In practise we investigate two separate marks consisting of the tidal
field and tidal torque as they are, respectively. Using these fields in this
way should give an insight on the suitability of the tidal field or tidal
torque to disentangle MG from GR.

In Figure 3 we present the different marked correlation functions introduced
in this section. For the cluster marked correlation function we see a strong
signal on very small scales which relates to the correlation between galaxies
insides clusters. The compensation feature on scales between 20$\,h^{-1}\,{\rm
Mpc}$ and 60$\,h^{-1}\,{\rm Mpc}$ comes from less clustered regions around
clusters and is similar, although reversed, to the compensation seen in the
void-galaxy cross-correlation function (e.g. Aubert et al., 2022; Hamaus et
al., 2022). The filament and wall marked correlation functions show
progressively less signal as the clustering of galaxies inside walls and
filaments are closer to the total clustering of all galaxies. Notably, if
voids are considered, the observed signal below unity implies that void
galaxies are less clustered compared to the total clustering. The large signal
of the cluster marked correlation function comes at the cost of larger errors
due to small amount of galaxies residing in clusters. In general we aim for
mark correlation functions with a signal different from unity over a wide
range of scales as this might lead to differences at those scales between MG
and GR. On the other side, if the marked correlation function stays very close
to unity on most scales, then any possible difference between MG and GR can
only originate from the clustering itself. The latter is matched between MG
and GR in the ELEPHANT simulations to fit observations as described in Section
4. For these reasons, the $\textrm{Void}_{\textrm{AC}}$ mark and particularly
the $\textrm{Wall}_{\textrm{AC}}$ mark are of strong interest as they exhibit
a signal up to large scales. Looking at the lower panel of Fig. 3, we can see
a strong signal when the tidal field is used, which extends to large scales.
The same, although with considerably less amplitude on small scales, is found
for the tidal torque. Hence, these marks are also interesting candidates to be
investigated to discriminate between MG and GR. Although an impact of the mark
is a necessary prerequisite, it is not sufficient to guarantee a
disentanglement of GR from MG because the signal could be the same in MG and
GR, nevertheless.

It has to be noted that the marked correlation functions shown in Figure 3 are
not corrected for a possible bias due to the estimation of the mark on a
discrete catalogue, as we discuss in the next section. Hence, the exact
amplitude of the measurements might be subject to changes if such a correction
is applied. For the marks based on the environmental classification we do not
expect this bias to be particularly strong because possible miss-
classifications, originating from a biased estimate of the density field,
should not affect every galaxy in a catalogue.

Figure 3: Summary of the different marked correlation functions using marks
based on the tidal field, tidal torque (both in the lower panel) or large-
scale environment (upper panel). The black dashed line indicates an amplitude
of 1. Curves represent the mean taken over 5 realisations and the error
corresponds to the mean standard deviation over 5 realisations. The
measurements have not been corrected for any form of bias due to shot-noise
effects.

###  5.2 Anti-correlating galaxies using local density

Until now, we have seen that marks based on the local density are particularly
simple and introducing an anti-correlation with the mark appears to be
promising regarding the discrimination between GR and MG. Therefore, we
propose the following mark function based on the hyperbolic tangent,
satisfying both aforementioned advantages,

$m(\mathbf{x})=M_{\tanh}[\delta_{R}(\mathbf{x})]\equiv\tanh(a\,(\delta_{R}(\mathbf{x})+b)),$
(48)

where $a$ and $b$ are parameters controlling how steeply the transition from
-1 to 1 takes place and where the transition happens, respectively. In
general, we could use a third parameter $c$ as an overall factor in front of
the hyperbolic tangent but constant factors can be pulled out of the mean and
hence are cancelled by the normalisation in Eq. (28). It is worth mentioning
the fact that theoretical modelling of marks based on the environmental
classification, such as the $\textrm{Void}_{\textrm{AC}}$ mark, might be
particularly challenging as it is not straightforward to express the mark in
terms of the density contrast. From a theoretical perspective, marked
correlation functions with marks based on the density are more tractable.
Furthermore, discreteness effects, arising in the density estimation itself
can be more easily corrected for in the measurement of the marked correlation
function as elucidated in the next section.

##  6 Propagation of discreteness effects of the mark estimation into weighted
correlation functions

When dealing with finite point sets we can assume the sampling process to be
locally of Poisson nature (Layzer, 1956). This means that the number of points
found in some small-enough grid cells appears as being drawn from a Poisson
distribution with some expectation value. However, the expectation value of
the local Poisson process does have a PDF on its own. The PDF from which the
expectation values is drawn is continuous and describes the density field
globally. If this PDF is a Dirac delta function then the expectation value of
the Poisson process is the same everywhere and the moments estimated from the
sample points coincide with the moments from the continuous PDF. However, if
the continuous PDF is not a Dirac delta function, as is the case for the
cosmological density field, then the estimated moments contain a bias with
respect to the true moments of the continuous PDF. This bias is usually called
shot noise or Poisson noise in the literature. In the power spectrum
estimation, the shot noise appears as an additive constant for all scales in
$k$-space. In the 2PCF instead, the shot noise emerges only at zero lag, that
is at a pair separation of zero. Hence, shot noise is inherently a problem of
correlating a point with itself. When we use the density field inside a mark
function, we use a smoothed version of the true field. We spread points over a
finite volume leading to self-correlations also at non-zero pair separation
and in turn to shot-noise effects. Intuitively, this can be understood in the
following manner: in the unsmoothed case all points are infinitely small dots,
while in the smoothed case the points are represented by circles with a non-
zero radius. Inside this radius one point can be correlated with itself.

In order to precisely understand how shot noise affects marked correlation
functions we have to do a small detour and carefully distinguish between the
statistical properties of the true density contrast $\delta(\mathbf{x})$, the
smoothed true density contrast $\delta_{R}(\mathbf{x})$, and the respective
quantities estimated from finite point sets, hereby denoted with an $f$ in the
subscript $\delta_{f}(\mathbf{x})$ and $\delta_{Rf}(\mathbf{x})$. The weighted
correlation function estimated from a finite point set can be written as

$1+W_{f}(\mathbf{r})=\frac{w_{f}(\mathbf{r})}{\bar{m}_{f}^{2}},$ (49)

where we defined the quantity $w_{f}(\mathbf{r})$ as

$w_{f}(\mathbf{r})\equiv\langle
M[\delta_{Rf}(\mathbf{x})](1+\delta_{f}(\mathbf{x}))M[\delta_{Rf}(\mathbf{x}+\mathbf{r})](1+\delta_{f}(\mathbf{x}+\mathbf{r}))\rangle,$
(50)

and $\bar{m}_{f}=\frac{1}{V}\int_{V}m_{f}({\bf x})\rho_{f}({\bf
x})/\bar{\rho}_{f}{\mathrm{d}}^{3}x$ is the mean mark taken over the points,
i.e. weighted by the density. In Eq. (49) both $w_{f}(\mathbf{r})$ and
$\bar{m}_{f}$ are expected to be sensitive to the noise induced by the auto-
correlation of objects with themselves, and we will denote the corresponding
shot-noise free signals $w$ and $\bar{m}$. Indeed, Eq. (50) shows that there
is a mark function $M$ of the smoothed density field that is multiplied by the
density field itself. This constitutes the main source of shot noise that is
expected to happen even at large separation $r$, where there is no overlap
between the smoothing kernels. In this section, we first show the effect of
shot noise on the marked correlation function for a specific mark and then
devise a general method to correct for shot noise.

###  6.1 A toy model

In order to understand how the shot noise propagates into $w_{f}(\mathbf{r})$
and $\bar{m}_{f}$, we focus on a very simple mark function defined by
$M[\delta_{Rf}({\bf x})]=\delta_{Rf}({\bf x})$ and $M[\delta_{Rf}({\bf x}+{\bf
r})]=1$. This corresponds to a marked correlation function in which only one
point of the pair is weighted by the density contrast and the other point
stays unweighted. In this instructive case we can split $w_{f}(\mathbf{r})$
into three terms

$w_{f}(\mathbf{r})=\sigma_{Rf}^{2}+\xi_{Rf}({\bf r})+\zeta_{Rf}({\bf r}),$
(51)

where the individual contributions are given by

$\begin{split}\sigma_{Rf}^{2}&=\langle\delta_{Rf}(\mathbf{x})\delta_{f}(\mathbf{x})\rangle\\\
\xi_{Rf}({\bf
r})&=\langle\delta_{Rf}(\mathbf{x})\delta_{f}(\mathbf{x}+\mathbf{r})\rangle\\\
\zeta_{Rf}({\bf
r})&=\langle\delta_{Rf}(\mathbf{x})\delta_{f}(\mathbf{x})\delta_{f}(\mathbf{x}+\mathbf{r})\rangle,\end{split}$
(52)

and consist of correlators between the smoothed and unsmoothed density field
estimated on a finite point set. Given that the smoothed density field
$\delta_{Rf}(\mathbf{x})$ is related to the density field
$\delta_{f}(\mathbf{x})$ through the convolution

$\delta_{Rf}(\mathbf{x})=\frac{1}{a^{3}}\int_{\mathbf{x}^{\prime}}F\left(\frac{\mathbf{x}-\mathbf{x}^{\prime}}{a}\right)\delta_{f}(\mathbf{x}^{\prime})\text{d}^{3}x^{\prime},$
(53)

one can immediately see that the three contributions in Eq. (51) are involving
integrals over two- and three-point correlation functions of the density
field. In general, $n$-point correlation functions are affected by shot noise
as (see Chan & Blot, 2017)

$\begin{split}\Xi_{f}(\mathbf{r}_{1},..,\mathbf{r}_{n-1})&=\langle\delta_{f}(\mathbf{x})\delta_{f}(\mathbf{x}+\mathbf{r}_{1})...\delta_{f}(\mathbf{x}+\mathbf{r}_{n-1})\rangle_{c}\\\
&=\Xi(\mathbf{r}_{1},...,\mathbf{r}_{n-1})+\sum_{m=1}^{n-1}\frac{1}{\bar{n}^{m}}\mathcal{A}^{(n)}_{m}(\mathbf{r}_{1},...,\mathbf{r}_{n-1})\end{split}$
(54)

where the function $\mathcal{A}^{(n)}_{m}$ contains all the scale dependency
of the shot-noise contribution to the $n$-point correlation function at the
respective order in $\bar{n}$, the mean density of points in the volume $V$.
Therefore, the shot noise takes the form of a power series in $1/\bar{n}$. In
particular for the 2PCF we have

$\xi_{f}(\mathbf{r})=\xi(\mathbf{r})+\frac{\delta_{D}(\mathbf{r})}{\bar{n}},$
(55)

and for the three-point correlation function

$\begin{split}\zeta_{f}(\mathbf{r},\mathbf{s})&=\langle\delta_{f}(\mathbf{x})\delta_{f}(\mathbf{x}+\mathbf{r})\delta_{f}(\mathbf{x}+\mathbf{s})\rangle_{c}\\\
&=\zeta(\mathbf{r},\mathbf{s})+\frac{1}{\bar{n}}\left[\delta_{D}(\mathbf{r})\xi(\mathbf{s})\right.\\\
&\left.\quad+\delta_{D}(\mathbf{r}-\mathbf{s})\xi(\mathbf{r})+\delta_{D}(\mathbf{s})\xi(\mathbf{r}-\mathbf{s})\right]\\\
&\quad+\frac{1}{\bar{n}^{2}}\delta_{D}(\mathbf{r})\delta_{D}(\mathbf{s}).\end{split}$
(56)

As a result, one can express each individual term of Eq. (51) in terms of the
true signal and a shot-noise contribution (depending on the number density of
objects) as

$\begin{split}\sigma_{Rf}^{2}&=\sigma_{R}^{2}+\frac{F({\bf 0})}{\bar{N}}\\\
\xi_{Rf}({\bf r})&=\xi_{R}({\bf r})+\frac{F({\bf r}/a)}{\bar{N}}\\\
\zeta_{Rf}({\bf r})&=\zeta_{R}({\bf r})+\frac{\xi({\bf
r})}{\bar{N}}\left[F({\bf 0})+F({\bf r}/a)\right],\end{split}$ (57)

where we introduce $\bar{N}=a^{3}\bar{n}$ that corresponds to the mean number
of objects per grid cell. These noise contributions are obtained by using Eq.
(53) in Eq. (52), inserting Eq. (55) and (56), and integrating out the Dirac
delta functions where applicable. It has to be noted that we only report noise
contributions in Eq. (57) that are not proportional to Dirac delta functions,
as these would appear at zero lag only and hence be irrelevant for our
considerations. By utilising the aforementioned splitting into signal and
noise we can write Eq. (51) as

$w_{f}(\mathbf{r})=w(\mathbf{r})+{\epsilon}_{w}(\mathbf{r}),$ (58)

where $w$ is the true signal and the shot-noise contribution is formally
expressed as

${\epsilon}_{w}(\mathbf{r})=\frac{1}{\bar{N}}\left(1+\xi({\bf
r})\right)\left[F({\bf 0})+F({\bf r}/a)\right].$ (59)

Equation (59) shows that even if on a scale $r$ larger than the smoothing
scale (when $F({\bf r}/a)=0$) there is still a large-scale contribution to the
shot noise due to $F({\bf 0})$. In addition, the large-scale contribution is
expected to decrease when increasing the order of the MAS ($F({\bf 0})$ is
decreasing). That is the reason why, in general, increasing the order of the
MAS is reducing the intrinsic shot-noise contribution to the signal.

Following the same reasoning, it is straightforward to show that with the toy
model the shot-noise affected mean mark $\bar{m}_{f}$, estimated from a
discrete set of objects, can be related to the true mean mark $\bar{m}$ via

$\bar{m}_{f}=\bar{m}+{\epsilon}_{\bar{m}},$ (60)

where

${\epsilon}_{\bar{m}}=\frac{F({\bf 0})}{\bar{N}}.$ (61)

Finally, by combining Eq. (58) and (60) we can show that the shot-noise-
corrected weighted correlation function $1+W(\mathbf{r})$ can be expressed as

$1+W({\bf
r})=\frac{w_{f}(\mathbf{r})-{\epsilon}_{w}(\mathbf{r})}{\bar{m}_{f}-{\epsilon}_{\bar{m}}}.$
(62)

This demonstrates that the shot-noise correction on the marked correlation
function implies to correct both the numerator $w_{f}(\mathbf{r})$ and
denominator $\bar{m}_{f}$ in order to properly extract the true signal. This
is at odd with usual shot-noise corrections on $n$-point correlation functions
that are only additive.

There is another subtlety due to the fact that we assign the mark field back
on the galaxies to measure the weighted correlation function. This back-
assignment is done with a specific scheme in the sense that we check in which
grid cell a galaxy is located and assign the mark corresponding to that grid
cell, thereby introducing another smoothing of the field with a NGP kernel. In
our computation this leads to an additional convolution for the field
$\delta_{Rf}(\mathbf{x})$. We show in Section B that this additional
convolution is equivalent to a single one with a kernel that is a convolution
of both a PCS and a NGP kernel, that is, a quartic kernel. Hence, in the
actual calculation of the analytic shot noise as in Eq. (59) and (61), a
quartic kernel has to be used, which explicit expression can be found in the
appendix of Chaniotis & Poulikakos (2004).

In order to validate the analytic prediction of the noise in
$w_{f}(\mathbf{r})$ and $\bar{m}_{f}$ we use five realisations of Covmos as
introduced in Section 4. The goal is to have realisations at different number
densities to assess the behaviour of shot noise as a function of $\bar{N}$.
Therefore, for each realisation we deplete the catalogue by randomly throwing
away points down to the desired density. The exact densities are motivated by
applying the shot-noise correction later to the ELEPHANT simulation suite,
which has much lower point densities compared to Covmos. Hence by depleting
the Covmos realisations down to
$\\{1.7\%,1.53\%,1.36\%,1.19\%,1.02\%,0.85\%,0.68\%,0.51\%\\}$ we generate
catalogues with the same $\bar{N}$ as in the ELEPHANT suite if they were
depleted down to $\\{100\%,90\%,80\%,70\%,60\%,50\%,40\%,30\%\\}$ with 64 grid
cells per dimension. The depletion is done to match the density of points in
grid cells and not in the full volume because the shot-noise behaviour is a
power series in $1/\bar{N}$. The depletion is repeated 100 times followed by a
mean to minimise the sample variance coming from the stochasticity of the
random depletion process. We need to carefully distinguish the five
independent realisations of Covmos from the depletion realisations used to get
a converged result for a depleted catalogue, which has to be done for each of
the five independent realisations.

As we have seen in Eq. (62) we need to measure
$w_{f}(\mathbf{r})=(1+W_{f}(\mathbf{r}))\bar{m}_{f}$ as well as $\bar{m}_{f}$
and those can be straightforwardly computed from the weighted correlation
function at each level of depletion. In the upper panel of Figure 4 we present
the measurements in one Covmos realisation of $w_{f}(\mathbf{r})$ (blue
points) and $\bar{m}_{f}$ (orange points) as a function of $1/\bar{N}$. The
scale at which we plot $w_{f}(\mathbf{r})$ is fixed to a bin of
$20\leavevmode\nobreak\ \,h^{-1}\,{\rm Mpc}$. We can already see that we there
is a linear relation with $1/\bar{N}$ as predicted by the expression in Eq.
(59) and Eq. (61). The solid and dashed curve refer to the analytical
prediction using the depletion case to 1.7% (the second data point from the
left) as an anchor. That anchor is needed to obtain a noiseless signal by
correcting for shot noise and then add to the true signal the noise
contribution, as a function of $1/\bar{N}$, to obtain the curve. By doing so,
the relative difference between the prediction and the measurement is exactly
zero by construction for this depletion as can be seen in the lower panel of
Figure 4. Moreover, even for the other data points at different levels of
depletion we can predict the expected signal with high accuracy. The relative
difference is at the sub-percent level for both $w_{f}(\mathbf{r})$ and
$\bar{m}_{f}$.

Figure 4: Analytical correction for $w_{f}(\mathbf{r})$ and $\bar{m}_{f}$ in
case of the toy model where only one galaxy in the pair is weighted by the
density contrast $\delta_{Rf}$. The points refer to the measured data from one
of the Covmos realisations, in blue for $w_{f}$ and in orange for
$\bar{m}_{f}$. The upper panel shows the measurements alongside the analytical
prediction and the bottom panel presents the relative difference between the
measurements and the theory. We show only one Covmos realisations for which we
used 100 realisations of depletions to obtain the depleted catalogues and the
scale bin in $r$ for $w_{f}$ is fixed to be close to 20$\,h^{-1}\,{\rm Mpc}$.
For the analytic correction we use the depletion down to 1.7% as an anchor and
computed from there the expected signal using Eq. (59) and (61).

Now that we have established the correctness of our analytical predictions for
$w_{f}(\mathbf{r})$ and $\bar{m}_{f}$ individually, we can check how well they
perform when combining them into $1+W_{f}(\mathbf{r})$, as shown in Figure 5.
In the upper panel we present the mean over the five Covmos realisations at
different levels of depletion as indicated with different colours in the
legend. As expected, since the kernel and 2PCF drop off at large separations,
differences in the curves are only evident on smalle scales. This is further
underlined by the lower panel where the relative difference between the
depletion down to 1.7% and the undepleted case is shown in black. This curve
refers to the difference between the two if we would not have applied any
correction and only data with a depleted number density of 1.7% would be
available. For $1+W_{f}(\mathbf{r})$ the relative difference can reach more
than 10% on small scales but at scales above around 60$\,h^{-1}\,{\rm Mpc}$
the depleted and undepleted case lay within 1% and the effect of shot noise
becomes negligible. This is somewhat expected due to the smoothing of the
density field, as the quartic kernel decreases down to zero over the course of
2.5 grid cells, which in Covmos corresponds to $\approx$40$\,h^{-1}\,{\rm
Mpc}$. It is important to note here that, although we expect shot noise to be
stronger when correlating within the volumes of the smoothing kernels, it is
peculiar to the toy model that the shot-noise contribution does only contain
linear factors of the kernel with and without the 2PCF. It can be shown in the
more general case that if one weights both galaxies in a pair by the
associated density field, then the shot noise will contain contributions from
a convolution of two quartic kernels resulting in a nonic kernel, which is
much more extended in configuration space. In contrast to the black curve that
has no correction, we show the relative difference of the analytical
correction for $1+W_{f}(\mathbf{r})$ to the undepleted case in red. To have a
fair comparison we corrected the 1.7% case down to the density of the
undepleted realisations, as these still have a finite, yet very high density.
The analytical correction reproduces the undepleted measurements to within 1%
relative difference on all scales. We conclude that for the toy model we are
able to analytically predict the shot noise. Moreover, we show that even with
this simple toy model the shot noise acquires a non-trivial scale dependency.
In the next section we extend this formalism to general weighted correlation
functions and describe a procedure to estimate the signal without having to
rely on an analytical model.

Figure 5: Results of the analytic shot-noise correction applied to the toy
model for which only one point in a pair is weighted by the density contrast
$\delta_{Rf}$. The upper panel shows the full weighted correlation function
$1+W_{f}$ as a mean over the 5 Covmos realisations, where different colours
denote different levels of depletion. Errorbars are computed by taking the
mean standard deviation over 5 realisations. To obtain the depleted catalogues
we took the mean over 100 depletions. The lower panel shows the relative
difference of the analytically corrected result to the undepleted case as a
red dashed line. The solid black line refers to the relative difference of the
depletion level $1.7\%$ to the undepleted case, which illustrates the effect
of no correction. Horizontal dashed lines in black indicate levels of relative
differences of $\pm 1\%$ and the vertical dashed line in grey refers to side
length of one grid cell. We used 64 grid cells per dimension and a PCS MAS to
obtain the density field on the grid.

###  6.2 A general model

Building on top of the results obtained with the toy model, we can devise a
general model that has a mark function expandable in powers of the density
contrast as

$M[\delta_{Rf}(\mathbf{x})]=\sum_{i=0}^{\infty}\frac{c_{i}}{i!}\delta_{Rf}^{i}(\mathbf{x}),$
(63)

where $c_{i}$ are the coefficients of the Taylor series. Plugging the series
expansion in Eq. (63) into Eq. (49) we arrive at

$w_{f}({\bf
r})=\sum_{i,j}\frac{c_{i}c_{j}}{i!j!}\langle\delta_{Rf}^{i}(\mathbf{x})(1+\delta_{f}(\mathbf{x}))\delta_{Rf}^{j}(\mathbf{x}+\mathbf{r})(1+\delta_{f}(\mathbf{x}+\mathbf{r}))\rangle$
(64)

and

$\bar{m}_{f}=\sum_{i}\frac{c_{i}}{i!}\langle\delta_{Rf}^{i}(\mathbf{x})\frac{\rho_{f}(\mathbf{x})}{\bar{\rho}_{f}}\rangle.$
(65)

It is evident from these expressions that by weighting both galaxies in a
given pair, the resulting marked correlation function will contain auto-
correlation contributions of the mark with itself. If for example the weight
is constructed by an external catalogue of voids then the weighted correlation
function will consist of a smoothed version of the void auto-correlation and
void-galaxy cross-correlation functions. In Section A, we give some further
insights in how the weighted correlation function can be split up into two
auto- and one cross-correlation function for certain weighting schemes.

In order to work out the shot-noise contribution to $w_{f}(\mathbf{r})$ and
$\bar{m}_{f}$ we can use, analogously to Eq. (54), the relation

$\begin{split}&\langle\delta_{Rf}^{i}(\mathbf{x})(1+\delta_{f}(\mathbf{x}))\delta_{Rf}^{j}(\mathbf{x}+\mathbf{r})(1+\delta_{f}(\mathbf{x}+\mathbf{r}))\rangle&\\\
&=\langle\delta_{R}^{i}(\mathbf{x})(1+\delta(\mathbf{x}))\delta_{R}^{j}(\mathbf{x}+\mathbf{r})(1+\delta(\mathbf{x}+\mathbf{r}))\rangle+\sum_{p=1}^{i+j+1}\frac{1}{\bar{N}^{p}}\mathcal{B}^{(i,j)}_{p}(\mathbf{r}),\end{split}$
(66)

where $\mathcal{B}^{(i,j)}_{p}(\mathbf{r})$ contains the shot-noise
contribution for a given $(i,j)$ proportional to the inverse of $\bar{N}$ to
the power of $p$. Inserting this expression into Eq. (64) we obtain

$\begin{split}w_{f}(\mathbf{r})=&\sum_{i,j}\frac{c_{i}c_{j}}{i!j!}\langle\delta_{R}^{i}(\mathbf{x})(1+\delta(\mathbf{x}))\delta_{R}^{j}(\mathbf{x}+\mathbf{r})(1+\delta(\mathbf{x}+\mathbf{r}))\rangle\\\
&+\sum_{i,j}\frac{c_{i}c_{j}}{i!j!}\sum_{p=1}^{i+j+1}\frac{1}{\bar{N}^{p}}\mathcal{B}^{(i,j)}_{p}(\mathbf{r})\\\
=&w(\mathbf{r})+{\epsilon}_{w}(\mathbf{r}),\end{split}$ (67)

where we identified in the second equality the first sum to be the desired
true signal $w(\mathbf{r})$ and the second sum to be the shot-noise
contribution ${\epsilon}_{w}(\mathbf{r})$. Similarly, for $\bar{m}_{f}$ we
obtain

$\bar{m}_{f}=\bar{m}+\sum_{i}\frac{c_{i}}{i!}\sum_{p=1}^{i}\frac{1}{\bar{N}^{p}}\mathcal{B}^{(i)}_{p}=\bar{m}+{\epsilon}_{\bar{m}}.$
(68)

At this point it is clear that the double sum can be written as a power series
of $1/\bar{N}$ such that

${\epsilon}_{w}(\mathbf{r})=\sum_{p=1}^{\infty}\frac{1}{\bar{N}^{p}}{\epsilon}_{w,p}(\mathbf{r})\qquad{\epsilon}_{\bar{m}}=\sum_{p=1}^{\infty}\frac{1}{\bar{N}^{p}}{\epsilon}_{\bar{m},p},$
(69)

with correspondingly defined ${\epsilon}_{w,p}(\mathbf{r})$ and
${\epsilon}_{\bar{m},p}$. The correction in the general case is therefore,
analogously to Eq. (62),

$1+W(\mathbf{r})=\frac{w_{f}(\mathbf{r})-{\epsilon}_{w}(\mathbf{r})}{(\bar{m}_{f}-{\epsilon}_{\bar{m}})^{2}}.$
(70)

The power series in $1/\bar{N}$ (which in principle extends to infinite order)
together with the fact that, following Eq. (55) and (56), shot noise of
$n$-point correlation functions is scale-dependent and contains $(n-1)$-point
correlation functions, makes an analytic correction for the general case
untractable. It would require in particular to compute higher-order
correlation functions, which are computationally expensive. One might think
that a simple truncation of the Taylor expansion would solve the problem, but
to avoid computing four-point correlators and above, the Taylor expansion
would need to be cut already at linear order. Moreover, the conversion of
moments into cumulants might lead to significant contributions from higher-
order correlators at low order in $1/\bar{N}$. In the following we outline an
approach to circumvent analytical computation and that uses the resummation of
contributions into a power series in $1/\bar{N}$.

The quantities $w_{f}=(1+W_{f}(\mathbf{r}))\bar{m}_{f}^{2}$ as well as
$\bar{m}_{f}$ are directly measurable from simulations for a given mark. We
propose therefore an algorithm consisting of a polynomial fit through
measurements of $w_{f}(\bar{N})$ and $\bar{m}_{f}(\bar{N})$ made at different
levels of depletion, that is, at different values of $1/\bar{N}$. For
$w_{f}(\mathbf{r})$, the fit is done with the same polynomial order for each
bin in $r$ but with separate coefficients, which is necessary since the shot
noise ${\epsilon}_{w}(\mathbf{r})$ is scale-dependent. With such a polynomial
we can simply read off the noiseless signal from the $y$-axis intersection as
this gives the extrapolation to $1/\bar{N}=0$, i.e. infinite densities. It is
important to note that truncating the fit at some polynomial order is not the
same as truncating the Taylor expansion in $\delta_{Rf}$ as the linear
coefficient in the power series contains the resummed contributions from all
higher-order correlators as well. To test this approach and find the best
order of polynomial to fit, we used the same depletion levels as described in
the previous section.

In Figure 6 we present the results from polynomial fits to the quantities
$w_{f}(\mathbf{r})$ and $\bar{m}_{f}$ as appearing in Eq. (70). It is evident
that for the general case, higher-order shot-noise contributions play an
important role resulting in a more curved shape due to quadratic and cubic
dependencies on $1/\bar{N}$. Therefore, a simple linear fit is not sufficient
anymore and at least 2nd- or 3rd-order polynomials are to be used. Going to
even higher orders, as we show with a 4th-order polynomial in purple, the
behaviour outside of the fitted range becomes more unstable and can lead to
severe over- or under-estimation of the true signal at the $y$-axis
intersection. Moreover, since we only employ 8 data points, it is crucial to
keep the polynomial order as low as possible since otherwise an overfitting of
the data might happen. In contrast to the behaviour for $w_{f}(\mathbf{r})$,
the dependency of $\bar{m}_{f}$ on $1/\bar{N}$ appears to be much more linear
and fits with 1st- or 2nd-order polynomials should be sufficient to recover
the true signal accurately.

Figure 6: Fitting procedure to obtain the shot-noise-corrected signal in case
of $w_{f}(\mathbf{r})$ (left panel) and $\bar{m}_{f}$ (right panel) for the
$\tanh$-mark in the configuration $(a,b)=(0.6,-0.5)$. Following Eq. (69) we
present the fits in dependency of $1/\bar{N}$, the reciprocal of the average
number of points per grid cell. We show only one realisation of Covmos and the
scale bin in $r$ for $w_{f}$ is fixed to be close to $20$$\,h^{-1}\,{\rm
Mpc}$. The depleted measurements were obtained by taking the mean over 100
realisations of depletions. The orange points refer to the fitted measurements
and the blue point is the undepleted reference (not included in the fit).
Differently coloured lines refer to different orders in the polynomials we
used to fit the data. Errorbars of the orange points are obtained by taking
the mean standard deviation over the 100 realisations. Since the depletion
down to 1.7% (first orange point from the left) mimics the undepleted ELEPHANT
density we use for this point an error that is 10% of the minimum uncertainty
over the remaining depletions.

In Figure 7 the performance of the different correction orders as well as the
weighted correlation function $1+W_{f}(\mathbf{r})$ and the quantity
$w_{f}(\mathbf{r})$ are shown. The diverging behaviour of
$1+W_{f}(\mathbf{r})$ in the upper left panel at the depletions down to
$0.68\%$ and $0.85\%$ (olive green and brown lines, respectively) can be
understood when looking at the corresponding points of the mean mark
$\bar{m}_{f}$ in Figure 6, that is the second and third to last data point at
$1/\bar{N}$ of around 2.0 and 1.6. The mean mark is very close to zero in that
case and therefore $1+W_{f}(\mathbf{r})$ acquires a very large amplitude due
to the division by $\bar{m}_{f}^{2}$. This behaviour is somewhat peculiar to
marks that can switch signs, as the $\tanh$-mark, because in certain cases
this can lead to a mean mark which is very close to zero. Moreover, this can
result in a turn-around in the dependency on $1/\bar{N}$ as seen for the last
two depletions, $0.68\%$ and $0.51\%$, both in Figure 6 and 7. While a very
small mean mark does not appear to be problematic for the fit, it can be an
issue if the true mean mark is very close to zero. Since we use only 8 data
points in the fit, we have a limited accuracy on the recovery of
$\bar{m}_{f}$. This can become problematic as soon as the amplitude of the
recovered mean mark approaches the accuracy of the fit, leading to very large
relative uncertainties on $\bar{m}_{f}$ and on $1+W_{f}(\mathbf{r})$. In this
case, the accuracy of the polynomial fit is not enough to properly recover
small mean marks and the results should not be trusted. Even though one could
try to mitigate this issue and improve the accuracy of the fit with better
estimates of the data points from larger sets of depleted catalogues, in
general, we advise against using marks with a recovered mean mark being very
close to zero.

In the upper panel to the right of Figure 7, we show the measurements of
$w_{f}(\mathbf{r})$. While the noise behaviour on small scales appears to be
less severe compared to $1+W_{f}(\mathbf{r})$, on large scales we can observe
a constant offset. Focusing on the relative differences in the lower left
panel of Figure 7, it is evident that for this particular mark the effect of
shot noise diminishes at scales of around 60$\,h^{-1}\,{\rm Mpc}$ as the
relative difference between the undepleted case and uncorrected measurement at
1.7% depletion shrinks to below 5%. This is in contrast to the toy model in
Figure 5 as here the 5% border is crossed already at scales of around
40$\,h^{-1}\,{\rm Mpc}$. This illustrates the fact that the shot-noise
behaviour can have different amplitudes and scale-dependency that is subject
to the chosen mark. The coloured lines in the lower panels refer to different
orders in the polynomial fit used for $w_{f}(\mathbf{r})$ and $\bar{m}_{f}$.
From this we can conclude that fitting the behaviour of $w_{f}(\mathbf{r})$
with a 3rd-order and $\bar{m}_{f}$ with a 2nd-order polynomial results in a
satisfactory performance and should be used hereafter as the adequate shot-
noise correction. With this choice, the relative difference in
$1+W_{f}(\mathbf{r})$, depicted by the orange line in Fig. 7, is within 5%
across all scales all the way up to 150$\,h^{-1}\,{\rm Mpc}$.

Figure 7: Results of the fitted shot-noise correction for the $\tanh$-mark
with $(a,b)=(0.6,-0.5)$. The left panels present $1+W_{f}$ while the right
panel shows $w_{f}=(1+W_{f})\bar{m}_{f}^{2}$. Different colours in the upper
panel refer to the mean over 5 Covmos realisations and errorbars correspond to
the mean standard deviation over those 5 realisations. Depleted measurements
were obtained by taking the mean over 100 realisations of depletions. The
bottom panels show the relative differences of the corrected signal from the
fit to the undepleted reference case. Different colours in the lower panels
refer to different orders used in the polynomial fit and e.g. ’3rd/2nd’
indicates that a 3rd- and 2nd-order polynomial fit was used for $w_{f}$ and
$\bar{m}_{f}$, respectively. The horizontal dashed lines in black corresponds
to a relative difference of $\pm 5\%$ and the vertical dashed line in grey
refers to the side length of one grid cell.

Now that we have an optimal choice for the polynomial orders to describe the
shot noise as a function of $1/\bar{N}$, we need to assess how many
realisations of depletions yield converged results for the polynomial fits.
Since the process of depletion consists in randomly throwing away a number of
points such that we end up with some desired percentage of the original
points, it is inherently noisy and should be repeated several times. The aim
is to obtain a representation of the original catalogue with a lower density
that looks like as if the simulation has been run with lesser points in the
first place. In Figure 8 we show the best correction as obtained from Figure
7, being the 3rd/2rd-order polynomial, and compute relative differences for
different amounts of realisations of depletion. As we can see, using 30
realisations or more, the curves do not differ substantially and results can
be considered converged. Even with only 10 or 20 realisations at hand the
performance is nevertheless acceptable and well within 5% except for the
lowest bin in $r$. As a conservative choice, we will use 30 realisations in
the following. This should allow mitigating sample variance and not affecting
resulting corrections.

Figure 8: Test of convergence of the shot-noise correction for different
amounts of depletion realisations. Black dashed lines refer to relative
differences of $\pm 5\%$ and colours denote the number of realisations for the
depletions. Relative differences are shown for the total $1+W_{f}(r)$ where
both $w_{f}(r)$ and $\bar{m}_{f}$ are corrected for the $\tanh$-mark with
$(a,b)=(0.6,-0.5)$.

###  6.3 Shot noise in redshift space

In real observations, the quantity of interest in clustering analyses are
usually the multipoles of the anisotropic 2PCF. In the following we show that
the correction to the marked correlation function in redshift space is similar
to that in real space. We indicate redshift-space quantities via their
explicit dependency on the separation $s$ and angle $\mu$, but also the mean
mark has to be understood as measured in redshift space. For the full
anisotropic marked correlation function $\mathcal{M}_{f}(s,\mu)$, the shot
noise enters in the following way

$\begin{split}\mathcal{M}_{f}(s,\mu)&=\frac{1+W_{f}(s,\mu)}{1+\xi(s,\mu)}=\frac{(1+W(s,\mu))\bar{m}^{2}+{\epsilon}_{w}(s,\mu)}{(1+\xi(s,\mu))(\bar{m}+{\epsilon}_{\bar{m}})^{2}}\\\
&=\mathcal{M}(s,\mu)\frac{\bar{m}^{2}}{(\bar{m}+{\epsilon}_{\bar{m}})^{2}}+\frac{{\epsilon}_{w}(s,\mu)}{(1+\xi(s,\mu))(\bar{m}+{\epsilon}_{\bar{m}})^{2}},\end{split}$
(71)

where we used the fact that the unweighted 2PCF is not affected by shot noise
and hence does not have the subscript $f$. Since the shot noise does contain
$n$-point correlators, it will acquire an angle dependency as well. After
decomposing the anisotropic marked correlation function into multipoles
$\mathcal{M}_{f,\ell}(s)$ we obtain

$\begin{split}\mathcal{M}_{f,\ell}(s)=\mathcal{M}_{\ell}(s)\frac{\bar{m}^{2}}{\bar{m}_{f}^{2}}+\frac{2\ell+1}{2\bar{m}_{f}^{2}}\int_{-1}^{1}P_{\ell}(\mu)\widetilde{{\epsilon}}_{w}(s,\mu)\text{d}\mu,\end{split}$
(72)

with

$\begin{split}\widetilde{{\epsilon}}_{w}(s,\mu)=\frac{{\epsilon}_{w}(s,\mu)}{(1+\xi(s,\mu))}.\end{split}$
(73)

Solving the former expression for the true signal $\mathcal{M}_{\ell}(s)$, the
shot-noise correction takes the same form analogous to Eq. (70)

$\mathcal{M}_{\ell}(s)=\frac{\mathcal{M}_{f,\ell}(s)\bar{m}_{f}^{2}-\widetilde{{\epsilon}}_{w,\ell}(s)}{(\bar{m}_{f}-{\epsilon}_{\bar{m}})^{2}},$
(74)

with

$\widetilde{{\epsilon}}_{w,\ell}(s)=\frac{(2\ell+1)}{2}\int_{-1}^{1}P_{\ell}(\mu)\widetilde{{\epsilon}}_{w}(s,\mu)\text{d}\mu,$
(75)

being the redefined shot-noise contribution to $w_{f}$ in redshift space. The
shot noise will still be, under the assumption of a mark that can be Taylor
expanded in powers of $\delta_{Rf}(\mathbf{x})$, a power series in $1/\bar{N}$
because the integral is additive. This means that the previously introduced
methodology of fitting polynomials to $w_{f}(\mathbf{r})$ and $\bar{m}_{f}$ is
applicable to redshift-space multipoles of the marked correlation function as
well. In the case of redshift-space multipoles of the weighted correlation
function the formulas are the same except the division by $1+\xi(s,\mu)$. It
is important to note here that, since the shot noise acquires a non-trivial
angle dependency and hence has to be corrected for in each multipole, also in
the marked power spectrum the shot noise will appear in higher multipoles.
This is in contrast to the ordinary power spectrum with constant shot noise
that will only affect the monopole.

In order to test our derived shot-noise correction in redshift space we
compute the analytical correction to the toy model in redshift space and
compare with our fitting methodology. As we have seen in Eq. (59) and Eq. (61)
the shot noise in the toy model is described by a linear polynomial in
$1/\bar{N}$. It has to be noted that in this case, analogous to the real-space
weighted correlation function, we correct $w_{f}$ with a linear factor of the
mean mark due to only one of the galaxies in each pair being actually
weighted. The results can be seen Figure 9, where we plot both the result from
the polynomial fit as well as the analytical correction for the monopole and
quadrupole. We find very good agreement between the two methods and the
relative difference in the monopole is below 1% over all scales up to
150$\,h^{-1}\,{\rm Mpc}$. For the quadrupole, the agreement is worse but still
within around 2% for most of the scales. There are specific spikes in the
relative difference caused by the quadrupole crossing zero at about
15$\,h^{-1}\,{\rm Mpc}$ and approaching zero on large scales.

Figure 9: Redshift-space multipoles of the toy model, measured in the GR
simulations of ELEPHANT, corrected for shot noise. The upper panel shows both
the monopole and quadrupole in different colours corrected via the polynomial
fit as solid lines and analytically corrected in dashed lines. In order for
better visualisation we offset the quadrupole by +1. Errorbars refer to the
mean standard deviation over 5 realisations. The lower panel shows the
relative difference between the analytic and polynomial correction in percent.

###  6.4 Limits of the shot-noise correction

It is important to assess where the shot-noise correction breaks down due to
the assumptions not being valid anymore. While an exhaustive investigation is
beyond the scope of this work, we discuss here several points in order to give
conservative limits on when it is safe to apply the proposed method. The most
crucial assumption comes from approximating the mark functional as a Taylor
expansion of the density contrast, as well as the derived power series in
$1/\bar{N}$ to be approximated by a low-order polynomial. A Taylor series of a
function $f(x)$ has a convergence radius $B$ for which the series converges to
the true function if it is evaluated inside the radius, that is, for $|x|\leq
B$. There is a straightforward way to compute the convergence radius of the
Taylor series of the $\tanh(x)$, for which techniques such as the ratio test
might not be applicable in certain cases. We can simply compute the closest
singularity to the point that we expand around ($x=0$), which gives us the
convergence radius. While the $\tanh(x)$ is non-singular on the axis of real
numbers, it has singularities in the complex plane. Since
$\tanh(x)=\sinh(x)/\cosh(x)$ singularities appear when $\cosh(x)=0$, which is
the case at $x=-b-1/a(1/2i\pi+2i\pi c)$ with $c\in\mathbb{Z}$, if we have both
a shift $b$ and factor $a$ as in the mark of Eq. (48). The closest singularity
to $x=0$ is therefore obtained by setting $c=0$, leading to a convergence
radius of $B=\sqrt{b^{2}+(\pi/(2a))^{2}}$. For the White mark the convergence
radius is simply given by $B=1+\rho_{*}$ in the case of a positive exponent
$p$ in Eq. (39).

In order to assess the effective validity of our correction based on the
convergence radius, we can reformulate it into a criterion involving a
measurable statistical quantity from the catalogues. For this, we follow a
similar approach as in Philcox et al. (2020). Starting from the mathematical
convergence criterion for the Taylor expansion $|\delta_{Rf}|\leq B$, we take
the density-weighted average222That is $\langle
g(\mathbf{x})\rangle_{\rho}=\frac{1}{V}\int
g(\mathbf{x})\frac{\rho_{f}(\mathbf{x})}{\bar{\rho}_{f}}\textrm{d}^{3}x$ of
the square on both sides. After taking the square-root we obtain

$\sqrt{\langle\delta_{Rf}^{2}\rangle_{\rho}}\leq B.$ (76)

The left-hand side quantity can be straightforwardly estimated by taking the
arithmetic mean of the square of $\delta_{Rf}$ at galaxy positions. For the
Covmos realisations using 64 grid cells per dimension
$\sqrt{\langle\delta_{Rf}^{2}\rangle_{\rho}}$ ranges from 0.61 to 1.01 over
the different levels of depletion (1.7% to 0.048%). In contrast, for the
ELEPHANT simulations of GR, $\sqrt{\langle\delta_{Rf}^{2}\rangle_{\rho}}$
takes values from 1.32 to 1.72 for no depletion down to 30%, respectively. By
choosing $(a,b)=(0.6,-0.5)$ for the $tanh$-mark we have a convergence radius
of $B\approx 2.67$, which satisfies the convergence criterion and thus we can
trust the Taylor expansion. Furthermore, the White mark with
$(\rho_{*},p)=(4.0,10.0)$ and $(\rho_{*},p)=(1.0,-1.0)$ have convergence radii
of $B=5$ and $B=\infty$ therefore they do also fulfil our convergence
criterion.

The differences between the catalogues used in this analysis affecting the
convergence criterion is illustrated in Figure 10 showing the density weighted
PDF of $\delta_{Rf}$. First thing that has to be noted is the fact that for
the black dashed curve, depicting the undepleted ELEPHANT PDF, there is not
depletion involved leading to more noise compared to the solid black curve for
Covmos. For the latter, we took the mean over 30 realisations of depletion
down to $1.7\%$ of the points of the full catalogue to match the undepleted
$\bar{N}$ of ELEPHANT. It is evident from the figure that the PDF for Covmos
is more peaked around zero with a less pronounced tail to high densities as
exhibited by the ELEPHANT PDF. That difference can be explained by the fact
that the Covmos realisations are meant to reproduce the distribution of dark
matter particles while the ELEPHANT catalogues are galaxies and hence contain
bias. Due to the higher skewness for the PDF in the ELEPHANT simulation, care
has to be taken to select appropriate marks not violating the convergence
criterion of the Taylor expansion.

Figure 10: Comparison of the PDF for the density contrast weighted by the
galaxie density as measured in Covmos and ELEPHANT (GR). Colours encode the
amount of depletion resulting in number densities per grid cell as indicated
in the legend. Dashed lines refer to the ELEPHANT simulation while solid lines
refer to the Covmos catalogues. We show the mean over 5 realisations and
depleted measurement were obtained by taking the mean over 30 realisations of
depletions.

The careful reader might have noticed that for the White mark with
$(\rho_{*},p)=(10^{-6},1.0)$ the radius of convergence is only around unity
and the above criterion would not be valid for the ELEPHANT simulation and
only partially valid for the Covmos catalogues. However, we do find good
recovery of the undepleted signal via applying the shot-noise correction to
this configuration in the Covmos catalogues. Moreover, we find good recovery
for the $\tanh$-mark with parameters $(a,b)=(10.6,-0.5)$ with a convergence
radius of only $B\approx 0.52$ therefore breaking the criterion completely.
This shows that even if the convergence criterion is not completely fulfilled,
it does not directly imply a failure of the shot-noise correction. Due to the
strong skewness of the distributions as shown in Figure 10, it might make
sense to look directly at the percentage of points with assigned densities
located inside the convergence radius. While this percentage ranges from 70%
down to 46% in the Covmos realisations for the $\tanh$-mark with
$(a,b)=(10.6,-0.5)$, for the White mark with $(\rho_{*},p)=(10^{-6},1.0)$ in
ELEPHANT it is still ranging from around 64% down to 51%, including at least
half the amount of points. Therefore, we would still trust the results of the
applied shot-noise correction in this configuration as we analyse them in the
next section. After all, even without a Taylor expansion, the shot-noise
behaviour might be well described with a polynomial since almost any function
can be locally fitted with a polynomial.

As mentioned earlier, next to the Taylor expansion, another limiting factor
are possible contributions of higher powers of $1/\bar{N}$, which are not
captured by e.g. 3rd-order polynomial fits. This is connected to the amplitude
of Taylor coefficients, which do only grow if the convergence radius is
smaller than 1. In general, the coefficients of the polynomial are expected to
decrease at some power as otherwise the shot-noise would look very noisy as a
function of $1/\bar{N}$. However, the coefficients might not be decreasing
fast enough such that a low-order polynomial is sufficient. This situation is
expected to worsen when the Taylor coefficients are growing. Non-accounted for
higher-order polynomials in the fit should manifest as a bias in the recovered
signal. A way to circumvent this issue is by extending the polynomials to
higher order, which, however, is not recommended in the case of only 8 data
points due to overfitting. Although we do find good performance of the shot-
noise correction for some marks with growing coefficients we also find worse
performance if e.g. the $b$ parameter is set to zero in the $\tanh$-mark. In
the latter case the Taylor coefficients are only non-zero for odd powers of
$\delta_{Rf}$ hence we suspect that higher-order correlators are more
important leading to higher-order polynomial contributions. While a thorough
assessment of the impact of Taylor expansion coefficients on the polynomial
fit is not done in this work we conservatively advise to use only marks with a
convergence radius larger than one to have decreasing Taylor coefficients.
Furthermore, the shift parameter $b$ in the $\tanh$-mark should be non-zero.
In Section 8 we elaborate further on the connection between polynomial
amplitudes and the smoothing induced by the MAS.

###  6.5 Shot noise in the White mark

With the previously developed methodology for correcting for shot-noise
effects in the weighted correlation function, we can evaluate the shot-noise
contributions to marks used in the literature so far. One particular widely-
used mark is the White mark given in Eq. (39), which has been used with the
parameter combinations $(\rho_{*},p)=(4.0,10.0)$ and $(\rho_{*},p)=(1.0,-1.0)$
in the work of Alam et al. (2021) and in the combination
$(\rho_{*},p)=(10^{-6},1.0)$ in the work of Hernández-Aguayo et al. (2018).
Both studies used a count-in-cells scheme to compute a density field on the
grid, which is the same as using a NGP MAS. In contrast, we used for most of
our analysis a PCS MAS, which is much wider in configuration space and
produces a smooth continuous and differentiable field.

Before we come to shot-noise-corrected differences between GR and MG in the
White mark, it is instructive to check how the shot-noise behaviour changes
when a different MAS is used. This can be investigated even without a shot-
noise correction by studying marked correlation functions for the undepleted
point set and the depleted point set down to 1.7% in the Covmos catalogues as
depicted in Figure 11. For this particular case, we use 60 grid cells per
dimension, instead of 64, to mimic closer the procedure in Alam et al. (2021)
and Hernández-Aguayo et al. (2018). By comparing the upper and middle panel in
Figure 11 it can be immediately seen that a more extended MAS, meaning a more
smooth density field, decreases the amplitude of the marked correlation
function. Intuitively this makes sense as we can expect a possibly stronger
small-scale correlation of the mark if the density field, used for the mark,
is less smooth. However, this comes at the price of stronger shot-noise
effects on small separations as can be seen in the lowest panel. Below scales
of around 10$\,h^{-1}\,{\rm Mpc}$, the relative difference between the
depleted and undepleted catalogue is larger if a NGP MAS is used as compared
to a PCS scheme. However, while the PCS scheme certainly reduces shot-noise
effects on the smallest scales it does show extended effects up to larger
scales, which can be seen as the solid lines (PCS) crossing the 5% limit at
larger scales compared to the dashed lines (NGP). We have shown such a feature
already in the toy model in Eq. (59) and (61), where the shot noise is
regulated to some extend by the MAS kernels. In our test on the Covmos
realisations, we find the shot-noise correction to give biased results when a
NGP MAS is used due to a very low smoothing size and shape of the NGP MAS. In
Section 8 we give a more extended discussion on this aspect and in the
following we refrain ourselves from using the NGP scheme in the White mark and
use the PCS MAS throughout, unless otherwise indicated.

Figure 11: Marked correlation functions for the White mark measured in the
Covmos catalogues. The upper and mid panel show the marked correlation
functions for different configurations of the parameters $(\rho_{*},p)$ using
NGP and PCS MAS, respectively. The dashed line refers to the undepleted case
whereas the solid line shows the measurement for a depleted catalogue down to
1.7%, which corresponds to the same mean density of points per grid cell as in
the undepleted ELEPHANT simulations. The lowest panel displays the relative
difference between the marked correlation function as measured in the full
data set and the depleted one. We used 60 grid cells per side length for the
density field and the vertical dashed line in grey corresponds to the side
length of a grid cell.

Let us now focus on the impact of shot noise on differences between MG and GR
in the ELEPHANT simulations. We present in Figure 12 the uncorrected marked
correlation function as measured in ELEPHANT for different configurations of
the White mark alongside with the shot-noise-corrected version. The MAS is
fixed to a PCS and we use again 60 grid cells per dimension. The shot noise
appears to only significantly affect the measurement below around 20
$h^{-}1$Mpc, but relative differences to the case with no correction can reach
up to 40% at these scales (middle panel). This is similar to what we reported
for the effect of shot noise in the Covmos simulations presented in Figure 11.
The shot noise has the smallest contributions for the configuration
$(\rho_{*},p)=(1.0,-1.0)$ where galaxies in high density regions get
upweighted as compared to the configurations with positive $p$ for which
galaxies in low-density regions are getting upweighted. In general, the effect
of shot noise as shown in Figure 12 is expected to be even stronger on the
smallest scales when a NGP MAS is used. In the studies of Alam et al. (2021)
and Hernández-Aguayo et al. (2018) relative differences between MG and GR are
particularly pronounced for scales below 20$\,h^{-1}\,{\rm Mpc}$ where we see
that shot noise has a significant effect. However, our findings do not nullify
those claimed differences as they might still be present after correcting both
GR and MG for shot noise. Rather the amplitude of the marked correlation
function itself will be different when correcting for shot noise, which is
particularly important for the modelling of marked statistics as e.g. done in
Aviles et al. (2020); Philcox et al. (2020, 2021). Indeed, as can be seen in
the lowest panel of Figure 12 the relative differences are largely unaffected
regarding a correction for shot noise. Although caution is advised as this
does not have to be universally valid for every mark. Furthermore, as we have
shown in the last paragraph about shot-noise effects in the Covmos catalogues
(see Figure 11), if a NGP MAS is used, the larger contributions at small
scales might impact relative differences stronger on those scales.

Figure 12: Shot noise in the White mark for the ELEPHANT suite. The different
configurations of parameters are colour coded and the solid line refers to the
corrected case while the dashed line refers to the uncorrected case. The upper
panel displays the marked correlation function both corrected and uncorrected.
The middle panel shows the relative difference between the corrected and the
uncorrected marked correlation function in percent. The MG model is fixed to
F4 in the upper and the middle panel. In the lower panel we show the relative
differences for the different configurations between the F4 model and GR in
percent. We used 60 grid cells per dimension and the vertical dashed line in
grey refers to the side length of one grid cell.

##  7 Results

###  7.1 Performance of marks based on large-scale environment

It has been already argued in Section 5 that, even though we do not have a
method for correcting the bias introduced by shot noise in this case, it is
instructive to assess the performance of marks based on the environmental
classification and tidal field/torque regarding distinguishing GR from MG.

In Figure 13 we present the SNR as defined in Eq. (37) for marks based on the
environmental classification introduced in Section 5. While marks like the
Void and Wall correlation function only exhibit significant differences on
small scales below 10$\,h^{-1}\,{\rm Mpc}$ and 40$\,h^{-1}\,{\rm Mpc}$,
respectively, the marks introducing anti-correlation produce significant
differences up to scales of around 80$\,h^{-1}\,{\rm Mpc}$. Particularly the
weaker modifications of gravity like F5 and F6 appear to profit from anti-
correlation as the $\textrm{Void}_{\textrm{AC}}$ mark has an SNR of 3 up to
around 40$\,h^{-1}\,{\rm Mpc}$ and the $\textrm{Wall}_{\textrm{AC}}$ mark
shows a similar behaviour for F5. Furthermore N1 can be well distinguished
with the $\textrm{Wall}_{\textrm{AC}}$ mark up to scales of around
60$\,h^{-1}\,{\rm Mpc}$. The $\textrm{Void}_{\textrm{LEM}}$ mark performs well
on the F4 simulations up to scales around $\sim$60$\,h^{-1}\,{\rm Mpc}$ but
the SNR is only around 3 from 30$\,h^{-1}\,{\rm Mpc}$ onwards. N1 and F5 show
a SNR of around 3 or larger only for scales up to around 20$\,h^{-1}\,{\rm
Mpc}$ and 30$\,h^{-1}\,{\rm Mpc}$, respectively.

The situation is different for marks using the tidal field/torque as depicted
in Figure 14. Using the tidal field as it is does not seem to yield any
significant difference for the investigated MG theories. Interestingly, by
taking just a linear function of the tidal torque, we can report significant
differences for F6 up to scales of $\sim$60$\,h^{-1}\,{\rm Mpc}$. As the tidal
torque is small for symmetric large-scale environments we basically upweight
filaments and walls by taking the tidal torque as a mark. Since there are many
galaxies located in walls this appears to compensate for the fact that MG
effects are expected to be more screened in walls compared to voids and lead
to significant differences to GR when used as a mark.

Although shot noise is expected to decrease at higher scales and that they
might not affect differences between two measurements strongly, there is no
guarantee that the observed significant differences seen in Figure 13 and 14
are still present after a correction for shot noise. In principle we could
assess the overall impact of shot noise on these marks by looking at
differences in the corresponding measurements for the Covmos catalogues as we
have done it in Figure 11 for the White mark. However, due to the peculiar way
the Covmos catalogues have been set up we do not expect that they also
represent large-scale structure environments in a realistic way. Furthermore,
qualitative differences are not enough to assess precisely how the SNR will
change after a proper correction. We relegate, therefore, a thorough
investigation of shot-noise effects for these marks to future work for which
high resolution full N-body simulations of MG are necessary.

Figure 13: Summary of the SNR measured in the ELEPHANT suite for the marked
correlation functions using marks based on the large-scale environment as
introduced in Section 5. Colours refer to MG simulations and the panels show
different marks as indicated on the labels. The horizontal dashed lines in
black indicate a SNR of $\pm 3$.  Figure 14: Same as Figure 13 but for marks
based on the tidal field and tidal torque.

###  7.2 Performance of the White mark

In Figure 15 we present the performance of the White mark, corrected for shot
noise, to be compared with the other marks to follow. We used 64 grid cells
per dimension and a PCS MAS to obtain the density field on the grid and the
parameters were fixed to $(\rho_{*},p)=(10^{-6},1.0)$. As described in Section
6, we are using third and second order polynomials to fit the shot-noise
dependency of $w_{f}(\mathbf{r})$ and $\bar{m}_{f}$, respectively. Overall,
the amplitude of the marked correlation function does not differ much from
unity implying a minor impact of the mark. Although we can see differences
from unity up to scales of around 70$\,h^{-1}\,{\rm Mpc}$, the signal is very
similar between GR and MG for most of the scales except below
20$\,h^{-1}\,{\rm Mpc}$. This can be deduced from the lower plot as well where
we show the SNR, directly quantifying the difference between GR and MG. The
SNR lays inside the $3\sigma$ region with only occasional peaks outside that
range, which can be accounted to sample variance. Only for F4 we can report
significant differences for the first four bins in $r$ ranging up to $\sim
20$$\,h^{-1}\,{\rm Mpc}$. The SNR for the case $(\rho_{*},p)=(4.0,10.0)$,
although not shown here, exhibit similar overall structure as the presented
case. This makes the White mark with those configurations in real space not
particularly powerful in distinguishing MG from GR as possible differences
only show up at very low scales.

Figure 15: Performance of the White mark for fixed parameters
$(\rho_{*},p)=(10^{-6},1.0)$. The upper panel displays the mean marked
correlation function taken over 5 realisations together with errorbars
estimated as the mean standard deviation. Colours refer to the different
gravity simulations. The lower panel shows the SNR as introduced in Eq. (37)
and the blue shaded region refers to the error of a single measurements
divided by the error of the mean difference (see Eq. (38)) for F4.

In Figure 16 we present the monopole and quadrupole of the White mark in
redshift space with the same parameter configuration as before, meaning
$(\rho_{*},p)=(10^{-6},1.0)$. The left panel depicts the monopole, exhibiting
a very similar amplitude and shape as the real space measurements in Figure
15. The monopole does converge to unity at higher scales which is expected as,
even in redshift space, the marks should become uncorrelated at high scales
leading to $\mathcal{M}(s,\mu)=1$. The SNR in the monopole shows no
significant difference over all scales except for F4 below $20\,h^{-1}\,{\rm
Mpc}$ that we have also seen in real space. The quadrupole (right panel) has
much weaker signal compared to the monopole and converges to zero at high
scales. This can be explained with the same reasoning as to why the monopole
converges to one. $\mathcal{M}(s,\mu)$ is independent of $\mu$ on large scales
and therefore does not possess higher-order multipoles. Even though the
amplitude is very small we nevertheless see interesting differences for N1 in
the SNR on moderate to high scales. However, the SNR is barely above 3 and the
bin-to-bin variance is fairly high. Lastly, we did not find significant
differences for the configuration $(\rho_{*},p)=(4.0,10.0)$ over extended
scales. Similar to what we have seen in real space, the White mark in these
configurations is overall not really promising in redshift space.

Figure 16: Redshift-space monopole (left side) and quadrupole (right side) of
the marked correlation function for the White mark with parameters fixed to
$(\rho_{*},p)=(10^{-6},1.0)$ and corrected for shot noise. Upper panels show
the multipoles themselves with colours referring to different gravity
simulations. Displayed is the mean over five realisations with respective mean
standard deviations as the error bars. For the monopole the horizontal dashed
line in black marks an amplitude of 1 while for the quadrupole it marks an
amplitude of 0. Lower panels show the signal to noise ratio, as in Eq. (37),
with according colour for the different MG models. Shaded regions refer to the
error of a single measurement divided by the mean error of the difference as
introduced in Eq. (38). Dashed black lines in the lower panels indicate a SNR
of $\pm 3$.

###  7.3 Performance of the $\tanh$-mark

In Figure 17 we present both the marked correlation function before and after
shot-noise correction in the left and right panel, respectively, for the
$\tanh$-mark as introduced in Section 5.2. The configuration is set to
$(a,b)=(0.6,-0.5)$ for which we showed in Section 6 that our correction
algorithm can be safely applied. We use polynomials of order three and two for
correcting $w_{f}(\mathbf{r})$ and $\bar{m}_{f}$, respectively. The marked
correlation function is largely featureless and converges to 1 on large scales
regardless if a correction for shot noise is applied or not. This convergence
is also present in the White mark in Figure 15 and is due to the mark getting
uncorrelated at large scales. It is striking how the amplitudes differ on
smaller scales between the uncorrected (left panel) and corrected case (right
panel) underlining again the importance of applying a shot-noise correction in
order to measure correct amplitudes. In the lower panels we are showing the
SNR as defined in Eq. (37). The general trend of the SNR for the different MG
models appears to be similar regardless if a shot-noise correction is applied
or not. However, e.g. F4 shows much larger significance on small scales in the
case of no correction. Most importantly, the $\tanh$-mark in this
configuration leads to significant differences for F6 in the corrected case,
up to scales of around 80$\,h^{-1}\,{\rm Mpc}$. Worth to notice here that
significant differences are also found for F4 and N5 but only for scales
smaller than roughly 30 and 40$\,h^{-1}\,{\rm Mpc}$, respectively. While
differences on these scales might be sufficient to be grasped by theoretical
models appropriately, it has yet to be tested as modelling becomes
increasingly more challenging at small scales. In Section 8 we are discussing
our results on the $\tanh$-mark in the light of current constraints in the
literature on the $f_{R0}$ parameter.

Since our SNR as defined in Eq. (37) computes the error over 5 realisations
any significant differences can only be claimed for a volume of 5
realisations. It is therefore instructive to elaborate on how the difference
compares to the error of a single measurement as indicated by the shaded area
in the plot. It appears that particularly at higher scales the difference is
of the same size as the error itself rendering a detection at the current
volume of around $1\,h^{-3}\textrm{Gpc}^{3}$ with only one measurement at hand
impossible. However, this could be alleviated by considering simulations in
larger volumes. Assuming that the error of the single measurement is Gaussian,
it scales like $\propto 1/\sqrt{V}$, where $V$ is the volume of the
survey/simulation. An increase in volume by a factor of 9 for F6 would be
necessary in order to detect the difference with just one measurement at
scales between $60\,h^{-1}\,{\rm Mpc}$ and $80\,h^{-1}\,{\rm Mpc}$. This
increase in volume translates in a larger box side length by a factor of
around 2.1. On smaller scales, below around $40\,h^{-1}\,{\rm Mpc}$, the
reported differences would be significant even with only a single measurement.

Figure 17: Marked correlation function for the $\tanh$-mark with
$(a,b)=(0.6,-0.5)$. The left panels depict the case where we do not apply any
correction for shot noise while in the right panels we do apply our shot-noise
correction methodology as described in Section 6. Upper panels show the marked
correlation functions with colours encoding the different MG models. Displayed
is the mean over 5 realisations and errorbars are obtained by taking the mean
standard deviation. The lower panel displays the SNR where the black-dashed
line indicates a difference of $\pm 3$. Shaded areas mark the error of a
single measurement divided by the mean error of the difference. The black
dashed line in the upper panels indicates an amplitude of 1.

In order to better compare our results with the literature where often only
relative differences between GR and MG are reported (Hernández-Aguayo et al.,
2018; Armijo et al., 2018; Alam et al., 2021) we show a corresponding plot in
Figure 18. Large relative differences beyond 5% are reached for F4 only on
smaller scales below around 20$\,h^{-1}\,{\rm Mpc}$ while N5 shows larger
relative differences all the way up to around 50$\,h^{-1}\,{\rm Mpc}$. This is
topped by F6 exhibiting relative differences above 15% over almost all scales
decreased only above $\sim$80$\,h^{-1}\,{\rm Mpc}$. However, the relative
error for the GR simulation increases towards those scales rendering the
detection with only one realisation at those high scales unfeasible. As shown
in Figure 17 more volume is needed to shrink the uncertainties to a level at
which the large reported differences between GR and F6 at high scales can be
taken advantage of. These results can be compared with the Fig. 5 in the work
of Hernández-Aguayo et al. (2018) and Fig. 15 in Alam et al. (2021) where
marks based on the Newtonian gravitational potential were used and appear to
be the most performant in terms of relative difference between GR and MG.
While certainly performing very good for F4 and to some extend for F5, we can
report much larger relative differences for F6, reaching up to higher scales
if a $\tanh$-mark in the configuration $(a,b)=(0.6,-0.5)$ is used.
Furthermore, our mark is very easy to compute and does not need information
from halos as is the case if the Newtonian potential is to be computed in the
way defined in the mentioned studies.

Figure 18: Relative differences between GR and MG for the $\tanh$-mark with
parameters fixed to $(a,b)=(0.6,-0.5)$. We show the mean relative difference
over five realisations and the shaded area corresponds to the relative
standard deviation (relative error of single measurement) for the GR
realisations. Grey dashed lines indicate relative differences of $\pm$5%.

Finally, in Figure 19 we present the shot-noise-corrected monopole and
quadrupole of the marked correlation function in redshift space for the
$\tanh$-mark with parameters fixed to $(a,b)=(0.6,-0.5)$. In general, compared
to the White mark in Figure 16, the amplitude of both the monopole and the
quadrupole is much larger but the large-scale behaviour is very similar.
Looking at the monopole in the left panel, similarities with the real space
measurements in Figure 17 are striking both in the shape of the measurements
as well as the SNR. However, N5 has a reduced SNR in redshift space while
conversely N1 has now significant differences at scales lower than
40$\,h^{-1}\,{\rm Mpc}$. F6 and F4 look largely the same as in real space,
most importantly the former still showing significant differences up to around
80$\,h^{-1}\,{\rm Mpc}$. The shaded region seems to increase in size for F6
making an even larger volume necessary to detect the difference with a single
observation only. The amplitude in the quadrupole in the upper right panel is
smaller and also the SNR in the lower right panel stays within 3$\sigma$
rendering the quadrupole not suitable as a statistic to detect MG with this
mark.

Figure 19: Redshift-space monopole (left side) and quadrupole (right side) of
the marked correlation function for the $\tanh$-mark with parameters fixed to
$(a,b)=(0.6,-0.5)$ and corrected for shot noise. Upper panels show the
multipoles themselves with colours referring to different gravity simulations.
Displayed is the mean over five realisations with respective mean standard
deviations as the errorbars. For the monopole the horizontal dashed line in
black indicates an amplitude of 1 and for the quadrupole an amplitude of 0.
Lower panels show the SNR, as in Eq. (37), with according colours for the
different MG models. Shaded regions refer to the error of a single measurement
divided by the mean error of the difference and dashed lines in the lower
panels indicate a SNR of $\pm 3$.

##  8 Discussion

We have found the $\tanh$-mark to be particularly promising regarding
distinguishing MG from GR. Of course, studying possible differences between
$f(R)$ theories and GR has to be done in the light of constraints on $f_{R0}$
in the literature. A somewhat older compilation of constraints can be found in
Tab. 1 in the work of Lombriser (2014), where the strongest limit comes from
dwarf galaxies and the solar system imposing $|f_{R0}|\leq 10^{-7}-10^{-6}$.
In a more recent analysis, Liu et al. (2021) found similar limits using Fisher
forecasts on cluster abundances and galaxy clustering. Even tighter
constraints, $f_{R0}<1.4\times 10^{-8}$, are found in the work of Desmond &
Ferreira (2020) using galaxy morphology. Although F6 might not be a viable MG
theory after all and has to be replaced with weaker modifications like F7 or
F8, finding significant differences for F6 makes the $\tanh$-mark promising to
distinguish also weaker models.

In Section 6 we presented a robust technique to correct for shot-noise effects
for general marks, where the mark function can be expressed as a Taylor
expansion in powers of the density contrast. Since the error on the
measurements plays a crucial role in our performance metric in Eq. (37), it is
instructive to discuss how the relative error of the measurements is impacted
when the shot-noise correction is applied. In particular due to the
approximations made in estimating the shot-noise-corrected signal, additional
uncertainties might be introduced. In Figure 20 we present the relative error
for the undepleted case and the corrected case, both for the $\tanh$-mark with
$(a,b)=(0.6,-0.5)$ and toy model. Since we can capture the shot-noise
behaviour in the toy model very accurately, the relative error stays almost
the same and does not significantly change. However, in the case of the