# Testing modified gravity theories with marked statistics

Alejandro Aviles  Consejo Nacional de Ciencia y Tecnología,
Av. Insurgentes Sur 1582, Colonia Crédito Constructor, Del. Benito Juárez,
03940, Ciudad de México, México,
and
Departamento de Física, Instituto Nacional de Investigaciones Nucleares
Apartado Postal 18-1027, Col. Escandón, Ciudad de México,11801, México.
E-mail<EMAIL_ADDRESS>
###### Abstract

In the last two decades, Modified Gravity (MG) models have been proposed to
explain the accelerated expansion of the Universe. However, one of the main
difficulties these theories face is that they must reduce to General
Relativity (GR) at sufficiently high energy densities, such as those found in
the solar system. To achieve this, MG theories typically employ so-called
screening mechanisms: nonlinear effects that bring them to GR at the
appropriate limits. For this reason, low-energy regions where the screenings
do not operate efficiently, such as cosmic voids, are identified as ideal
laboratories for testing GR. Hence, the use of marked statistics that up-
weight low energy densities have been proposed for being implemented with data
from future galaxy surveys. In this proceeding note, we show how to construct
theoretical templates for such statistics and test their accuracy with the use
of N-body simulations.

###### keywords:

Modified Gravity; Cosmic Large Scale Structure; Testing General Relativity.

##  1 Modified gravity in Cosmology: a very brief introduction

The standard model of Cosmology, the so-called $\Lambda$CDM, is the most
successful theory to date for describing the large-scale behaviour of our
Universe. The outstanding way this model fits the cosmic background radiation
(CMB) anisotropies measurements from the Planck satellite[1]—and prior to that
from the WMAP and COBE experiments—leaves small room for modifications of it.
However, among other drawbacks in Cosmology, nowadays cosmic acceleration is
still not well understood, and in models that modify the background history of
the Universe at late times, the CMB differs mainly because the journey of
photons from the last-scattering surface to us is affected. This is a
secondary effect called the integrated Sachs-Wolfe effect and its influence
over the CMB is mainly localized at small multipoles, where the error on
measurements is nonetheless large because of the cosmic variance—the fact that
we have only one universe and hence a small number of long wave-length modes
to average over. Such a liberty has led to a lot of proposals (literally
hundreds of them) as alternatives to the cosmological constant to explain the
speeding up of the expansion of the Universe. Unfortunately, there is much
more literature proposing new gravity models or that makes simple tests to
them (the vast majority focusing only on the homogeneous and isotropic
cosmology; I will come back to this in a moment), than literature that explore
methods to probe the clustering due to gravity at cosmological scales. Of
course, this is a highly unbalanced situation and observational cosmologists
cannot test each of the proposed models. These alternatives can be categorized
into two big branches: dark energy (DE) and Modified Gravity (MG); though, one
should be careful since there is not a well defined boundary between them. In
this proceedings note, we will focus on the latter assuming well defined, and
to the particular case of chameleon fields[2, 3] although the techniques
described here have been applied to other kinds of theories.

Perhaps the most difficult aspect of modifying gravity is the following.
Observations indicate that the true model of the universe should be close to
$\Lambda$CDM, at least from the epoch of the primordial (“Big-Bang”)
Nucleosynthesis until now. In particular the expansion history should be close
to that of $\Lambda$CDM, perhaps with small differences with the aim of, for
example, resolving the Hubble tension.[4] However, as a rule of thumb, MG
models that slightly depart from $\Lambda$CDM at the background level, become
quite different at the perturbative level.[5, 6] That is, two models that have
almost the same expansion history, typically predict a completely different
clustering. This is because in MG the two scalar gravitational potentials of
the metric, thinking of in Newtonian gauge, are different even in the absence
of anisotropic stresses, which makes GR a very special theory and not easy to
mimic. This would mean that the gravitational potential appearing in the
geodesic potential is not the same to the gravitational potential sourced by
the energy density fields in the Poisson equation; and what is particularly
important to this discussion is that the homogeneous and isotropic expansion
is completely insensitive to this. For example, take a look at figure 1, where
we show the power spectrum of matter density fluctuations in GR and in the F5
gravitational model (we will define what is F5 in appendix Acknowledgments).
The F5 gravitational model is indistinguishable for any practical purpose from
GR at the background level, but their predicted matter power spectra differ by
a lot. This is the reason why it does not make too much sense to test MG
models at the background level only, for example using supernovae data only.
However, for the sake of fairness we must say there are exceptions to this
rule, for example the recently proposed No-Slip gravity,[7] tailored with the
specific purpose that the two scalar gravitational potentials become equal.

When gravity is modified, in particular in the infrared, one should be careful
to not spoil the well-tested regimes of validity of GR, as for example in the
solar system or regions with much larger densities,111It is perhaps surprising
to consider the solar system a high energy density region, but it is very
large compared to the mean background energy density $3H_{0}^{2}/8\pi G\sim
10^{-29}\,\text{g}/\text{cm}^{3}$. The latter is the one that is important at
the largest cosmological scales, and hence drives the background history of
the Universe.  but at the same time provide the accelerated expansion of the
Universe. In order to deal with this unquestionable observational
constriction, MG models often invoke mechanisms that yield the theory to GR in
the appropriate limits. These are generically called screening mechanisms:
non-linear effects that drive MG theories to GR when either the environmental
energy density is large, or there are large gravitational tidal fields, among
other scenarios. The most popular is the chameleon mechanism,[2, 3] in which a
scalar field mediates a universal fifth force, but its mass depends on the
environmental density where it resides: as larger is the ambient density, the
lighter is the mass of the scalar field and smaller the range of the force.
For example, in successful theories, and for human-like densities ($\sim$ 1
gram per cubic centimeter) the range of the fifth-force becomes so
ridiculously small that it becomes effectively invisible to any experiment. On
the other hand, at extremely low ambient densities the range of the fifth-
force becomes quite large, with cosmological scope.

The most studied MG theory in Cosmology, one can affirm with high confidence,
is the $f(R)$ gravity,[8, 9, 10] and in particular the Hu-Sawicki model.[11]
Although, these theories were first understood as higher-order gravity
theories, it was soon recognized that under a few assumptions they can be
described as second order differential equations for the metric plus a scalar
field[12, 13] which propagates only an additional massive mode. In this
proceeding note we will lead with the Hu-Sawicki model to exemplify our
analytical results and to take advantage of the fact that we have good, state-
of-the-part N-body simulations for them—the Extended LEnsing PHysics using
ANalaytic ray Tracing Elephant cosmological suite of simulations.[14, 15] Even
though $f(R)$ is not a representative model of the vast universe of MG models
(or more narrowly of Hordenski theories), the techniques reviewed here can be
applied to different models and can be beneficial to study MG that rely on
screening mechanisms. The material presented in this note is extracted mainly
from Refs. 16, 17, 18, 19; we refer the reader to these works for further
details.

Figure 1: Matter power spectrum in the models GR (blue line) and Hu-Sawicki F5
(red lines). The latter is further shown with and without screenings with
solid and dashed lines, respectively. To generate this analytical results we
use the Standard Perturbation Theory with effective field theory counterterms
and infrared-resummations of Refs. 20, 21. The mock data is extracted from the
Elephant simulations.

##  2 The importance of up-weight low energy density regions

In figure 1 we show the chameleon screening in action. The black line is the
analytical power spectrum of the $\Lambda$CDM. The red coloured lines show the
matter power spectrum in the Hu-Sawicki F5 model with screenings (solid line)
and without screenings (dashed line), together with data from the Elephant
simulations. We note the screening mechanism drives the theory to GR, more
evident at small scales (high wave-numbers $k$), while without the screenings
the differences between the two models are larger. This behaviour raises the
question if the power spectrum, and its counterpart in configuration space,
the correlation function, are the most clever way to test the gravitational
theories from its cosmological clustering. After all, by construction these
statistics give more weight to high density regions, since they sum over all
tracers evenly (either particles, halos, galaxies…), and tracers tend to clump
together by gravitational attraction. Meanwhile, cosmological voids are
underrepresented by the power spectrum—as well as by the correlation function
and their N-points generalizations—even though they are the largest regions in
our universe, and as such, they contain a significant share of the total
energy budget. But clumpy regions of space, being denser, behave very close to
GR due to the screenings mechanism. Hence, the standard statistics used in
Cosmology may not be the best option to test infrared modifications of GR.

With the above ideas in mind, M. White proposed in Ref. 16 to use marked
correlation functions222Marked statistics have a long history;[22] they have
been used to assign properties to objects, such as the luminosity, color, and
morphology of galaxies,[23, 24] and to break degeneracies between Halo
Occupation Distribution (HOD; a semi-analytical method to populate halos from
N-body simulations with galaxies) and cosmologies.[25] to test gravity at
cosmological scales. The proposal consists in weighting the tracer field
density fluctuations with a function (the mark) whose value depends on the
environmental density, such that objects residing in low density regions
become up-weighted in 2-point statistics, while objects residing in high
density become down-weighted. In this way, statistics of such marked density
fields will probe low density regions where gravitational fifth-force
screenings are inefficient and the effects of hypothetical MG theories would
be more pronounced.

Figure 2: 1-dimensional cartoon for the construction of the marked density
field.

The construction of a marked density field $\delta_{M}({\text{\bf x}},t)$ is a
three step process, depicted in figure 2. It goes as follows,

  1. 1.

One starts with the matter fluctuations field $\delta({\text{\bf x}},t)$ and
convolve it with a window function $W_{R}({\text{\bf x}})$ of width $R$ (a
Gaussian is used in this work) to obtain the smoothed density

$\delta_{R}({\text{\bf x}},t)=\int d^{3}y\;W_{R}(|{\text{\bf x}}-{\text{\bf
y}}|)\delta({\text{\bf y}},t).$ (1)

The scale $R$ is chosen from the beginning and it determines the size of the
region which is used to define the environmental density of the tracers. In
this sense, it is good to have sufficiently large values for $R$, typically
comparable to, or a little bit smaller than, the size of cosmic voids. Here we
use $R=10\,h^{-1}\,\text{Mpc}$.

  2. 2.

Construct the mark function $m(\delta_{R}({\text{\bf x}}))$. We choose a
function that decays with the environmental (or smoothed) density, such that
low density regions are up-weighted, as required. For example, a simple
exponential $m=\exp[-\alpha(1+\delta_{R}({\text{\bf x}}))]$ can do the job.
However it is better to have more free parameters such that we have more
liberty to find an optimal mark for a particular MG model. The most used
weight is the White-mark[16]

$m(\delta({\text{\bf
x}}),R)=\left(\frac{1+\delta_{*}}{1+\delta_{*}+\delta_{R}({\text{\bf
x}})}\right)^{p},$ (2)

with $p>0$ and $\delta_{*}>0$ dimensional parameters. It was find with
numerical simulations that the choices $p=10$ and $\delta_{*}=4$ maximize both
the Fisher information and the signal-to-noise.[26] However, these parameters
are not good for perturbative analysis, so we will use the parameters $p=7$
and $\delta_{*}=10$, unless otherwise is stated. Other authors working with
simulations have used different parameters and mark functions.[26, 27, 28, 18]

  3. 3.

Weight the density field of tracers $1+\delta_{X}({\text{\bf x}})$ with the
mark function $m(\delta({\text{\bf x}}),R)$ to obtain the marked density field

$1+\delta_{M}({\text{\bf x}},t)=\big{(}1+\delta_{X}({\text{\bf
x}},t)\big{)}m[\delta_{R}({\text{\bf x}})].$ (3)

The density field can be that of any tracer of the underlying matter density,
it can even be the matter density field itself. But we observe tracers, as
galaxies for example, in the sky. These are related through a biasing function
$F$ to the dark matter as

$1+\delta_{X}({\text{\bf x}})=F_{\text{\bf x}}[\delta({\text{\bf
x}}),\nabla^{2}\delta({\text{\bf x}}),...;{\text{\bf x}}]$ (4)

where the stochasticity of this relation is emphasized by writing x as an
additional argument. The dependence on the non-local curvature operator
$\nabla^{2}\delta({\text{\bf x}})$ is expected to be small in GR, but a well
defined bias expansion within chameleon theories must contain higher-order
derivative operators.[29, 30]

##  3 Marked correlation function

A marked correlation function (mCF) [22, 31, 23, 24, 32, 25, 16] is defined as
the sum of pairs of objects ($i,j$) separated by a distance $r$, and weighted
by the ratio of the mark function value to the mean mark $m_{i}/\bar{m}$ at
each point and divided by the number of pairs $n(r)$,

$\mathcal{M}(r)=\sum_{i,j|r_{ij}=r}\frac{m_{i}m_{j}}{n(r)\bar{m}^{2}},$ (5)

so it quantifies the deviation of a pair of objects to have the mean mark
$\bar{m}$ subjected to be located at a separation $r$. By construction it can
be written as

$\mathcal{M}(r)=\frac{1+W(r)}{1+\xi(r)}$ (6)

where $\xi(r)$ is the standard 2-point correlation function of tracers

$1+\xi(r)=\left\langle\big{(}1+\delta_{X}({\text{\bf x}}+{\text{\bf
r}})\big{)}\big{(}1+\delta_{X}({\text{\bf x}})\big{)}\right\rangle,$ (7)

and $W(r)$ is the 2-point correlation function of marked fields:

$1+W(r)=\frac{1}{\bar{m}^{2}}\left\langle\big{(}1+\delta_{M}({\text{\bf
x}}+{\text{\bf r}})\big{)}\big{(}1+\delta_{M}({\text{\bf
x}})\big{)}\right\rangle.$ (8)

At large scales both $\xi$ and $W$ are small and one gets

$\mathcal{M}(r)=1+W(r)-\xi(r).$ (9)

Now, one may show that $\xi(r)\in W(r)$,[17] hence the mark function becomes a
2-point statistic of the clustering of marks with the clustering of objects
effectively factorized out. The presence of the mean mark is not necessary to
estimate the mCF, but it plays an important role in the analytical
renormalization of the theory.

##  4 Perturbative treatment of the marked correlation function

It is well known that a good modeling of N-point statistics of density fields
in configuration space requires the use of Lagrangian Perturbation Theory
(LPT),[33] since the spread and degradation of the BAO oscillations is due to
large scale bulk flows, long wave-length Lagrangian displacement fields that
are well described in the Lagrangian framework, even at the linear order.[34]
However the derivation of the marked correlation function within the
Lagrangian framework is quite large and technical, so we prefer to omit it in
this note; the interested reader can find it in Ref. 17, and later on we will
display the main results. However, in this section we develop the Standard
Perturbation Theory (SPT) approach to the marked correlation function, which
is easier, but yet shows the main ingredients and insights of the theory.

We start by expanding the mark in a power series of $\delta_{R}$

$m(\delta_{R};C_{i})=C_{0}+C_{1}\delta_{R}+\frac{1}{2}C_{2}\delta^{2}_{R}+\cdots$
(10)

with $C_{i}=m^{(i)}[0]$, the $i$-th derivative of $m[\delta_{R}]$ evaluated at
$\delta_{R}=0$. The fewer parameters are needed to model the mark, the better
the convergence will be. We notice that $C_{1}<0$ enhances low density
regions.

For the sake of simplicity we will assume a local relation between matter and
tracer overdensities

$1+\delta_{X}({\text{\bf x}})=c_{0}+c_{1}\delta({\text{\bf
x}})+\frac{c_{2}}{2}\delta^{2}({\text{\bf x}})+\cdots,$ (11)

with $c_{i}=\partial F_{\text{\bf
x}}^{(i)}[\delta_{R}=0]/\partial\,\delta_{R}$ the bare bias parameters and
$F_{\text{\bf x}}$ given by eq. (4). A Taylor expansion is an unnatural and
restrictive assumption for biased tracers, but being pragmatic we note that it
works at the level of precision required by our simulations.

The mean mark is given by the mark weighted by the tracer density field

$\displaystyle\bar{m}$ $\displaystyle=\langle m[\delta_{R}({\text{\bf
x}})](1+\delta_{X}({\text{\bf x}}))\rangle$
$\displaystyle=\Big{\langle}\Big{(}C_{0}+C_{1}\delta_{R}({\text{\bf
x}})+\frac{C_{2}}{2}\delta^{2}_{R}({\text{\bf
x}})+\cdots\Big{)}\Big{(}c_{0}+c_{1}\delta({\text{\bf
x}})+\frac{c_{2}}{2}\delta^{2}({\text{\bf x}})+\cdots\Big{)}\Big{\rangle}$
$\displaystyle=c_{0}C_{0}+c_{1}C_{1}\sigma^{2}_{R}+\frac{1}{2}c_{0}C_{2}\sigma^{2}_{RR}+\frac{1}{2}C_{0}c_{2}\sigma^{2}+\cdots,$
(12)

with zero-lag correlators defined as

$\sigma^{2}\equiv\langle(\delta(0))^{2}\rangle,\quad\sigma^{2}_{R}\equiv\langle\delta_{R}(0)\delta(0)\rangle,\quad\sigma^{2}_{RR}\equiv\langle(\delta_{R}(0))^{2}\rangle.$
(13)

We will use also the correlation and cross-correlation functions

$\displaystyle\xi(r)$ $\displaystyle\equiv\langle\delta({\text{\bf
x}})\delta({\text{\bf x}}+{\text{\bf r}})\rangle,$ (14)
$\displaystyle\xi_{R}(r)$ $\displaystyle\equiv\langle\delta_{R}({\text{\bf
x}})\delta({\text{\bf x}}+{\text{\bf r}})\rangle,$ (15)
$\displaystyle\xi_{RR}(r)$ $\displaystyle\equiv\langle\delta_{R}({\text{\bf
x}})\delta_{R}({\text{\bf x}}+{\text{\bf r}})\rangle.$ (16)

We have to write eq. (4) in terms of renormalized bias parameters,[35]

$b_{n}^{E}=\int\frac{d\lambda}{2\pi}e^{-\lambda^{2}\sigma^{2}/2}\tilde{F_{\text{\bf
x}}}(\lambda)(i\lambda)^{n},$ (17)

where $\tilde{F}_{\text{\bf x}}(\lambda)$ is the Fourier transform of
$F_{\text{\bf x}}$ with spectral parameter $\lambda$ (dual to $\delta$). We
are using the label “$E$” to distinguish Eulerian from Lagrangian biases.
Equation (17) leads to the relation between bare bias $c_{n}$ and renormalized
bias $b_{n}^{E}$ parameters:[36, 37]

$b_{n}^{E}=\sum_{k=0}^{\infty}\frac{\sigma^{2k}}{2^{k}k!}c_{n+2k}.$ (18)

Analogously, we introduce the “resummed” expansion parameters $B_{n}$ as

$B_{n}=\frac{B_{n}^{*}}{B_{0}^{*}}\qquad\text{with}\qquad
B_{n}^{*}=\int\frac{d\Lambda}{2\pi}e^{-\Lambda^{2}\sigma^{2}_{RR}/2}\tilde{m}(\Lambda)(i\Lambda)^{n},$
(19)

where $\tilde{m}(\Lambda)$ is the Fourier transform of $m(\delta_{R})$, and
$\Lambda$ is a spectral parameter, dual to $\delta_{R}$. One easily finds[17]
the following relation between the expansion parameters $B_{n}$ and the Taylor
coefficients of the mark function $C_{n}$

$B_{n}(C_{n},\sigma_{RR}^{2})=\frac{\sum_{k=0}^{\infty}C_{n+2k}\sigma^{2k}_{RR}/(2^{k}k!)}{\sum_{k=0}^{\infty}C_{2k}\sigma^{2k}_{RR}/(2^{k}k!)}.$
(20)

Inserting the $b_{n}^{E}$ and $B_{n}$ parameters in eq. (4) we get

$\bar{m}=B_{0}^{*}\big{[}1+b_{1}^{E}B_{1}\sigma^{2}_{R}+\cdots\big{]}.$ (21)

We compute now

$\displaystyle\bar{m}^{2}(1+W(r))=\langle m[\delta_{R}({\text{\bf
x}}_{1})]\big{(}1+\delta_{X}({\text{\bf x}}_{1})\big{)}m[\delta_{R}({\text{\bf
x}}_{2})]\big{(}1+\delta_{X}({\text{\bf x}}_{2})\big{)}\rangle$
$\displaystyle=\int\frac{d\lambda_{1}d\lambda_{2}d\Lambda_{1}d\Lambda_{2}}{(2\pi)^{4}}\langle
e^{i(\lambda_{1}\delta_{1}+\Lambda_{1}\delta_{R,1}+\lambda_{2}\delta_{2}+\Lambda_{2}\delta_{R,2})}\rangle\tilde{F}_{\text{\bf
x}}(\lambda_{1})\tilde{F}_{\text{\bf
x}}(\lambda_{2})\tilde{m}(\Lambda_{1})\tilde{m}(\Lambda_{2})$
$\displaystyle=\int\frac{d\lambda_{1}d\lambda_{2}d\Lambda_{1}d\Lambda_{2}}{(2\pi)^{4}}\tilde{F}_{\text{\bf
x}}(\lambda_{1})\tilde{F}_{\text{\bf
x}}(\lambda_{2})\tilde{m}(\Lambda_{1})\tilde{m}(\Lambda_{2})e^{-\frac{1}{2}(\lambda_{1}^{2}+\lambda_{2}^{2})\sigma^{2}-\frac{1}{2}(\Lambda_{1}^{2}+\Lambda_{2}^{2})\sigma^{2}_{RR}}$
$\displaystyle\quad\times\Big{[}1-(\lambda_{1}\Lambda_{1}+\lambda_{2}\Lambda_{2})\sigma^{2}_{R}-\lambda_{1}\lambda_{2}\xi(r)-\Lambda_{1}\Lambda_{2}\xi_{RR}(r)-(\lambda_{1}\Lambda_{2}+\lambda_{2}\Lambda_{1})\xi_{R}(r)+\cdots\Big{]}$
$\displaystyle=\bar{m}^{2}\Big{[}1+(b_{1}^{E})^{2}\xi(r)+B_{1}^{2}\xi_{RR}(r)+2b_{1}^{E}B_{1}\xi_{R}(r)+\cdots\Big{]},$
(22)

Figure 3: Smoothed correlation functions defined in defined in eqs. (14)-(16)
and appearing in the Eulerian marked correlation function [Eq. (24)]. We use a
Gaussian damping kernel $W_{R}$ of width $R=10\,h^{-1}\text{Mpc}$ to smooth
the density fields.

where the ellipsis denotes second order terms in $\xi_{{},R,RR}(r)$. In the
second equality, we have shifted to Fourier space ($\delta\rightarrow\lambda$,
$\delta_{R}\rightarrow\Lambda$) and in the third equality we used the cumulant
expansion theorem,

$\langle e^{iX}\rangle=\exp\left[\sum_{N=1}^{\infty}\frac{i^{N}}{N!}\langle
X^{N}\rangle_{c}\right],$ (23)

and expanded out of the exponential all terms but those containing
$\sigma^{2}$ and $\sigma^{2}_{RR}$, such that we can use Eqs. (17) and (19) to
get biases from spectral parameters, as we did in the last equality. The mCF
becomes

$\mathcal{M}^{E}(r)=\frac{1+W(r)}{1+\xi_{X}(r)}=\frac{1+(b_{1}^{E})^{2}\xi(r)+B_{1}^{2}\xi_{RR}(r)+2b_{1}^{E}B_{1}\xi_{R}(r)+\cdots}{1+(b_{1}^{E})^{2}\xi(r)+\cdots}.$
(24)

In figure 3 we show the correlation functions defined in eqs. (14)-(16) and
used to construct the above mCF. We use a Gaussian smoothing kernel $W_{R}$
with width $R=10\,h^{-1}\text{Mpc}$.

We notice that the zero-lag correlators $\sigma^{2}$ and $\sigma^{2}_{RR}$ do
not appear in the mCF as is guaranteed because we are using renormalized $b$
and $B$ parameters. Meanwhile, the cross-covariance $\sigma^{2}_{R}$ is
canceled out by the mean mark squared appearing in the definition of the mCF
in eq. (5). We remark that a process of renormalization of the Taylor
expansion coefficients $C$ is not strictly necessary because the scale $R$ is
physical, and chosen from the beginning by the observer to mark the tracers,
so there is no impediment for $\sigma_{RR}^{2}$ to appear in the final
expressions. However, the use of parameters $B$ instead of the Taylor
coefficients leads to simpler final equations. Furthermore, in real
applications, we have to use a mark computed by the number density of
galaxies. Another advantage of using the renormalized $B$ parameters is that
the effect of this reassignment can be included simply by re-scaling the $B$
parameters. For the application to dark matter halos and galaxies, one should
treat them as free parameters and fit them with simulations.

##  5 Lagrangian space

In Lagrangian space one considers regions of space, at an initial early time
with spatial coordinates q. We relate matter and tracers overdensities, in an
analogous way to eq. (4), by

$1+\delta_{X}({\text{\bf q}})=F[\delta({\text{\bf q}})],$ (25)

with the difference that the overdensity fields are taken at an initial early
time.333This initial time is taken early enough that all scales of interest
are still in the linear regime and one can safely take $\delta({\text{\bf
q}})\ll 1$, but late enough that the background universe is already in the
Einstein-de Sitter phase. Assuming number conservation of tracers,
$(1+\delta_{X}({\text{\bf x}}))d^{3}x=(1+\delta_{X}({\text{\bf q}}))d^{3}q$,
one can relate the function $F_{\text{\bf x}}$, introduced in eq. (4), and $F$
through

$F_{\text{\bf x}}[\delta({\text{\bf x}})]=\int\frac{d^{3}k}{(2\pi)^{3}}\int
d^{3}q\,e^{i{\text{\bf k}}\cdot({\text{\bf x}}-{\text{\bf
q}})}\int\frac{d\lambda}{2\pi}\tilde{F}(\lambda)e^{i\lambda\delta({\text{\bf
q}})-i{\text{\bf k}}\cdot\mathbf{\Psi({\text{\bf q}},t)}},$ (26)

with $\mathbf{\Psi}({\text{\bf q}})$ the Lagrangian displacement vector, that
maps Lagrangian coordinates q to Eulerian coordinates x as ${\text{\bf
x}}({\text{\bf q}},t)={\text{\bf q}}+\mathbf{\Psi}({\text{\bf q}},t)$.
Equivalently to eq. (17), we introduce the renormalized Lagrangian local
biases [35] with

$b_{n}=\int\frac{d\lambda}{2\pi}e^{-\lambda^{2}\sigma^{2}/2}\tilde{F}(\lambda)(i\lambda)^{n}.$
(27)

The next step is to evolve initially biased tracers with overdensity
$\delta_{X}({\text{\bf q}})$ using Convolution Lagrangian perturbation
theory[38, 39] and thereafter assign them a mark $m[\delta_{R}({\text{\bf
x}})]$ at the moment of observation. The use of the Lagrangian approach has
some advantages with respect to the Eulerian. In the first place it is well
known that the two-point correlation function is poorly modeled within the
Eulerian approach, particularly at the BAO peak position; second, the
(renormalized) Lagrangian bias parameters are obtained through the peak
background split prescription [40, 41, 42, 43], and hence are physically
appealing. But the price to pay is that the mCF is cumbersome to compute,
since it becomes a double Gaussian convolution[17]:

$\displaystyle 1+W(r)$
$\displaystyle=\int\frac{d^{3}q\,e^{-\frac{1}{2}({\text{\bf r}}-{\text{\bf
q}})^{T}\mathbf{A}{-1}({\text{\bf r}}-{\text{\bf
q}})}}{(2\pi)^{3/2}|\mathbf{A}|^{1/2}}\int\frac{d^{3}Q\,e^{-\frac{1}{2}({\text{\bf
R}}-{\text{\bf Q}})^{T}\mathbf{C}^{-1}({\text{\bf R}}-{\text{\bf
Q}})}}{(2\pi)^{3/2}|\mathbf{C}|^{1/2}}$ (28)
$\displaystyle\quad\qquad\times\Big{(}1+\mathcal{I}^{0}({\text{\bf
r}},{\text{\bf q}},{\text{\bf R}},{\text{\bf Q}};b_{i},B_{i})\Big{)},$ (29)

with $\mathcal{I}^{0}$ a function, whose expression is quite large and not
very illuminating, and contains the information of the mark (through the
expansion parameters $B_{i}$) and of the tracers (through the bias parameters
$b_{i}$). It can be found in eq. (B.8) of Ref. 17. The matrices A and C are
given by

$A_{ij}({\text{\bf q}})=2\int\frac{d^{3}p}{(2\pi)^{3}}\big{(}1-e^{i{\text{\bf
p}}\cdot{\text{\bf q}}}\big{)}\frac{p_{i}p_{j}}{p^{4}}P_{L}(p),$ (30)

and

$C_{ij}({\text{\bf q}})=-\frac{1}{4}A_{ij}({\text{\bf
q}})+\frac{1}{6\pi^{2}}\int_{0}^{\infty}P_{L}(p)dp\,\,\delta_{ij},$ (31)

with $P_{L}(p)$ the linear matter power spectrum. We notice that if
$1+\mathcal{I}$ is not a function of Q, it can be pulled out of the Q integral
and due to that $C_{ij}$ depends only on q, the integration over Q gives 1.
This is the case of the “standard” correlation function in Convolution
Lagrangian Perturbation Theory (CLPT),[38] reducing the double Gaussian
convolution in eq. (28), to a single three-dimensional convolution.

Figure 4: Gaussian kernels in eq. (28). GK-$q$ is the kernel of the
$q$-integral (black curves) and GK-$Q$ the kernel of the $Q$-integral (red
curves). Dashed and solid curves show the regions that enclose the 68% and 95%
of the volume respectively. A similar plot can be found in Ref. 44. Figure
adapted from Ref. 17.

In figure 4 we show level plots for the Gaussian kernels appearing in eq.
(28). GK-$q$ is the kernel of the $q$-integral and GK-$Q$ the kernel of the
$Q$-integral. These are shown as a function of $r_{\parallel}$
($R_{\parallel}$) and $r_{\perp}$ ($R_{\perp}$), the components of r (R)
parallel and perpendicular to the Lagrangian coordinate q (Q) with
$Q=q=100\,h^{-1}\,\text{Mpc}$ fixed. To have a sense of the meaning of this
plot, consider that the kernel GK-$q$ is the probability distribution of
finding two dark matter particles separated by a distance r at the redshift of
evaluation ($z=0.5$ in the plot), given that they were separated by a distance
q at the early initial time.[44] We notice that both kernels have their
maximum value at $Q=q=100\,h^{-1}\,\text{Mpc}$, but the GK-$Q$ is more sharply
peaked because $|\mathbf{C}|<|\mathbf{A}_{L}|$. This observation suggests to
approximate

$\frac{e^{-\frac{1}{2}({\text{\bf R}}-{\text{\bf
Q}})^{T}\mathbf{C}^{-1}({\text{\bf R}}-{\text{\bf
Q}})}}{(2\pi)^{3/2}|\mathbf{C}|^{1/2}}\approx\delta_{\text{D}}({\text{\bf
R}}-{\text{\bf Q}}).$ (32)

By doing so, one arrives at

$1+W^{\text{W16}}(r)=\int\frac{d^{3}q}{(2\pi)^{3/2}|A_{L}|^{1/2}}e^{-\frac{1}{2}({\text{\bf
r}}-{\text{\bf q}})^{T}\mathbf{A}^{-1}_{L}({\text{\bf r}}-{\text{\bf
q}})}\Big{(}1+J^{0}({\text{\bf q}},{\text{\bf r}};b_{i},B_{i})\Big{)}$ (33)

which is the generalization of the Zeldovich approximation (ZA) result
obtained in Ref. 16. The expression for $J^{0}({\text{\bf q}},{\text{\bf
r}};b_{i},B_{i})$ is given by eq. (4.18) of Ref. 17.

Figure 5: Marked correlation function modeling comparison. The left panel
shows the mCF with parameters $b_{1}=1.2$, $b_{2}=0.2$, $B_{1}=-0.7$,
$B_{2}=B_{1}^{2}$, for the different methods: solid red curve is the hybrid
model presented in Ref. 17; dashed blue the method of Eq. (33), which is the
method introduced in Ref. 16, plus 1-loop nonlinear contributions; in dotted
gray, the ZA; and in dot-dashed green, the Eulerian linear model of eq. (24).
The right panel shows the relative differences with respect to the ZA. The
gray dot-dashed horizontal line denotes the $1\%$ differences. The solid black
curve in the right panel shows the relative differences between the CLPT
1-loop and ZA standard correlation functions. Figure adapted from Ref. 17.

Figure 5 shows the mCF using different analytical methods for the MG Hu-
Sawicki F5 model: solid red curve is the one presented above and developed in
Ref. 17—called hybrid from now on, because it evolves the tracers using LPT
but marks the tracers using the Eulerian smoothed matter overdensity fields at
the moment of observation. The dashed blue is the method of eq. (33), which is
that introduced in Ref. 16 with the addition of leading non-linear (loop)
corrections; in dotted gray, the Zeldovich Approximation, which in this work
means the linear model of Ref. 16; and, in dotdashed green, the Eulerian
linear model of eq. (24). We are using Lagrangian local biases $b_{1}=1.2$ and
$b_{2}=0.2$, and mark parameters $B_{1}=-0.7$, $B_{2}=B_{1}^{2}$. In the right
panel we show the relative difference between the models and the ZA. The
agreement between them is very good, even for the linear Eulerian theory,
below 1% for scales above $\sim 20\,h^{-1}\,\text{Mpc}$. The relative
differences between the CLPT standard correlation function with and without
loop contributions ($\xi_{\text{CLPT}}$ and $\xi_{\text{ZA}}$, respectively)
is shown in the right panel (black curve); by comparing it with the dashed
blue curve (which shows the relative differences of the W16+1-loop and ZA
mCFs) we confirm that mCFs which efficiently enhance low density regions in
the sky are indeed more linear than the standard correlation function.

##  6 Comparison to simulations

We use the Elephant suite of simulations run with the ECOSMOG code. [14] This
suite contains five realisations of the initial conditions and for each
realisation we have one simulation of $\Lambda$CDM (GR) together with three
simulations of the Hu-Sawicki $f(R)$ model with parameters $f_{R0}=-10^{-4}$
(F4), $f_{R0}=-10^{-5}$ (F5) and $f_{R0}=-10^{-6}$ (F6). It also contains
galaxy mock catalogs that were made with a Halo Occupation Distribution (HOD)
method. The HOD parameters for the $\Lambda$CDM model are the best-fit
parameter values from the CMASS data.[45] For the $f(R)$ models, we tune the
HOD parameters so that we reproduce the correlation function in the
$\Lambda$CDM model. The simulations were ran in a box of size
$L=1024\,h^{-1}\,\text{Mpc}$ with $N=1024^{3}$ particles and the cosmological
parameters used to make the initial conditions were $\Omega_{b}=0.046$,
$\Omega_{\Lambda}=0.719$, $\Omega_{m}=0.281$, $h=0.697$, $\sigma_{8}=0.82$ and
$n_{s}=0.971$, and for our analysis we choose snapshots at redshift $z=0.5$.
To compute the mark for each of our tracers (dark matter particles, halos or
mock galaxies) we binned the particles/halos/mock galaxies to a grid with
gridsize of $20\,h^{-1}\,\text{Mpc}$ (corresponding to a $N=52^{3}$ grid)
using a Nearest Grid Point assignment scheme to get an estimate for the
density for which the mark depends on.

Both, the standard and weighted 2-point correlation functions were computed
using the Correlation Utilities and Two-point Estimates (CUTE)[46]
code.444https://github.com/damonge/CUTE From this the marked correlation
function follows simply as $\mathcal{M}(r)=(1+W(r))/(1+\xi(r))$.

We consider the White-mark with $\rho_{*}=10$, $p=7$, corresponding to
coefficients $C_{0}=1$, $C_{1}=-0.64$, $C_{2}=0.46$ in the Taylor expansion of
the mark function in eq. (10). In our analytical models we smooth the matter
fields that assign the mark with a top-hat filter $W_{R}$ of radius
$R=10\,h^{-1}\,\text{Mpc}$.

Figure 6: Ratios of different mCFs analytical models to the ZA mCF for the F5
gravity (we obtain similar results for F4, F6 and GR). We show the W16 model
plus 1-loop corrections (dashed blue); the hybrid model (solid red); and the
linear Eulerian model of eq. (24) (dot-dashed green). Figure adapted from Ref.
17.

We first confront different analytical methods against matter data, since in
this case the mCF does not suffer from biasing-marking degeneracies and we can
observe more neatly the effects of the weights. The expansion parameters $B$
are obtained from eq. (20): $B_{1}^{\text{GR}}=-0.6495$,
$B_{1}^{\text{F6}}=-0.6497$, $B_{1}^{\text{F5}}=-0.6505$,
$B_{1}^{\text{F4}}=-0.6522$, $B_{2}^{\text{GR}}=0.4628$,
$B_{2}^{\text{F6}}=0.4502$, $B_{2}^{\text{F5}}=0.4495$,
$B_{2}^{\text{F4}}=0.4480$; which are numerically very close to $C_{1}$ and
$C_{2}$ values because $C_{0}=1$ and the variances $\sigma^{2}_{RR}$ are
small.

In figure 6 we show the ratios of the different analytical methods to the ZA
model. We do it for the Hu-Sawicki F5 gravity model, and we have checked also
the gravity theories GR, F4 and F6, and obtained similar results. The dashed
blue curve shows the analytical result with the W16 model plus 1-loop
corrections; solid red, the hybrid method; and dot-dashed green, the linear
Eulerian model of eq. (24). The differences between the analytical methods are
apparent but small, being lesser than the $1\%$. At scales
$r>40\,h^{-1}\,\text{Mpc}$, all analytic models are indistinguishable and they
are within the errors of the simulation data. At smaller scales, the method
hybrid outperforms the other perturbative approaches and captures reasonably
the trend of the data all the way up to the smoothing scale
$R=10\,h^{-1}\,\text{Mpc}$.

Figure 7: Halos and galaxies mCFs. Left panel is for halo masses
$12.5<\log_{10}[M/(M_{\odot}h^{-1})]<15$ and gravitational models F4, F5, F6
and GR. Right panel is for HOD and shows only the F4 and GR models. The linear
local biases are fitted by using the ZA correlation function. Figure adapted
from Ref. 17.

Finally, in figure 7 we use the hybrid model to fit the data for halos and
galaxies (HOD), finding very good agreement between simulations and theory for
$r>30\,h^{-1}\,\text{Mpc}$. In the left panel we show that the MG models are
differentiated by the data and the analytical curves in the interval
$30<r<50\,h^{-1}\,\text{Mpc}$, with the exception of the F5 and F6 models,
which become indistinguishable within the error bars. However, we notice that
the HOD mCFs for the different MG models do not look differently enough to be
distinguished by our data. We think that the use of different parameters
$\delta_{*}$ and $p$ or different mark functions can help to discern between
different gravitational patterns; see, for example, Ref. 18.

##  7 Marked power spectrum

Recently, it was discovered that the marked statistics can be very useful for
the estimation of cosmological parameters, particularly the mass scale of the
primordial neutrinos.[47] The reason is that given the large velocity
dispersion of neutrinos, they do not contribute to the small scale clustering
(see, e.g. Ref. 48). But instead, they free-stream over large distances,
allowing them to populate the cosmic voids considerably. On the other hand,
the influence of neutrinos on high density regions, for example halos, is
quite marginal because their contribution to the total matter density is very
small. Thus, statistics that up-weight low density regions can be very
beneficial for parameter estimation. In fact, the Fisher formalism results are
quite impressive: the marked power spectrum of the matter field outperforms
the constraining power of the standard power spectrum by at least a factor of
2 on all cosmological parameters, and in particular of up to 80 times for the
sum of the neutrino masses, when a combination of two different mark functions
is used.

The analytical description of the marked power spectrum has been studied in
detail in a couple of recent works. First, for matter in real space[49] and
later for tracers in redshift space.[19] The main complication of such studies
is that the mark parameters optimal for improving the cosmological parameter
estimation are not able for a perturbative treatment, and the linear theory
fails at all scales. The large-scale theory contains non-negligible
contributions from all perturbative orders; in Ref. 19 a reorganization of the
theory that contains all terms relevant on large-scales is proposed, with an
explicit form at one-loop and the following structure at infinite-loop:

$\displaystyle\mathcal{M}(k,\mu)$
$\displaystyle=\big{[}C_{0}-C_{1}W_{R}(k)\big{]}\big{[}(\tilde{a}_{0}+\tilde{a}_{2}\mu^{2})P_{L}(k)+\tilde{b}_{0}\big{]}+\text{EFT}+\text{shot
noise}$

with $\tilde{a}_{0}$, $\tilde{a}_{2}$ and $\tilde{b}_{0}$ free parameters of
the theory. EFT refers to the standard 1-loop corrections to the linear power
spectrum and the effective field theory terms, while the shot noise includes
the stochastic contributions.

##  8 Summary

We have reviewed the analytical study of marked statistics that up-weight low
density regions in the Universe.[16, 17, 19, 18] The idea behind marked
statistics is to assign a value (the mark) to each entity in a catalogue and
perform statistics over the resulting weighted objects. An efficient way to
suppress non-linearities is by choosing a mark that gives more weight to
objects that reside in low density regions, which become less important as one
consider higher density regions. This marking process can be particularly
useful for testing general relativity with cosmological probes, since modified
gravity models that provide cosmic acceleration, often rely on screening
mechanisms to hide its range and strength on high density environments, and
the screenings switch off in regions that are more depleted of matter, like
cosmic voids. Therefore, the effects of an hypothetical MG are expected to be
more pronounced in low density regions. It is then natural to search for such
marked statistics, as proposed in Ref. 16. Further, these statistics can
potentially improve the constraining power of all cosmological parameters, in
particular the mass of the primordial neutrinos, as it was shown in Ref. 47.

## Acknowledgments

I want to thank Alfredo Macias and Carlos Herdeiro for inviting me to give a
talk in the BS1 parallel session “Scalar fields in Cosmology” of the Sixteenth
Marcel Grossmann Meeting – MG16, July 5-10 2021. I also want to thank Jorge L.
Cervantes-Cota, Baojiu Li, Kazuya Koyama, Hans Winther, Elisa Massara, Oliver
Philcox and Martin White for their collaboration in developing the material
that led to this work. I acknowledge support by CONACyT project 283151 and
CONACyT Ciencia de Frontera grant No. 102958. I would like to thank the DiRAC
Data Centric system at Durham University (www.dirac.ac.uk) for computational
facilities.

Hu-Sawicki MG model

The Hu-Sawicki model [11] is a particular realisation of $f(R)$ gravity that
is able to evade the strong constraints coming from local test of gravity and
still give rise to interesting observable signatures on cosmological scales.
$f(R)$ theories consist on replace the Einstein-Hilbert Lagrangian density
$\sqrt{-g}R$ by a general function of the Ricci scalar
$\sqrt{-g}\big{(}R+f(R)\big{)}$; see Ref. 50 for a review. The Hu-Sawicki
model is defined by the function

$f(R)=-M^{2}\frac{c_{1}(R/M^{2})^{n}}{c_{2}(R/M^{2})^{n}+1}$ (34)

where the energy scale is chosen to be $M^{2}=\Omega_{m}H^{2}_{0}$. In this
parametrized model, at high curvature ($R\ll M^{2}$) the function $f(R)$
approaches a constant, recovering GR with cosmological constant, while at low
curvature it goes to zero, recovering GR; the manner these two behaviors are
interpolated is dictated by the free parameters. Given that
$d^{2}f(R)/dR^{2}>0$ for $R>M^{2}$, the solutions are stable and the scalar
tensor gravity description is possible.[12] In order to be have a similar
background evolution than the $\Lambda$CDM model it is also necessary that
$c_{1}/c_{2}=6\Omega_{\Lambda}/\Omega_{m}$, thus leaving two parameters to fix
the model. One choose these parameters to be $n$ and $f_{R0}\equiv
df(R)/dR|_{R=R_{0}}$, with $R_{0}=3H_{0}^{2}(\Omega_{m}+4\Omega_{\Lambda})$
the Ricci scalar of the background cosmological metric evaluated nowadays.

In this work we have considered a fixed value $n=1$, and hence the model has
only one free parameter $f_{R0}$. Taking this parameter to zero we recover GR.
The three choices for the parameters we are considering in this paper (F4, F5
and F6) are such that they lie in the region around where the best constraints
lie today. The F4 model corresponds to $|f_{R0}|=10^{-4}$ (is in tension with
local experiments), the F5 model to $|f_{R0}|=10^{-5}$ (agrees with most
experiments and observations, but are in tension with others) and the F6 model
to $|f_{R0}|=10^{-6}$ (which is still allowed). All these models are
completely compatible with observations of the expansion history of the
Universe, but present large deviations from the $\Lambda$CDM in the clustering
and formation of structures.

## References

  * [1] N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters (2018).
  * [2] J. Khoury and A. Weltman, Chameleon fields: Awaiting surprises for tests of gravity in space, Phys. Rev. Lett. 93, p. 171104 (2004).
  * [3] J. Khoury and A. Weltman, Chameleon cosmology, Phys. Rev. D69, p. 044026 (2004).
  * [4] L. Verde, T. Treu and A. G. Riess, Tensions between the Early and the Late Universe, Nature Astron. 3, p. 891 (7 2019).
  * [5] E. Bertschinger, On the Growth of Perturbations as a Test of Dark Energy, Astrophys. J. 648, 797 (2006).
  * [6] E. Bertschinger and P. Zukin, Distinguishing Modified Gravity from Dark Energy, Phys. Rev. D78, p. 024015 (2008).
  * [7] E. V. Linder, No Slip Gravity, JCAP 03, p. 005 (2018).
  * [8] A. A. Starobinsky, A New Type of Isotropic Cosmological Models Without Singularity, Phys. Lett. B 91, 99 (1980).
  * [9] S. M. Carroll, V. Duvvuri, M. Trodden and M. S. Turner, Is cosmic speed - up due to new gravitational physics?, Phys. Rev. D70, p. 043528 (2004).
  * [10] S. Capozziello, Curvature quintessence, Int. J. Mod. Phys. D11, 483 (2002).
  * [11] W. Hu and I. Sawicki, Models of f(R) Cosmic Acceleration that Evade Solar-System Tests, Phys. Rev. D76, p. 064004 (2007).
  * [12] G. Magnano and L. M. Sokolowski, On physical equivalence between nonlinear gravity theories and a general relativistic selfgravitating scalar field, Phys. Rev. D 50, 5039 (1994).
  * [13] L. G. Jaime, L. Patino and M. Salgado, Robust approach to f(R) gravity, Phys. Rev. D 83, p. 024039 (2011).
  * [14] B. Li, G.-B. Zhao, R. Teyssier and K. Koyama, ECOSMOG: an Efficient COde for Simulating MOdified Gravity, JCAP 2012, p. 051 (Jan 2012).
  * [15] M. Cautun, E. Paillas, Y.-C. Cai, S. Bose, J. Armijo, B. Li and N. Padilla, The Santiago–Harvard–Edinburgh–Durham void comparison – I. SHEDding light on chameleon gravity tests, Mon. Not. Roy. Astron. Soc. 476, 3195 (2018).
  * [16] M. White, A marked correlation function for constraining modified gravity models, JCAP 1611, p. 057 (2016).
  * [17] A. Aviles, K. Koyama, J. L. Cervantes-Cota, H. A. Winther and B. Li, Marked correlation functions in perturbation theory, JCAP 01, p. 006 (2020).
  * [18] S. Alam, C. Arnold, A. Aviles et al., Testing the theory of gravity with DESI: estimators, predictions and simulation requirements (11 2020).
  * [19] O. H. E. Philcox, A. Aviles and E. Massara, Modeling the Marked Spectrum of Matter and Biased Tracers in Real- and Redshift-Space, JCAP 03, p. 038 (2021).
  * [20] A. Aviles and J. L. Cervantes-Cota, Lagrangian perturbation theory for modified gravity, Phys. Rev. D96, p. 123526 (2017).
  * [21] A. Aviles, G. Valogiannis, M. A. Rodriguez-Meza, J. L. Cervantes-Cota, B. Li and R. Bean, Redshift space power spectrum beyond Einstein-de Sitter kernels, JCAP 04, p. 039 (2021).
  * [22] C. Beisbart and M. Kerscher, Luminosity- and morphology-dependent clustering of galaxies, Astrophys. J. 545, p. 6 (2000).
  * [23] R. K. Sheth, The halo-model description of marked statistics, Mon. Not. Roy. Astron. Soc. 364, p. 796 (2005).
  * [24] R. K. Sheth, A. J. Connolly and R. Skibba, Marked correlations in galaxy formation models, Submitted to: Mon. Not. Roy. Astron. Soc. (2005).
  * [25] M. White and N. Padmanabhan, Breaking Halo Occupation Degeneracies with Marked Statistics, Mon. Not. Roy. Astron. Soc. 395, p. 2381 (2009).
  * [26] G. Valogiannis and R. Bean, Beyond $\delta$: Tailoring marked statistics to reveal modified gravity, Phys. Rev. D97, p. 023535 (2018).
  * [27] C. Hernández-Aguayo, C. M. Baugh and B. Li, Marked clustering statistics in $f(R)$ gravity cosmologies, Mon. Not. Roy. Astron. Soc. 479, 4824 (2018).
  * [28] J. Armijo, Y.-C. Cai, N. Padilla, B. Li and J. A. Peacock, Testing modified gravity using a marked correlation function, Mon. Not. Roy. Astron. Soc. 478, 3627 (2018).
  * [29] V. Desjacques, D. Jeong and F. Schmidt, Large-Scale Galaxy Bias, Phys. Rept. 733, 1 (2018).
  * [30] A. Aviles, M. A. Rodriguez-Meza, J. De-Santiago and J. L. Cervantes-Cota, Nonlinear evolution of initially biased tracers in modified gravity, JCAP 1811, p. 013 (2018).
  * [31] S. Gottloeber, M. Kerscher, A. V. Kravtsov, A. Faltenbacher, A. Klypin and V. Mueller, Spatial distribution of galactic halos and their merger histories, Astron. Astrophys. 387, p. 778 (2002).
  * [32] R. Skibba, R. K. Sheth, A. J. Connolly and R. Scranton, The luminosity-weighted or ‘marked’ correlation function, Mon. Not. Roy. Astron. Soc. 369, 68 (2006).
  * [33] M. White, The Zel’dovich approximation, Mon. Not. Roy. Astron. Soc. 439, 3630 (2014).
  * [34] T. Baldauf, M. Mirbabayi, M. Simonović and M. Zaldarriaga, Equivalence Principle and the Baryon Acoustic Peak, Phys. Rev. D92, p. 043514 (2015).
  * [35] T. Matsubara, Nonlinear perturbation theory with halo bias and redshift-space distortions via the Lagrangian picture, Phys. Rev. D78, p. 083519 (2008), [Erratum: Phys. Rev.D78,109901(2008)].
  * [36] A. Aviles, Renormalization of Lagrangian bias via spectral parameters, Phys. Rev. D98, p. 083541 (2018).
  * [37] A. Eggemeier, R. Scoccimarro and R. E. Smith, Bias Loop Corrections to the Galaxy Bispectrum, Phys. Rev. D99, p. 123514 (2019).
  * [38] J. Carlson, B. Reid and M. White, Convolution Lagrangian perturbation theory for biased tracers, Mon. Not. Roy. Astron. Soc. 429, p. 1674 (2013).
  * [39] Z. Vlah, M. White and A. Aviles, A Lagrangian effective field theory, JCAP 1509, p. 014 (2015).
  * [40] N. Kaiser, On the Spatial correlations of Abell clusters, Astrophys. J. 284, L9 (1984).
  * [41] H. J. Mo, Y. P. Jing and S. D. M. White, High-order correlations of peaks and halos: A Step toward understanding galaxy biasing, Mon. Not. Roy. Astron. Soc. 284, p. 189 (1997).
  * [42] R. K. Sheth and G. Tormen, Large scale bias and the peak background split, Mon. Not. Roy. Astron. Soc. 308, p. 119 (1999).
  * [43] F. Schmidt, D. Jeong and V. Desjacques, Peak-Background Split, Renormalization, and Galaxy Clustering, Phys. Rev. D88, p. 023515 (2013).
  * [44] S. Tassev, Lagrangian or Eulerian; Real or Fourier? Not All Approaches to Large-Scale Structure Are Created Equal, JCAP 1406, p. 008 (2014).
  * [45] M. Manera, R. Scoccimarro, W. J. Percival, L. Samushia, C. K. McBride, A. J. Ross, R. K. Sheth, M. White, B. A. Reid, A. G. Sánchez, R. de Putter, X. Xu, A. A. Berlind, J. Brinkmann, C. Maraston, B. Nichol, F. Montesano, N. Padmanabhan, R. A. Skibba, R. Tojeiro and B. A. Weaver, The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: a large sample of mock galaxy catalogues, Mon. Not. Roy. Astron. Soc. 428, 1036 (Jan 2013).
  * [46] D. Alonso, CUTE solutions for two-point correlation functions from large cosmological datasets, arXiv e-prints , p. arXiv:1210.1833 (Oct 2012).
  * [47] E. Massara, F. Villaescusa-Navarro, S. Ho, N. Dalal and D. N. Spergel, Using the Marked Power Spectrum to Detect the Signature of Neutrinos in Large-Scale Structure, Phys. Rev. Lett. 126, p. 011301 (2021).
  * [48] J. Lesgourgues and S. Pastor, Massive neutrinos and cosmology, Phys. Rept. 429, 307 (2006).
  * [49] O. H. E. Philcox, E. Massara and D. N. Spergel, What does the marked power spectrum measure? Insights from perturbation theory, Phys. Rev. D 102, p. 043516 (2020).
  * [50] K. Koyama, Cosmological Tests of Modified Gravity, Rept. Prog. Phys. 79, p. 046902 (2016).


# Factor Graphs for Heterogeneous Bayesian Decentralized Data Fusion

Ofer Dagan  Smead Aerospace Engineering Sciences Dept.
University of Colorado Boulder
Boulder, USA
<EMAIL_ADDRESS>Nisar R. Ahmed  Smead Aerospace Engineering Sciences
Dept.
University of Colorado Boulder
Boulder, USA
<EMAIL_ADDRESS>
###### Abstract

This paper explores the use of factor graphs as an inference and analysis tool
for Bayesian peer-to-peer decentralized data fusion. We propose a framework by
which agents can each use local factor graphs to represent relevant partitions
of a complex global joint probability distribution, thus allowing them to
avoid reasoning over the entirety of a more complex model and saving
communication as well as computation cost. This allows heterogeneous multi-
robot systems to cooperate on a variety of real world, task oriented missions,
where scalability and modularity are key. To develop the initial theory and
analyze the limits of this approach, we focus our attention on static linear
Gaussian systems in tree-structured networks and use Channel Filters (also
represented by factor graphs) to explicitly track common information. We
discuss how this representation can be used to describe various multi-robot
applications and to design and analyze new heterogeneous data fusion
algorithms. We validate our method in simulations of a multi-agent multi-
target tracking and cooperative multi-agent mapping problems, and discuss the
computation and communication gains of this approach.

###### Index Terms:

Bayesian decentralized data fusion (DDF), factor graphs, heterogeneous multi-
robot systems, sensor fusion.

##  I Introduction

In the context of multi-agent systems, there is a growing interest to allow
for the coordinated and cooperative operation of _heterogeneous systems_ in
uncertain environments. For example a team of air and ground vehicles can vary
in their sensing/computing/communication capabilities and/or might have a
different set of models and objectives. The question is how to share
information in a scalable and modular way, despite this heterogeneity?

If we assume that the underlying task is to locally infer some quantity
(state) of the system or environment, such as the temperature field or an
environment map and individual robot positions in Fig. 1(a)-(b), respectively,
then the problem translates to a data fusion problem. While there different
methods for data fusion, this paper focuses on Bayesian decentralized data
fusion (DDF). More specifically, this paper builds the foundations for DDF on
factor graphs, to enable large scale multi-agent heterogeneous systems by
reducing local communication and computation costs. The main idea is to allow
each agent to only reason about its local task and communicate only the
information relevant to its neighbor.

For example, when cooperatively estimating the local biases and temperature
field in a lab [1] (Fig. 1(a)), it is clear that effect of far away sensors is
negligible, while the influence of neighborhood sensors can’t be ignored. Thus
allowing each agent to process only locally relevant states instead of the
full field might be crucial for scalability. Similarly, in the smoothing and
mapping (SAM) application in [2] (Fig. 1(b)), allowing each agent to only
reason about relevant parts of the map and its own position scales with the
size of the map and network. But in both of these examples, for a correct
estimate, we must properly account for the hidden correlations induced by the
local states.

Figure 1: a) Application examples for cooperative heterogeneous tasks (figures
adapted with permission from [3], [2]). b) A full system inference problem
represented by a factor graph and its subdivision into agent’s local graphs.
Arrows show the tree-structured communication graph.

Proper treatment of the problem requires revisiting DDF theory, as it builds
on the assumptions that all agents have the same monolithic/homogeneous model.
Further, identifying the conditional independence structure of the problem, to
enable heterogeneous fusion, might not be trivial for other real world
applications and complex dependencies. This raises the question, _What is an
appropriate representation for real world, mission and task oriented
heterogeneous multi-robot systems, to allow for the design and analysis of
decentralized data fusion algorithms?_

Probabilistic graphical models (PGMs) are a powerful and useful tool for
inference problems [4], where it is assumed that the full system graph,
representing the joint distribution over all variables, is known and can be
reasoned about. However, in decentralized systems this requires each agent to
hold a copy of the full graph, irrespective to its local task or
variables/states of interest. Instead, as shown in Fig. 1(c), we suggest
“slicing” the graph into local sub-graphs, such that each agent only reasons
about the part of the graph relevant to its task.

The goal of this paper is to develop the theory and analyze the limitations of
this approach using factor graphs. We show that factor graphs are a natural
choice for this aim, as they can represent a variety of robotics applications
with the focus on heterogeneous DDF problems. We present _FG-DDF_ \- a
formulation of the DDF problem in terms of factor graphs. We set the
foundations by formulating the problem and defining graph operations in
general terms when possible. We show that the fusion operation becomes almost
trivial, as messages between agents can be treated as “factors” simply added
to the local graph at the receiving agent. Further we show how common
information can be identified (when it can be defined), especially in
heterogeneous fusion problems. To explicitly demonstrate the algorithms
presented, we focus on static linear Gaussian problems in tree-structured
networks and leave conservative marginalization of past states in dynamic
problems and ad-hoc network topologies to future work. Finally, we present
simulation results which show significant computational ($95\%$) and
communications ($97\%$) savings compared to homogeneous fusion.

##  II Problem Statement

Let $V$ be a global set of random variables describing states of interest
monitored by a set $N_{a}$ of $n_{a}$ autonomous agents. The states of
interest are distributed between the agents such that each agent $i$ monitors
some “local states of interest”, which are a subset of the global set
$\chi^{i}\subseteq V$. The set $\chi^{i}$ can be further divided into a set of
local states $\chi^{i}_{L}$, observed locally by agent $i$, and a set of
common states between agent $i$ and its neighbors $\chi^{i}_{C}=\bigcup_{j\in
N_{a}^{i}}\chi^{ij}_{C}$, where $N^{i}_{a}\subseteq N_{a}$ is the set of
agents communicating with agent $i$ and
$\chi^{i}=\chi^{i}_{L}\bigcup\chi^{i}_{C}$.

Each agent is able to gather information about its local states of interest,
$\chi^{i}$, by (i) taking independent local measurements $Y_{k}^{i}$,
described by the measurement likelihood $p(Y^{i}_{k}|\chi^{i})$ and (ii) by
fusing messages from a neighboring agent $j\in\\{1,...,N_{a},i\neq j\\}$
regarding common states $p^{j}(\chi^{ij}_{C}|Z^{j}_{k})$, where
$\chi^{ij}_{C}\subseteq\chi^{i}_{C}$ are the common states to agents $i$ and
$j$ and $Z^{j}_{k}$ is all the information available at agent $j$ up to and
including time step $k$. In the above definitions $\chi$ can represent a
static or dynamic state, and we neglect states time indexing for ease of
notation.

By allowing the pdfs to represent distributions over identical or overlapping
sets of states, along with the flexible definition of the local and common
sets of states, we enable the analysis of a variety of robotics and data
fusion applications, such as static/dynamic, homogeneous/heterogeneous
problems, under one framework. For instance, in the temperature sensing static
example of Fig. 1(a), sensor $i$ bias state would belong to the local set
$\chi^{i}_{L}$ and a neighboring sensor $j$ temperature measurement would
belong to the common set $\chi_{C}^{ij}$. Similarly, in the dynamic SAM
example of Fig. 1(b), agent $i$ position states would be $\chi^{i}_{L}$ and
landmarks communicated from agent $j$ in $\chi^{ij}_{C}$.

The DDF problem is defined as locally inferring the posterior distribution
based on all available information at agents $i$ and $j$ up to and including
time step $k$. Here we look for a fusion function $\mathbb{F}$, that takes as
input the agents prior distributions $p^{i}(\chi^{i}|Z^{i,-}_{k})$ and
$p^{j}(\chi^{j}|Z^{j,-}_{k})$ and returns the locally fused posterior
distributions $p^{i}_{f}(\chi^{i}|Z^{i,+}_{k})$ and
$p^{j}_{f}(\chi^{j}|Z^{j,+}_{k})$:

$\begin{split}(p^{i}_{f}(\chi^{i}|Z^{i,+}_{k}),p^{j}_{f}(\chi^{j}|Z^{j,+}_{k}))=\mathbb{F}(p^{i}(\chi^{i}|Z^{i,-}_{k}),p^{j}(\chi^{j}|Z^{j,-}_{k})),\end{split}$
(1)

where $Z^{j,-}_{k}\equiv Y^{j}_{k}\cup Z^{j}_{k-1}$ is the local information
at agent $j$ before fusion at time step $k$, and $Z^{i,+}_{k}$, $Z^{j,+}_{k}$
is local information at the agents after fusion. Note that in general
$Z^{i,+}_{k}\neq Z^{j,+}_{k}$ as in the case of heterogeneous fusion, where
information regarding local states is not shared.

In this paper we suggest solving the decentralized inference problem with
factor graphs, where each agent holds a factor graph that represents the local
posterior distribution, which is the result of a fusion rule given in (1).
Fusion of new information, whether it is the result of local observations
($Y^{i}_{k}$) or of information from a neighboring agent ($Z_{k}^{j}$) is
performed by adding a factor to the local graph, and estimation of the states
of interest is done by local inference on the graph. We claim that this is a
shift in paradigm for how to analyze and design DDF fusion algorithms, where
instead of maintaining a graph over the full system variables $V$, each agent
holds a smaller graph, only over the subset $\chi^{i}$.

###  II-A Related Work

Algorithms for heterogeneous fusion, as detailed in [5], are advantageous as
they scale with a subset of locally relevant variables instead of the number
of agents in the system. However this is not without cost, as it requires
approximations in order to remove indirect dependencies between variables not
mutually monitored by both agents. The problem becomes acute in filtering
scenarios, where marginalization results in hidden correlations. In [5] a
family of heterogeneous fusion algorithms, conditionally factorized channel
filters (_CF 2_), is developed by utilizing the conditional independence
structure in a multi-agent multi-target tracking application. It is shown that
even in this limited setting of [5] (tree-structured network and full
communication rate), challenges arise in identifying the common information
and in properly accounting for it, as in the case of fixed-lag smoothing or
filtering problems.

In [1] Paskin et al. present a distributed inference problem in a network of
static sensors. It is solved by using a robust message passing algorithm on a
junction tree, but this is limited to static variables and requires the full
construction of the tree before performing inference. In [6] Makarenko et al.
extend Paskin’s algorithm for a dynamic state and formulate it as a DDF
problem, however, the algorithm is limited to a single common state, i.e., the
focus is on a homogeneous problem. Assuming static network topology for a
static model they show equivalence between the channel filter (CF) [7] and the
Hugin [8] algorithms. In [9] Chong and Mori use both Bayesian networks and
information graphs to identify conditional independence and track common
information, respectively. While they use conditional independence to reduce
state dimension, it is only for communication purposes and not for local state
reduction. Further, their graphical model representation is used for analysis
purposes and not for inference on the graph itself.

Thus, PGMs for decentralized data fusion is not yet utilized to the full
extent: for both analyzing problem structure by identifying conditional
independence and for efficient inference on the graph. Here we suggest using
factor graphs to leverage the sparsity induced by the decentralized structure
of the problem.

Factor graphs are used in a variety of robotic applications [10] and are
beneficial in (i) providing insight into the structure of problems across
robotics, (ii) enabling efficient inference algorithms and (iii) analysing the
correctness of local inference. Despite these benefits, to the best of our
knowledge, previous work only looked at limited/niche cases of heterogeneous
fusion [2], or when agents are considered as neighbors only if they share a
relative measurement [11], so present work endeavors for a wider and more
general framework.

##  III Bayesian Decentralized Data Fusion

Given two agents i and j, we are interested in forming an estimate of the
state $\chi$ given all the information available locally, defined by $Z^{i}$
and $Z^{j}$, respectively. Assuming all measurements are conditionally
independent given the state $\chi$ and using a distributed variant of Bayes’
rule, [12] shows that the exact state estimate can be achieved using:

$p_{f}(\chi|Z^{i,-}_{k}\cup
Z^{j,-}_{k})\propto\frac{p^{i}(\chi|Z^{i,-}_{k})p^{j}(\chi|Z^{j,-}_{k})}{p^{ij}_{c}(\chi|Z^{i,-}_{k}\cap
Z^{j,-}_{k})}$ (2)

where $p^{ij}_{c}(\chi|Z^{i,-}_{k}\cap Z^{j,-}_{k})$ is the distribution of
the common information (when it makes sense to define one) of the two agents
which needs to be removed in order to avoid double counting. Notice that the
above fusion rule satisfies (1) for the homogeneous case
$\chi=\chi^{i}=\chi^{j}$.
In this paper we assume that common information is explicitly tracked, by
either keeping a pedigree of all information transformed in fusion instances
or by restricting network topology to a tree structure, as in the channel
filter [7].

For heterogeneous fusion, $\chi^{i}\neq\chi^{j}$, assuming the local states
are conditionally independent given common states:

$\chi^{i}_{L}\perp\chi^{j}_{L}\mid\chi^{ij}_{C},$

Dagan et al. show in [5] that agents can fuse only a subset of common
variables without the need to communicate “irrelevant” local variables. The
fusion rule at agent $i$ can be written as,

$\begin{split}p^{i}_{f}(\chi^{i}|Z^{i,+}_{k})\propto\frac{p^{i}(\chi^{ij}_{C}|Z^{i,-}_{k})p^{j}(\chi^{ij}_{C}|Z^{j,-}_{k})}{p^{ij}_{c}(\chi^{ij}_{C}|Z^{i,-}_{k}\cap
Z^{j,-}_{k})}\cdot p^{i}(\chi^{i}_{L}|\chi^{ij}_{C},Z^{i,-}_{k}).\end{split}$
(3)

In [5], assuming tree topology, this is named the _heterogeneous-state channel
filter_ (HS-CF). As we show later, (3) amounts to sending/incorporating a
factor in the local graph of each of the agents $i/j$.

We move our attention to factor graphs, we describe and define them for
general pdfs and for the special case of Gaussian distributions (Sec. IV). We
then put it all together in Sec. V to formalize factor graphs for DDF.

##  IV Factor Graphs

Figure 2: a) Example factor graph representing a tracking scenario, where the
variables of interest are the target states from time $t=0$ to $t=k$. Shown
are the different types of factors (prior, transition and measurement). b-d)
show the graph operations, namely prediction, roll-up and estimation, with new
factors and variables marked by red empty nodes.

A factor graph is an undirected bipartite graph $\mathcal{F}=(U,V,E)$ with
factor nodes $f_{l}\in U$, random variable nodes $v_{m}\in V$ and edges
$e_{lm}\in E$ connecting factor node _l_ to variable node _m_. Each factor
$f_{l}$ has one set of neighboring variable nodes $V_{l}\in V$ such that
$f_{l}(V_{l})$ is a function of only those variables in $V_{l}$. The joint
distribution over the graph is then proportional to the global function
$f(V)$, factorized by $f_{l}(V_{l})$:

$p(V)\propto f(V)=\prod_{l}f_{l}(V_{l})$ (4)

For a more complete description of factor graphs and algorithms for different
operations on the graphs see [13]. For factor graphs in the context of
robotics and graph-SLAM see [14].

An example of a factor graph describing a simple tracking problem, where the
variables $v_{m}$ describe the target position states in time $x_{k}$, can be
seen in Fig. 2a. Here we have three types of factors: a prior, a measurement
likelihood and a dynamic transition factors. The posterior distribution, given
the measurements, is then proportional to the factorization of the global
function:

$\begin{split}p(x_{k:0}|y_{k:1})\propto
f(x_{0})\prod_{k=1}^{K}f(x_{k};y_{k})f(x_{k}|x_{k-1}),\end{split}$ (5)

where we use $f(x_{k}|x_{k-1})=f(x_{k},x_{k-1})\propto p(x_{k}|x_{k-1})$ and
$f(x_{k};y_{k})=f(x_{k})\propto l(x_{k};y_{k})$ to express and emphasize
conditional probabilities between variables and between variables to
measurements, respectively. This is based on the idea in [15] of extending
factor graphs to represent directed (Bayesian networks) and undirected (Markov
random fields), thus making dependencies explicit when reading the factor
graph.
As we show in more detail later, in DDF, fusion is done by incorporating a new
factor into the local graph, representing the new information that agent $i$
received from agent $j$.

###  IV-A Factor Graphs in Gaussian Canonical Form

A factor graph $F=(U,V,E)$ represents the distribution $p(V)$ by a product of
its factors (see (4)). This can be transformed to a summation by converting
the factors to log-space as described in [16]. A special case is the canonical
(information) form of linear Gaussian distributions; factors then consists of
two elements, namely the information vector and matrix, that in turn can be
summed over the graph to describe the multivariate Gaussian distribution in
canonical form,

$p(V)\sim\mathcal{N}(\zeta,\Lambda)=\sum_{l}f_{l}(V_{l}).$ (6)

Here $\zeta$ and $\Lambda$ are the information vector and matrix,
respectively. In this case all factors are Gaussian distributions represented
in canonical form $f_{l}(V_{l})\sim\mathcal{N}(\zeta_{l},\Lambda_{l})$.

###  IV-B Probabilistic Operations on Factor Graphs

In [17] a graphical model perspective of the SLAM problem is discussed. In
order to reduce the size of the graph a filtering approach is taken by
defining three types of operations on the dynamic Bayes net (DBN) type graph.
The three operations are: _prediction_ , _roll-up_ (marginalization of past
states) and _estimation_ (measurement update). We follow [17] and the
information augmented state (_iAS_) filter presented in [5] to define these
operations on a factor graph and show their translation into new factors.
_Prediction:_
In the prediction step three factors are added to the graph: two unary factors
$f^{-}(x_{k}),f^{-}(x_{k+1})$, connected to the variable nodes $x_{k}$ and
$x_{k+1}$, respectively, and a binary factor $f^{-}(x_{k+1}|x_{k})$ connected
to both variables and describes the correlation between the two variables. The
$(-)$ superscript denotes prediction. Each factor is defined by an information
vector and matrix $\\{\zeta,\Lambda\\}$:

$\begin{split}f^{-}(x_{k})&=\\{-F_{k}^{T}Q_{k}^{-1}G_{k}u_{k},\
F_{k}^{T}Q_{k}^{-1}F_{k}\\}\\\ f^{-}(x_{k+1})&=\\{Q_{k}^{-1}G_{k}u_{k},\
Q_{k}^{-1}\\}\\\ f^{-}(x_{k+1}|x_{k})&=\left\\{\begin{pmatrix}0_{n\times 1}\\\
0_{n\times 1}\end{pmatrix},\begin{pmatrix}0_{n\times n}&-Q_{k}^{-1}F_{k}\\\
-F_{k}^{T}Q_{k}^{-1}&0_{n\times n}\end{pmatrix}\right\\}.\end{split}$ (7)

Here $F_{k}$ and $G_{k}$ are the state transition and control matrices at time
step $k$, respectively. $u_{k}$ is the input vector and $Q_{k}$ is a zero mean
white Gaussian process noise covariance matrix. The graphical description of
the prediction step is given in Fig. 2b.
_Roll-up (marginalization):_
It is known that upon marginalization of a random variable $x$ in a directed
graphical model, all variables in the Markov blanket of $x$ are moralized,
i.e. “married” by adding an edge. The effect in a factor graph is similar, and
variables are moralized by adding a factor connecting all variables in the
Markov blanket of the marginalized variable. If we denote the Markov blanket
of $x$ by $\bar{x}$, then the new factor $f(\bar{x})$ is computed in two
steps:
1\. Sum all factors connected to $x$ to compute the following information
vector and matrix:

$\\{\zeta,\ \Lambda\\}=f(x)+\sum_{i\in\bar{x}}f(x,\bar{x}_{i})$ (8)

2\. Use Schur complement to compute the marginal and find the new factor
$f(\bar{x})$:

$f(\bar{x})=\\{\zeta_{\bar{x}}-\Lambda_{\bar{x}x}\Lambda_{xx}^{-1}\zeta_{x},\
\Lambda_{\bar{x}\bar{x}}-\Lambda_{\bar{x}x}\Lambda_{xx}^{-1}\Lambda_{x\bar{x}}\\}.$
(9)

Notice that as a result, conditionally independent variables become
correlated.
We demonstrate marginalization in Fig. 2c, marginalizing out $x_{k:n}$ induces
a new factor $f(x_{k+1})$ over $x_{k+1}$. Here, since the only variable in the
Markov blanket of the marginalized variables is $\bar{x}=x_{k+1}$, the new
factor is unary over $x_{k+1}$ alone.
_Estimation (measurement update):_
Adding a measurement in the canonical (information) form of the Kalman filter
is a simple task. In the factor graph this translates to adding a factor
$f(x_{k+1};y_{k+1})$ connected to all measured variables,

$f(x_{k+1};y_{k+1})=\\{H_{k+1}^{T}R_{k+1}^{-1}y_{k+1},\
H_{k+1}^{T}R_{k+1}^{-1}H_{k+1}\\}.$ (10)

Where $H_{k+1}$ is the sensing matrix, $R_{k+1}$ is a zero mean white Gaussian
measurement noise covariance matrix and $y_{k+1}$ is the noisy measurement
vector. Figure 2d shows the addition of a unary measurement factor
$f(x_{k+1};y_{k+1})$.

##  V Factor Graphs for DDF - FG-DDF

The shift in our approach to PGM based decentralized data fusion is that now
each agent holds its own _local_ factor graph (Fig. 3a and 3b) and does not
need to reason over a full (global) version of the graph. The fusion problem
can be viewed as exchanging factors between agents, and thus the challenge is
to account for common information, i.e. guarantee that mutual information is
not introduced more than once into the local graph. While this is a common
problem in DDF, the advantage of the graph approach, especially in the case of
heterogeneous fusion, is in utilizing the structure evident from the graph to
account for the common information to design and analyze fusion algorithms. We
now show how the fusion equations translate to a factor graph representation
and then demonstrate its applicability in analyzing the HS-CF algorithm [5].

###  V-A DDF in Factor Graphs

Taking natural logarithm and rearranging, the heterogeneous fusion given in
(3) can be written as,

$\begin{split}&\log p^{i}_{f}(\chi^{i}|Z^{i,+}_{k})=\log
p^{i}(\chi^{i}_{L}|\chi^{ij}_{C},Z^{i,-}_{k})+\log
p^{i}(\chi^{ij}_{C}|Z^{i,-}_{k})\\\ &+\log
p^{j}(\chi^{ij}_{C}|Z^{j,-}_{k})-\log p_{c}(\chi^{ij}_{C}|Z^{i,-}_{k}\cap
Z^{j,-}_{k})+\tilde{C}\end{split}$ (11)

We can divide the terms in the above equation into two contributions; the
first row corresponds to the information agent $i$ holds prior to fusion while
the second row (ignoring the normalization constant) describes the new
information agent $i$ should receive from agent $j$ regarding their common
states $X_{C}^{ij}$, assuming the common information is explicitly known. This
message can be viewed as sending a factor $f^{ji}(X_{C}^{ij})$ from agent $j$
to agent $i$, which is added into agent $i$’s factor graph as seen in Fig. 3d
and 3e.

The pseudo code for the FG-DDF Algorithm is given in Algorithm 1. This
represents the local fusion algorithm for each each agent $i$ communicating
with an agent $j$ in its neighborhood $N_{a}^{i}$. In the case of tree
structured network, if the agents are using a channel filter to track the
common information, then an additional factor graph is instantiated over the
common states for each communication channel agent $i$ has with its neighbors
(line 5). The main steps of the algorithm are the send / fuse message (lines
13-14), since these steps are specific to each fusion algorithm, we choose the
HS-CF as an example and give pseudo code in Algorithms 2 and 3, respectively.

Figure 3: DDF in factor graphs: a and b showing the local factor graphs over
each agent’s variables of interest. c) is the common graph describing the
factors common to both agents. d-e) demonstrates the fusion operation; the
message sent from $j$ with new information over the common variables (d) then
integrated into $i$’s local graph with a simple factor addition (e). Algorithm
1 FG-DDF

1:Define: $\chi_{i}$, Priors, Fusion algorithm

2:Initialize local factor graph $\mathcal{F}^{i}$

3:if CF algorithm then

4: for all $j\in N_{a}^{i}$ do

5: Initialize CF factor graph $\mathcal{F}^{ij}$ over $\chi_{C}^{ij}$

6: end for

7:end if

8:for all time steps do

9: for all $x\in\chi_{i}$ do

10: Measurement update step $\triangleright$ Eq.10

11: end for

12: for all $j\in N_{a}^{i}$ do

13: Send message to $j$ $\triangleright$ e.g., Algorithm 2

14: Fuse message from $j$ $\triangleright$ e.g., Algorithm 3

15: end for

16:end for

17:return

###  V-B Inference

In the previous sections we described how an agent builds a local factor graph
over its variables of interest through a series of predictions, local
measurements and fusion of data from neighboring agents. At this point it is
important to describe how, given the local factor graph, inference is
performed. The goal of inference in this case is to deduce the marginal pdfs
of each variable of interest from the joint function described by the factor
graph. When the factor graph is tree-structured the _sum-product algorithm_
[13], [18] can be used to directly work on the factor graph. However, when
there are more than one variable in common between two communicating agents,
the post-fusion graph has cycles due to the marginalization of the local
variables at the “sending” agent $j$.

To solve this problem we suggest transforming the factor graph into a tree by
forming cliques over the cycles. Note that we are not transforming the graph
into a clique or a junction tree, but only forming the minimum number of
cliques that will result a tree graph as we explain next. It is worth
mentioning here that Kaess et al. [19] use the elimination algorithm to
transfer a factor graph into a Bayes tree for incremental updates. This
approach might be useful in the future but is not advantageous in the scope of
this work. Instead, we keep the graph as a factor graph and only join into
cliques variables which are now not d-separated by the “local” variables and
then summarize their factors to new joint factors, connected to the clique.
This is demonstrated in Fig. 4, where assuming a tree-structured network
$j-i-k$, the common variables are separated into three different sets
$X_{C}^{ijk}$, $X_{C}^{ij\setminus k}$ and $X_{C}^{ik\setminus j}$,
representing variables common to the three agents, variables common to $i$ and
$j$ but not to $k$ and similarly variables common to $i$ and $k$ but not to
$j$, respectively. We can see that fusion results loops in the graph. We
restore the tree structure by forming a clique over the separation sets
$X_{L}^{i}\cup X_{C}^{ijk}$. Note that this needs to be done only before
inference and not for fusion or filtering.

Figure 4: Transitioning the graph into a clique factor graph demonstrated on a
network $j-i-k$. a) Local graph before fusion - tree structure. b) Messages
sent from $k$ and $j$ to $i$ results in loops in the local graph. c) a tree
structure is regained by forming cliques.

###  V-C HS-CF in Factor Graphs

To demonstrate the applicability in analyzing algorithms and track the
implicit common information in heterogeneous DDF problems we show how the
heterogeneous state channel filter (HS-CF) algorithm [5] translates to a
factor graph. Briefly, the HS-CF builds on the homogeneous channel filter (CF)
[7] and exploits conditional independence between states local to each agent,
given the common states between the agents as described in (3), to perform
heterogeneous fusion and only fuse information relating to common states of
interest. When the network is tree-structured, a channel filter is placed on
the communication channel between any two agents, tracking the information
passed through the channel. In a factor graph this translates to each agent
maintaining a _“locally full”_ factor graph over all its states of interest
$\chi^{i}$ and a separate _partial_ factor graph over the common states of
interest $\chi_{C}^{ij}$ for each communication channel, as shown in Fig. 3a-c
respectively.

Fusion over common states of interest is done in two steps, sending a message
and fusing incoming messages. When sending a message, as described in
Algorithm 2, an agent first marginalizes out local “irrelevant” variables out
of its graph (Algorithm 4), then subtracts the channel filter factors to
synthesize new information factors to send to agent $j$ and updates the CF
graph with the new factors.

Fusing an incoming message is a simple process, as described in Algorithm 3.
An agent $i$ that received a message from a neighbor agent $j$ need only to
update its local graph and CF graph with the set of incoming factors over the
common variables of interest.

Algorithm 2 Send Message to $j$ (HS-CF)

1:Input: $j$ \- agent to communicate with

2:Marginalize out local variables $\chi^{i}\setminus\chi_{C}^{ij}$
$\triangleright$ Algorithm 4

3:Subtract factors of CF factor graph from factors of marginalized graph to
receive a set of new information factors
$f\in\boldsymbol{f}^{ij}(\chi_{C}^{ij})$

4:Update CF factor graph with all $f_{n}\in\boldsymbol{f}^{ij}(\chi_{C}^{ij})$

5:return Message $\boldsymbol{f}^{ij}(\chi_{C}^{ij})$

Algorithm 3 Fuse Message from $j$ (HS-CF)

1:Input: Message from agent $j$ \- $\boldsymbol{f}^{ji}(\chi_{C}^{ij})$

2:for all $f\in\boldsymbol{f}^{ji}(\chi_{C}^{ij})$ do

3: Add $f$ to local factor graph $\mathcal{F}^{i}$

4: Add $f$ to CF factor graph $\mathcal{F}^{ij}$

5:end for

6:return

Algorithm 4 Marginalize variable $x$ from a factor graph $\mathcal{F}$

1:Input: variable to remove $x$, factor graph $\mathcal{F}$

2:Sum $f(x)$ and all factors $f(x,\bar{x}_{i})$ adjacent to $x$

3:Create new marginal factor $f(\bar{x})$ using (9)

4:Add edges from every $\bar{x}_{i}\in\bar{x}$ to $f(\bar{x})$

5:Remove $f(x,\bar{x}_{i})$ factor and $x$ from $\mathcal{F}$

6:return marginal factor graph $\mathcal{F}$

##  VI Simulation

To validate the algorithms and demonstrate our approach for decentralized
fusion problems we perform two different simulations. First we validate the
algorithms with a multi-agent multi-target tracking simulation and compare it
to previous simulation results, which were implemented using conventional KF-
based filtering and HS-CF equations as described in [5]. Then we test the
algorithm on a multi-agent mapping problem and show its applicability and
advantages in computation and communication reduction, thus making it scalable
and attractive for large networks. The algorithms are written in Python and
are based on the _fglib_ [20] and _NetworkX_ [21] open source libraries

###  VI-A Example 1 - Target tracking

As a validation case, we repeat and compare to the simulation scenario
described in [5], and so for brevity, we will not repeat all the details here.
We consider a 2D static five-agent six-target tracking problem, where the
agents are connected bi-laterally in ascending order. Each agent $i$ inference
task is to estimate it own constant relative position measurement bias
(assuming perfect local position information) and the position of a set of
$n_{t}^{i}$ targets. The agent’s target assignments and sensor measurement
error covariances are given in Table I. In the PGM heterogeneous formulation
of the problem, each agent $i$ maintains a local factor graph over
$\chi^{i}=\cup_{t\in T^{i}}x_{t}^{i}\bigcup s^{i}$, where
$x_{t}^{i}=[e_{t},n_{t}]^{T}$ are the east and north positions of target $t$
in its target assignments $T^{i}$, and $s^{i}=[b_{e}^{i},b_{n}^{i}]^{T}$ is
the agent’s own east and north biases. In the homogeneous formulation, each
agent has to process and communicate the full set of variables (22 in this
example). In heterogeneous fusion, it processes and communicates only a subset
of the full system variables: maximum 8 and 4, respectively (e.g., agent 3).
This leads to about $2.6\%$ and $4.8\%$ of the communication and computation
requirements, respectively, of the full state homogeneous fusion.

Fig. 5 validates our implementation by comparing the results of 50 Monte-Carlo
simulations with the results from [5]. We show that the mean RMSE (full line)
and the $2\sigma$ lines of our new code FG-DDF, which has all its operations
defined and performed on a factor graph, exactly matches the results of the
HS-CF algorithm. Also plotted is a comparison with the centralized solution,
showing the estimate is conservative, as it produces estimation error
covariances which are larger in the positive semi-definite (psd) sense than
the centralized (optimal) error covariance.

TABLE I: Local platform target assignments and sensor measurement error covariances, taken from [5].  Agent | Tracked Targets | $R_{i}^{1}[m^{2}]$ | $R_{i}^{2}[m^{2}]$
---|---|---|---
1 | $T_{1},T_{2}$ | diag([1,10]) | diag([3,3])
2 | $T_{2},T_{3}$ | diag([3,3]) | diag([3,3])
3 | $T_{3},T_{4},T_{5}$ | diag([4,4]) | diag([2,2])
4 | $T_{4},T_{5}$ | diag([10,1]) | diag([4,4])
5 | $T_{5},T_{6}$ | diag([2,2]) | diag([5,5])
Figure 5: Validation of the FG-DDF algorithm. Shown are the RMSE results and
2$\sigma$ bounds of agents 2 and 3 over their local variables (with only
target 3 in common), compared to a centralized estimate and the HS-CF
algorithm from [5].

###  VI-B Example 2 - Cooperative Mapping

We test our algorithm with a simulation of a 2D multi-agent cooperative
mapping problem. In this scenario, four agents are moving in a
circular/elliptical trajectory in a $180m\times 160m$ field. The field is
subdivided into four separate but overlapping sections and each agent is
tasked with mapping one of them (Fig. 6a). Here mapping refers to the task of
inferring the $2D$ position of a subset of 25 known landmarks, scattered
across the field. We assume each agent has perfect self position information,
but with constant agent-landmark relative position measurement bias in the
east and north direction. Further, we assume that an agent takes a noisy
relative measurement only if the landmark is within $45m$ radius from the
agent and that the data association problem is solved. In addition, each agent
takes a noisy measurement to estimate its own bias.

In homogeneous DDF, each agent would need to reason about the full set of 25
landmarks, regardless of its relative distance and relevance to the local task
as well as the other agents local bias state, resulting in 58 variables. Those
variables then need to be processed and communicated repeatedly, incurring
high overhead cost. In the heterogeneous FG-DDF approach, each agent only
processes the locally relevant landmarks and biases (maximum 20) and
communicates only information over the common variables (maximum 6), which in
turns lead to about $96\%$ and $99\%$ computation and communication reduction,
respectively.

Figure 6: Multi-agent cooperative mapping illustration and simulation results.
Agent 2 estimate and $2\sigma$ covariance bounds of a) landmark 14, and b)
local bias. c) Global map and task subdivision with shaded overlapping
sections. Black x’s mark true landmark position and colored x’s mark agent’s
estimate (at $t=150s$) with $2\sigma$ covariance ellipses.

Fig. 6a illustrates agent’s trajectories and landmarks estimates after 150s
with their $2\sigma$ confidence ellipses. Each agent only estimates its local
map (subset of 6-9 landmarks). Fig. 6b and 6c show agent 2 estimate vs. time
of its local bias and landmark 14, respectively. It is worth exploring the
dynamics of the landmark’s estimate of this case as we can observe three
regions of the covariance bounds: 1) in the first 5s the agent is taking
relative measurements to the landmark, allowing it to reduce its uncertainty
to about 13.5m (two standard deviation). 2) In the time frame 5s-27s the agent
has no direct new information on the landmark as it is out of measurement
range, but it can nevertheless marginally improve its estimate to about 12.5m
due to bias measurements and its correlation to the landmark state. 3) Agent
1, which has landmark 14 in common with agent 2, starts taking relative
measurements to the landmark, agent 2 is then able to further improve its
estimate by fusion with agent 1 to about 3m at $t=85$. This cycle repeats as
the two agents get in and out of measurement range to landmark 14, but the
changes are small.

##  VII Conclusions

In heterogeneous robotic systems, heterogeneous decentralized data fusion is
key to enable data sharing and scalability for multi-agent real world task
oriented applications such as cooperative mapping, target tracking and more.
In this paper we merge the Bayesian decentralized data fusion (DDF) approach,
which sees each agent as an independent entity, inferring the variable of
interest based on its local information, with probabilistic graphical models,
which frequently assumes knowledge of the full system graph and correlations.
We show that factor graphs are a natural representation for heterogeneous
fusion in robotic applications as it gives insight into the problem structure.
This enables exploitation of conditional independence between sets of local
and common variables of interest. We then suggest a shift in paradigm by
slicing the global factor graph into smaller local graphs representing its
individual, though cooperative, task.

We assume linear Gaussian models and constant tree structured networks to
develop the theory for factor graph based DDF (FG-DDF), including basic
filtering, fusion and inference operations, performed solely on the graph. We
show how, by using this framework, fusion of information from a neighboring
agent reduces into simply integrating a new factor into the local graph. We
validate our algorithm by comparing to previous results of a static multi-
agent multi-target simulation and further demonstrate its applicability with a
cooperative mapping problem. In both cases, the results are consistent and
conservative ((i.e. produce estimation error covariances which are larger in
psd sense than the true error covariances)) and the advantages of
heterogeneous fusion in terms of computational and communication scalability
are highlighted.

Our approach enables further development and analysis of new fusion algorithms
in more complicated scenarios, based on this relatively simple and general
representation. While we used some simplifying assumptions in order to develop
the theory of FG-DDF, it is not limited to Gaussian pdf nor to tree-structured
networks (e.g. covariance intersection could be used in place of the CF here
[22]). These avenues, as well as dynamic problems are the subject of future
research.

## References

  * [1] M. Paskin and C. E. Guestrin, “Robust Probabilistic Inference in Distributed Systems,” _arXiv:1207.4174 [cs]_ , 2004, arXiv: 1207.4174. [Online]. Available: http://arxiv.org/abs/1207.4174
  * [2] A. Cunningham, V. Indelman, and F. Dellaert, “DDF-SAM 2.0: Consistent distributed smoothing and mapping,” in _2013 IEEE International Conference on Robotics and Automation_ , May 2013, pp. 5220–5227, iSSN: 1050-4729.
  * [3] C. Castello, J. Fan, A. Davari, and R.-X. Chen, “Temperature Control Framework Using Wireless Sensor Networks and Geostatistical Analysis for Total Spatial Awareness,” in _2009 10th International Symposium on Pervasive Systems, Algorithms, and Networks_ , Dec. 2009, pp. 717–721, iSSN: 2375-527X.
  * [4] J. Pearl, _Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference_. Elsevier, Jun. 2014, google-Books-ID: mn2jBQAAQBAJ.
  * [5] O. Dagan and N. R. Ahmed, “Exact and Approximate Heterogeneous Bayesian Decentralized Data Fusion,” _arXiv:2101.11116 [cs, eess]_ , Jan. 2021, arXiv: 2101.11116. [Online]. Available: http://arxiv.org/abs/2101.11116
  * [6] A. Makarenko, A. Brooks, T. Kaupp, H. Durrant-Whyte, and F. Dellaert, “Decentralised data fusion: A graphical model approach,” in _2009 12th International Conference on Information Fusion_ , Jul. 2009, pp. 545–554.
  * [7] S. Grime and H. Durrant-Whyte, “Data fusion in decentralized sensor networks,” _Control Engineering Practice_ , vol. 2, no. 5, pp. 849–863, Oct. 1994.
  * [8] S. K. Andersen, K. G. Olesen, F. V. Jensen, and F. Jensen, “HUGIN*— a Shell for Building Bayesian Belief Universes for Expert Systems,” _IJCAI_ , vol. 89, p. 6, 1989.
  * [9] C.-Y. Chong and S. Mori, “Graphical Models for Nonlinear Distributed Estimation,” in _Proceedings of the 7th International Conference on Information Fusion_ , Stockholm, Sweden, 2004, pp. 614–621.
  * [10] F. Dellaert, “Factor Graphs: Exploiting Structure in Robotics,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 4, no. 1, pp. 141–166, 2021, _eprint: https://doi.org/10.1146/annurev-control-061520-010504. [Online]. Available: https://doi.org/10.1146/annurev-control-061520-010504
  * [11] B. Etzlinger, F. Meyer, F. Hlawatsch, A. Springer, and H. Wymeersch, “Cooperative Simultaneous Localization and Synchronization in Mobile Agent Networks,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 14, pp. 3587–3602, Jul. 2017, conference Name: IEEE Transactions on Signal Processing.
  * [12] C. Y. Chong, E. Tse, and S. Mori, “Distributed Estimation in Networks,” in _1983 American Control Conference_ , Jun. 1983, pp. 294–300.
  * [13] B. J. Frey, F. R. Kschischang, H.-A. Loeliger, and N. Wiberg, “Factor graphs and algorithms,” in _In Proc. 35th Allerton Conf. Communications, Control, and Computing_ , 1997, pp. 666–680.
  * [14] F. Dellaert and M. Kaess, _Factor Graphs for Robot Perception_ , ser. Foundations and Trends(r) in Robotics Series. Now Publishers, 2017. [Online]. Available: https://books.google.com/books?id=BAFbswEACAAJ
  * [15] B. J. Frey, “Extending factor graphs so as to unify directed and undirected graphical models,” in _Proceedings of the Nineteenth conference on Uncertainty in Artificial Intelligence_ , ser. UAI’03. Acapulco, Mexico: Morgan Kaufmann Publishers Inc., Aug. 2002, pp. 257–264.
  * [16] D. Koller and N. Friedman, _Probabilistic graphical models: principles and techniques_ , ser. Adaptive computation and machine learning. Cambridge, MA: MIT Press, 2009.
  * [17] M. A. Paskin, “Thin Junction Tree Filters for Simultaneous Localization and Mapping,” University of California, Berkeley, Tech. Rep. UCB/CSD-02-1198, 2002.
  * [18] F. Kschischang, B. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” _IEEE Transactions on Information Theory_ , vol. 47, no. 2, pp. 498–519, Feb. 2001, conference Name: IEEE Transactions on Information Theory.
  * [19] M. Kaess, V. Ila, R. Roberts, and F. Dellaert, “The Bayes Tree: An Algorithmic Foundation for Probabilistic Robot Mapping,” in _Algorithmic Foundations of Robotics IX: Selected Contributions of the Ninth International Workshop on the Algorithmic Foundations of Robotics_ , ser. Springer Tracts in Advanced Robotics, D. Hsu, V. Isler, J.-C. Latombe, and M. C. Lin, Eds. Berlin, Heidelberg: Springer, 2011, pp. 157–173. [Online]. Available: https://doi.org/10.1007/978-3-642-17452-0_10
  * [20] D. Bartel, “danbar/fglib,” Mar. 2021, original-date: 2013-11-04T08:43:18Z. [Online]. Available: https://github.com/danbar/fglib
  * [21] A. A. Hagberg, D. A. Schult, and P. J. Swart, “Exploring Network Structure, Dynamics, and Function using NetworkX,” in _Proceedings of the 7th Python in Science Conference (SciPy2008)_ , Pasadena, CA, 2008, p. 5.
  * [22] I. Loefgren, N. Ahmed, E. Frew, C. Heckman, and S. Humbert, “Scalable Event-Triggered Data Fusion for Autonomous Cooperative Swarm Localization,” in _2019 22th International Conference on Information Fusion (FUSION)_ , Jul. 2019, pp. 1–8.