# Improved analytical modeling of the non-linear power spectrum in modified
gravity cosmologies

Suhani Gupta  [gupta, hellwing<EMAIL_ADDRESS>Wojciech A. Hellwing
Maciej Bilicki  Center for Theoretical Physics, Polish Academy of Sciences,
Al. Lotników 32/46, 02-668 Warsaw, Poland

###### Abstract

Reliable analytical modeling of the non-linear power spectrum (PS) of matter
perturbations is among the chief pre-requisites for cosmological analyses from
the largest sky surveys. This is especially true for the models that extend
the standard general-relativity paradigm by adding the fifth force, where
numerical simulations can be prohibitively expensive. Here we present a method
for building accurate PS models for two modified gravity (MG) variants: namely
the Hu-Sawicki $f(R)$, and the normal branch of the Dvali-Gabadadze-Porrati
(nDGP) braneworld. We start by modifying the standard halo model (HM) with
respect to the baseline Lambda-Cold-Dark-Matter ($\mathrm{\Lambda\text{CDM}}$)
scenario, by using the HM components with specific MG extensions. We find that
our $P(k)_{\text{HM}}$ retains 5% accuracy only up to mildly non-linear scales
($k\lesssim 0.3$ $h/\,\mbox{Mpc}$) when compared to PS from numerical
simulations. At the same time, our HM prescription much more accurately
captures the ratio $\Upsilon(k)=P(k)_{\text{MG}}/P(k)_{\Lambda\text{CDM}}$ up
to non-linear scales. We show that using HM-derived $\Upsilon(k)$ together
with a viable non-linear $\mathrm{\Lambda\text{CDM}}$ $P(k)$ prescription
(such as halofit), we render a much better and more accurate PS predictions in
MG. The new approach yields considerably improved performance, with modeled
$P(k)_{\text{MG}}$ being now accurate to within 5% all the way to non-linear
scales of $k\lesssim 2.5-3$ $h/\,\mbox{Mpc}$. The magnitude of deviations from
GR as fostered by these MG models is typically $\mathcal{O}(10\%)$ in these
regimes. Therefore reaching 5% PS modeling is enough for forecasting
constraints on modern-era cosmological observables.

##  I Introduction

The standard model of cosmology, the Lambda-Cold-Dark-Matter
($\mathrm{\Lambda\text{CDM}}$), has been remarkably well-tested
observationally in the last two decades. Presently, it is our best
approximation of the real Universe (Eisenstein  _et al._ , 2005;
Collaboration, 2020; Cyburt  _et al._ , 2016; Alam  _et al._ , 2021; Betoule
_et al._ , 2014). The precise observations of the Cosmic Microwave Background
(CMB) radiation Collaboration (2020); Hinshaw  _et al._ (2013), large-scale
galaxy clustering Alam  _et al._ (2021); Cole  _et al._ (2005); Blake  _et
al._ (2011); Pezzotta  _et al._ (2017), and the abundance of massive galaxy
clusters Planck Collaboration (2014) among others, form a long list where the
standard cosmological model predictions are successful.

So far, the bulk of these $\mathrm{\Lambda\text{CDM}}$ observational tests
concerns the linear regime, the large scales, and/or early times. But it is
the mildly non-linear and fully non-linear density fluctuation regimes where
the vast majority of the modifications to $\mathrm{\Lambda\text{CDM}}$ are
expected to deviate significantly from the standard model predictions Bean and
Tangmatitham (2010); Hellwing  _et al._ (2017); Winther  _et al._ (2017); Bose
_et al._ (2017); Barreira  _et al._ (2015); Koyama  _et al._ (2009). It is
also in this regime, stretching usually from hundreds down to a few
Megaparsecs, where the present and upcoming cosmological surveys like desi
DESI Collaboration (2016), lsst LSST Science Collaboration (2009), and euclid
Laureijs  _et al._ (2011) aim to measure various statistics concerning the
large-scale structure to a percent level accuracy. As a result, with the
influx of data from these surveys, the level of the statistical errors can get
so small that the measurements start to be more sensitive to systematic
effects. If both the new level of accuracy of base-level predictions as well
as the control of the known systematics will be successfully implemented,
these new large-scale surveys will yield new unprecedentedly accurate
estimates and constraints on cosmological parameters: like the DE equation of
state, the growth rate of structure, or parameters quantifying possible
departures from the standard General Relativity (GR)-based structure formation
scenario.

In this context, one of the most useful and widely used theoretical quantities
is the power spectrum (PS) of density fluctuations, $P(k)$. This statistic
generally characterizes the properties of large-scale structures across vast
cosmological epochs and scales. Not only it can be used as an end-goal model
prediction on its own, but it is also a basic quantity that is used to model
and forecast a number of other useful LSS observables, including galaxy
clustering measures, cluster abundance, weak-lensing shear and convergence,
the amplitude of the bulk peculiar galaxy motions, and many others Eisenstein
_et al._ (2005); Kilbinger  _et al._ (2013); Hikage  _et al._ (2019); LSST
Science Collaboration (2009); Alam  _et al._ (2021); To  _et al._ (2021);
Turner  _et al._ (2023).

Since the PS forms a basis for the predictions of many cosmological LSS
observational statistics, the accuracy, and scales to which we know the input
PS limits our abilities to forecast the derived observables. Thus, obtaining
robust estimates of the PS beyond the linear theory regime (i.e. scales of
$k>0.1$ $h/\,\mbox{Mpc}$) became of paramount importance to modern cosmology,
and has been a subject of extensive effort in recent years Cooray and Sheth
(2002); Smith  _et al._ (2003); Mead  _et al._ (2015); Lawrence  _et al._
(2010). A classical approach is to either directly use the results of a number
of N-body simulations of LSS evolution to predict PS, or use the simulation
results for calibration of more or less approximate models (Cooray and Sheth,
2002; Smith  _et al._ , 2003; Takahashi  _et al._ , 2012; Mead  _et al._ ,
2015; Smith and Angulo, 2019). Recently, machine-learning-based-emulators are
also employed for computing non-linear PS (e.g. Lawrence  _et al._ , 2010;
Angulo  _et al._ , 2021). This approach especially depends on the growing
computational power.

In recent years, the progress in modeling the PS has been truly significant.
The resulting current state-of-the-art PS models for
$\mathrm{\Lambda\text{CDM}}$ are already, or close to, attaining sub-percent
accuracy in the non-linear regime, as required for the success of the
cosmological tests offered by the incoming big survey data. However, this
amazing progress has been mostly limited to the $\mathrm{\Lambda\text{CDM}}$
alone. When it comes to many interesting extensions and modifications of the
standard model, such as the whole family of beyond-GR Modified Gravity (MG)
scenarios, the current accuracy, and versatility of PS modeling is still very
much lacking. The reasons for this are both higher theoretical complications
of such models, and their increased levels of non-linearity Oyaizu  _et al._
(2008); Alam  _et al._ (2021); Hu and Sawicki (2007a); Koyama  _et al._
(2009); Cataneo  _et al._ (2019); Ruan  _et al._ (2020). For MG models, N-body
simulations play an even more important role in fully assessing the effect of
the fifth force, and are crucial for disentangling pure MG effects from the
standard GR-based scenarios Baldi  _et al._ (2014); Winther  _et al._ (2015);
Puchwein  _et al._ (2013). This is connected with the richer phenomenology of
such models Clifton  _et al._ (2012); Baker  _et al._ (2015); Koyama (2016);
Nojiri  _et al._ (2017); Ferreira (2019); Sotiriou and Faraoni (2010); Baker
_et al._ (2021); Koyama (2018); Joyce  _et al._ (2016); Bull  _et al._ (2016).
Given the fact that the MG simulations are usually many times more expensive
than the standard $\mathrm{\Lambda\text{CDM}}$ case Barreira  _et al._ (2015);
Bose  _et al._ (2017), it becomes computationally prohibitive to obtain
simulation libraries of the same volume and precision for MG, as is possible
for $\mathrm{\Lambda\text{CDM}}$. However, such libraries are necessary to be
applied to the proven state-of-the-art emulating or fitting methods to achieve
the same precision, and success in modeling MG effects, as we have for the
case of $\mathrm{\Lambda\text{CDM}}$.

In this work, we attempt to remedy the deficit of accurate MG PS modeling. To
circumnavigate the problem of prohibitively expensive MG simulations, we
explore a different approach. Instead of trying to model the absolute MG PS
predictions, we take $\mathrm{\Lambda\text{CDM}}$ to always be our baseline,
and build a semi-analytical model for the relative MG effects on the
$\mathrm{\Lambda\text{CDM}}$ PS. We build our model on the basis of a more
general Halo Model (HM) approach (Seljak (2000), reviewed by (e.g. Cooray and
Sheth, 2002)). Next, we demonstrate how various degrees of modeling freedom
can be calibrated and constrained already by a relatively small library of
N-body simulations, to achieve an unprecedented level of PS modeling in the MG
scenarios studied here.

There are many models that can be considered beyond-GR structure formation
scenarios. Most of the viable, and at the same time cosmologically interesting
ones usually involve some extra couplings to the metric in the Einstein-
Hilbert action that manifests themselves as additional degrees of freedom
(d.o.f.). The propagation (gradient) of this d.o.f. induces an additional
gravitational force component, called as the fifth force, which acts on top of
the Newtonian gravitational force on the cosmological scales. However,
propagation of a significant fifth force both on small galaxy scales, and in
the strong field regime is tightly constrained observationally Will (2014);
Taylor and Weisberg (1982); Abbott  _et al._ (2017, 2016); Gravity
Collaboration  _et al._ (2020); Shapiro (1964). Thus, only MG models that
exhibit some kind of a fifth force screening mechanism, which, as the name
suggests, would screen the fifth force in these observationally tested regimes
are viable MG candidates Joyce  _et al._ (2015); Koyama (2016); Brax (2013);
Sakstein (2018); Koyama (2018).

The clockwork of MG models and their involved screening mechanisms can differ
in many ways. From our point of view, however, we can significantly simplify
the subject by focusing just on phenomenological effective modifications to
the density fluctuations PS. As our test-case models, we choose variants of
two popular MG set-ups: namely $f(R)$ Sotiriou and Faraoni (2010) and nDGP
gravity Dvali  _et al._ (2000), which will serve as a good representative for
their whole respective families. Further in the text, we offer a more detailed
description and definitions of these models.

Most of the works that have considered computing the non-linear PS in MG
models either rely on simulations (Oyaizu  _et al._ , 2008; Alam  _et al._ ,
2021), post-Friedmann (PPF) formalism Hu and Sawicki (2007a), or perturbation
theory focusing on quasi-linear scales (Koyama  _et al._ , 2009; Cataneo  _et
al._ , 2019). In Ref. Mead  _et al._ (2016), the non-linear PS is computed
using the hmcode Mead  _et al._ (2015) for a variety of extensions to the
standard cosmological model, including $f(R)$ and nDGP. The level of this
prediction is however significantly limited by a number of approximations. For
example, a simplified spherical collapse theoretical formalism is used there
to estimate DM halo properties. From another perspective, mg-halofit was
proposed in (Zhao, 2014) as an extension of standard halofit for $f(R)$
gravity models, but Tessore  _et al._ (2015) showed that the former has
limited applicability and accuracy.

The HM formalism has been used to model non-linear MG PS in (Achitouv  _et
al._ , 2016; Lombriser  _et al._ , 2014; Ruan  _et al._ , 2020; Barreira  _et
al._ , 2014, 2014; Cataneo  _et al._ , 2019; Schmidt  _et al._ , 2010), which
is mainly based on the theoretical spherical collapse model, and is explicitly
solved for each MG variant. In our approach, however, we rely on the
calibration of phenomenological components of HM to N-body simulations. An
additional strength of our approach is that it is general enough to be quite
straightforwardly extended, not only to a wider part of the model parameters
space but also, in principle, to other modified structure formation models.

This paper is organized as follows: In Section II, we describe the MG models,
numerical data sets, and simulations. In Section III, we elaborate on the HM
formalism and describe the empirical halo properties: halo mass function
(III.1), halo bias (III.2) and halo density profile (III.3). In Section IV, we
discuss the results obtained from extending the standard HM predictions to the
MG models considered in this work (IV.1), and from our new approach (IV.2). In
(IV.3), we test our approach on another suite of MG simulations, and the final
Section V includes our conclusions, discussion, and future work prospects.
Details of the Appendices are mentioned in the respective sub-sections.

##  II Modified gravity models, numerical data sets and tools

As our main data for calibration of the non-linear PS amplitude, we take the
elephant (Extended LEnsing PHysics using ANalytic ray Tracing) suite of N-body
simulations Alam  _et al._ (2021). These simulations provide a good test bed
to model the impact of $f(R)$ and nDGP physics on formation of the large-scale
structure.

In $f(R)$, the fifth force is manifested as a result of additional degrees of
freedom from the interaction between an auxiliary scalar field (or scalaron)
and matter. This additional force appears as a non-linear function of the
Ricci scalar, $R$ in the Einstein-Hilbert action, hence the term $f(R)$. We
work with the Hu-Sawicki form of $f(R)$ gravity, where Chameleon screening
screens the fifth force Hu and Sawicki (2007b). In this screening, the
scalaron becomes very massive in the high curvature (and high matter density)
regimes, and the fifth force exponentially decays above the length scale
determined by the inverse of the mass of the scalaron. This length scale is
termed the Compton wavelength. As a result of this decay, the scalar
interaction diminishes above the Compton wavelength, and GR is recovered
Khoury and Weltman (2004).

In the nDGP model, gravity, unlike other standard forces, mediates from 4D
brane to 5D Minkowski spacetime Dvali  _et al._ (2000); Deffayet  _et al._
(2002). In this model, the scalar is identified as the brane-bending mode
which describes the deformation of the 4D brane in the 5D bulk spacetime. The
brane bending mode has a second-order term in the equation of motion. On small
scales, this term dominates over the linear term. As a result, the coupling
between the scalar field and matter is suppressed, and the solutions for
metric perturbations approach GR. This is referred to as the Vainshtein
screening Vainshtein (1972).

In elephant, along with $\mathrm{\Lambda\text{CDM}}$, two $f(R)$ variants have
been employed, with their free parameter $|f_{R0}|$ (the strength of the
scalar field today), taken to be 10-6 and 10-5 (increasing order of deviation
from $\mathrm{\Lambda\text{CDM}}$) dubbed as F6 and F5, respectively. For nDGP
gravity, we have two variants with the model parameter rcH0 = 5 and 1 (which
is the dimensionless crossing-over scale characterizing transition from 4D to
5D gravity), marked consequently as N5 and N1, respectively.

The simulations were run from $z_{\text{ini}}=49$ to $z_{\text{fin}}=0$
employing the ecosmog code (Li  _et al._ , 2012; Bose  _et al._ , 2017; Li
_et al._ , 2013a), each using $1024^{3}$ $N$-body particles in a cubic box of
a size 1024 $\mbox{Mpc}/h$. The mass of a single particle is
$m_{p}=7.798\times 10^{10}$ M${}_{\odot}/h$, and the comoving force resolution
is $\varepsilon=15$ $\mbox{kpc}/h$. Each set of simulations has five
independent realizations, evolved from the same set of initial conditions. The
cosmological parameters of the fiducial background model are given as
$\Omega_{m}$= 0.281 (fractional matter density), $\Omega_{b}$ = 0.046
(fractional baryonic density), $\Omega_{\Lambda}$ = 0.719 (fractional
cosmological constant density), $\Omega_{\nu}$ = 0 (relativistic species
density), $h$ = 0.697 (dimensionless Hubble constant), $n_{s}$ = 0.971
(primordial spectral index), and $\sigma_{8}$ = 0.842 (power spectrum
normalization). These parameters apply to background cosmologies in the
simulations of all the gravity models. For further processing, we take
simulation snapshots saved at $z=0,0.3,0.5$ and $1$.

As indicated above, the elephant-suite will be our main calibration data set.
To test the accuracy of our PS modeling and the general quality of
extrapolation, we also use different N-body data. For these additional tests,
we take the MG simulations for F5 and N1, described in Naidoo  _et al._
(2022). These simulations have background cosmological parameters different
from our parent elephant simulations, with $\Omega_{m}$ = 0.3111, $\Omega_{b}$
= 0.049, $\Omega_{\Lambda}$ = 0.6889, $\Omega_{\nu}$ = 0, $h$ = 0.6766,
$n_{s}$ = 0.9665 and $\sigma_{8}$ = 0.8245. This simulation set is run using
mg-cola Winther  _et al._ (2017) in a 500 $\mbox{Mpc}/h$ box. For each model,
we build an ensemble based on five independent realizations.

Linear matter power spectra, $P(k)_{\text{lin}}$, used in this work were
calculated using a modified version of the camb cosmological code Lewis and
Challinor (2011), which includes a module implementing both the $f(R)$ and
nDGP models. The simulation power spectra, $P(k)_{\text{sim}}$, were computed
using powmes Colombi and Novikov (2011). In what follows, by $P(k)$ we will be
denoting the fully non-linear matter power spectrum, unless indicated
otherwise.

##  III Halo model formalism

As a baseline prediction and our starting point, we take the halo model (HM)
approach. It has been proposed as an attempt to analytically model the
variance of density fluctuations into the non-linear regime using the
properties and clustering of halos as main input parameters. HM describes the
statistics of the density field up to the mildly non-linear regimes (i.e.
$k\stackrel{{\scriptstyle<}}{{{}_{\sim}}}0.5$ $h/\,\mbox{Mpc}$). Despite its
inferior accuracy compared to heavy N-body simulations, the HM has been
successfully used for modeling observables and constraining cosmological
parameters (Cooray and Sheth, 2002; Seljak, 2000; Peacock and Smith, 2000).

In HM, the main presumption is that all contributions to the cosmic density
field variance come from the matter collapsed into halos. This allows for
moderately accurate modeling of the non-linear two-point clustering
statistics, although HM can be used to compute the density field at even
higher levels of the $n$-point hierarchy Cooray and Hu (2001).

Following HM, the total matter power spectrum $P(k)_{\text{HM}}$ can be
described as a sum of two contributions:

$P(k)_{\text{HM}}=P(k)_{1h}+P(k)_{2h}$ (1)

where P$(k)_{1h}$ models the contribution from the matter clustered inside
halos (called the one-halo term) and P$(k)_{2h}$ is the contribution from
clustering of separate halos (the two-halo term). In practice, the one-halo
term dominates at small scales (i.e.
$k\stackrel{{\scriptstyle>}}{{{}_{\sim}}}1$$h/\,\mbox{Mpc}$) and saturates to
a constant value at larger scales, where the two-halo term becomes the
dominant component of the power spectrum.

These contributions are further defined as:

$P(k)_{1h}=\int_{0}^{\infty}dM|\tilde{u}(k|M)|^{2}\left(\frac{M}{\bar{\rho}}\right)^{2}n(M)$
(2)

and

$P(k)_{2h}=I_{m}^{2}(k)P(k)_{\text{lin}}$ (3)

where,

$I_{m}(k)=\frac{1}{\bar{\rho}}\int_{0}^{\infty}dM|\tilde{u}(k|M)|Mn(M)b(M)$
(4)

and $I_{m}\rightarrow 1$ for $k\rightarrow 0$ in order to match the linear
theory predictions at large scales. The integrals in Eq. 2 and Eq. 4 should in
principle cover all possible halo mass ranges, but in practice, some
M${}_{\text{min}}$ and M${}_{\text{max}}$ are introduced (these mass limits
are discussed in more details in the next sub-sections).

Here, $\bar{\rho}$ corresponds to the mean density of the universe,
$P(k)_{\text{lin}}$ is the linear theory matter power spectrum, $n(M)$ is the
halo mass function, and $b(M)$ is the linear halo bias. The term
$|\tilde{u}(k|M)|$ is the normalized Fourier transform of the internal density
profile of a halo of mass $M$, such that $\tilde{u}(k\rightarrow
0,M)\rightarrow 1$. The above HM building blocks are intrinsically redshift-
dependent functions, which, in principle, allows one to obtain HM prediction
at any redshift for which the integrands are well-defined.

All the components of the HM can be varied independently from each other, and
each specific choice of fitting functions, formulae, or tabulated data creates
a unique realization. Thus, HM is a general framework under which one can
create many different families of PS models. Motivated by literature and our
own studies for each of our cosmological models (i.e.
$\mathrm{\Lambda\text{CDM}}$, and all MG variants), we find an optimal
combination of analytic formulae and fitting functions to describe the input
properties of halo mass function, halo bias, and halo concentrations. Below we
provide a more detailed description of the particular choices we make. For a
quick summary and look-up, we refer the reader to Table 1 which contains a
concise list and references of all the fitting functions for the halo
properties used in this work, and for each model.

Table 1:  Compilation of the fitting functions used in this work for the halo properties in HM build-up, for both $\Lambda$CDM and MG models. Halo properties | Fitting functions | Notes
---|---|---
Halo mass function, HMF |  |
$\Lambda$CDM : Watson et al. 2013 Watson  _et al._ (2013) | $f(\sigma)_{\Lambda\text{CDM}}=A\left[\left(\frac{\beta}{\sigma}\right)^{\alpha}+1\right]e^{-\gamma/\sigma^{2}}$ |  For MG(=$f(R)$, nDGP),
|  A = 0.282, $\alpha$ = 2.163, $\beta$ = 1.406 and $\gamma$ = 1.210. | $f(\sigma)_{\text{MG}}=\Delta_{\text{MG}}\times f(\sigma)_{\Lambda\text{CDM}}$
$f(R)$: Gupta et al. 2022 Gupta  _et al._ (2022) | $\Delta_{f(R)}=1+a\exp\left[-\frac{(X-b)^{2}}{c^{2}}\right]$ | Additional cut-off expression at
| For F5: a = 0.230, b = 0.100 and c = 0.360 |  low-mass scales for $f(R)$ (Eq. 9)
| For F6: a = 0.152, b = -0.583 and c = 0.375 |
| $X\equiv\ln(\sigma^{-1})$ |
nDGP: Gupta et al. 2022 Gupta  _et al._ (2022) | $\Delta_{\text{nDGP}}=p+q\arctan{(s\,X+r)}$ |
|  |  $\Xi(z)$: nDGP force enhancement
| For N1: p = 1.35, q = 0.258, r = 5.12, s = 4.05 |  w.r.t. GR Koyama and Silva (2007).
| For F6: p = 1.06, q = 0.0470, r = 11.8, s = 4.19 |
|  $X\equiv\ln(\widetilde{\sigma}^{-1})$, $\widetilde{\sigma}=\sigma/\Xi(z)$ |
Linear halo bias, $b(M)$ |  |
All models: Sheth et al. 2001 Sheth  _et al._ (2001) | $b(M)=\frac{1}{\sqrt{a}\delta_{c}(z)}(\sqrt{a}(a\nu^{2})+\sqrt{a}b(a\nu^{2})^{1-c}$ | This expression has been proposed
| $-\frac{(a\nu^{2})^{c}}{(a\nu^{2})^{c}+b(1-c)(1-c/2)})$ |  for $\Lambda$CDM. We extrapolated
|  $a=0.707,b=0.5$ and $c=0.6$. | the relation for MG.
Concentration-mass relation, $c(M)$ |  |
$\Lambda$CDM: Ludlow et al. 2016 Ludlow  _et al._ (2016) | $c(\nu)_{\Lambda\text{CDM}}=c_{0}\left(\frac{\nu}{\nu_{0}}\right)^{-\gamma_{1}}\left[1+\left(\frac{\nu}{\nu_{0}}\right)^{1/\beta}\right]^{-\beta(\gamma_{2}-\gamma_{1})}$ |  For MG(=$f(R)$, nDGP),
| $c_{0}=3.395\times(1+z)^{-0.215}$ | $c(M)_{\text{MG}}=\Delta_{c(M),\text{MG}}\times c(M)_{\Lambda\text{CDM}}$
| $\beta=0.307\times(1+z)^{0.540}$ |
| $\gamma_{1}=0.628\times(1+z)^{-0.047}$ |
| $\gamma_{2}=0.317\times(1+z)^{-0.893}$ |
| $\nu_{0}=(4.135-0.564a^{-1}-0.210a^{-2}$ |
| $+0.0557a^{-3}-0.00348a^{-4})\times D(z)^{-1}$ |
$f(R)$: Mitchell et al. 2019 Mitchell  _et al._ (2019) | $y(x)=\frac{1}{2}\left(\frac{\lambda}{\omega_{s}}\phi(x^{\prime})\left[1+erf\left(\frac{\alpha x^{\prime}}{\sqrt{2}}\right)\right]+\gamma\right)(1-tanh(\omega_{t}[x+\xi_{t}]))$ |
| $y=log_{10}(\Delta_{c(M),f(R)})$ |
| $x^{\prime}=(x-\xi_{s})/\omega_{s}$ |
| $x=log_{10}(M_{500}/10^{p_{2}})$ |
|  $p_{2}=1.5\log_{10}\left[\frac{\bar{f_{R}(z)}}{1+z}\right]+21.64$ Mitchell  _et al._ (2018) |
| $\lambda=0.55\pm 0.18$ |
| $\xi_{s}=-0.27\pm 0.09$ |  For M $\leq 10^{12}M_{\odot}/h$,
| $\omega_{s}=1.7\pm 0.4$ | $c(M)_{f(R),\text{nDGP}}=c(M)_{\Lambda\text{CDM}}$
| $\alpha=-6.5\pm 2.4$ |
| $\gamma=-0.07\pm 0.04$ |
| $\omega_{t}=1.3\pm 1.0$ |
| $\xi_{t}=0.1\pm 0.3$ |
nDGP: Mitchell et al. 2021 Mitchell  _et al._ (2021) | $\Delta_{c(M),\text{nDGP}}=[A-Blog_{10}(M_{200}M_{\odot}h^{-1})](H_{0}r_{c})^{-0.71\pm 0.05}+1$ |
| $A=(0.35\pm 0.01)(H_{0}r_{c})^{-0.71\pm 0.05}$ |
| $B=(0.0302\pm 0.0008)(H_{0}r_{c})^{-0.71\pm 0.05}$ |

###  III.1 Halo mass function

The halo mass function (HMF), $n(M)$, quantifies the number of halos per unit
mass per unit comoving volume. The most commonly adopted theoretical
formulation of the HMF is via the Extended Press–Schechter (EPS) formalism
(Press and Schechter, 1974; Bond  _et al._ , 1991), in which HMF is given by:

$n(M)\equiv\frac{dn}{dM}=\frac{\rho}{M^{2}}F(\sigma)\left|\frac{d\ln\sigma}{d\ln
M}\right|$ (5)

The halo multiplicity function, $F(\sigma)=\nu F(\nu)$ denotes the fraction of
matter collapsed into halos, in a logarithmic bin around the peak height,
$\nu=\delta_{c}(z)/\sigma(M,z)$. Here, $\delta_{c}(z)$ is the spherical
collapse density threshold, and $\sigma(M,z)$ is the linear variance in the
density fluctuation field smoothed using a top-hat filter. This scaling
relation has been modeled extensively in the literature and it has been shown
to be approximately universal across redshifts for
$\mathrm{\Lambda\text{CDM}}$ (Sheth  _et al._ , 2001; Jenkins  _et al._ ,
2001; Warren  _et al._ , 2006; Watson  _et al._ , 2013; Despali  _et al._ ,
2016). In our earlier work (Gupta  _et al._ , 2022), we have shown that after
simple re-scaling, the $F(\sigma){}$ in both $f(R)$ and nDGP also exhibits a
similar degree of universality as in the $\mathrm{\Lambda\text{CDM}}$-case.

Following our previous study, we will model MG HMF as a fractional deviation,
$\Delta_{\text{MG}}$ from the $\mathrm{\Lambda\text{CDM}}$ fiducial baseline.
We have shown that such an approach allows for achieving quite a good accuracy
($5-10\%$), which also holds for different background cosmologies. However, to
obtain such precision, a careful choice of the baseline
$\mathrm{\Lambda\text{CDM}}$ HMF model is paramount.

Thus, for our baseline $\mathrm{\Lambda\text{CDM}}$, we tested various HMF
models in the literature (e.g. Sheth  _et al._ (2001); Jenkins  _et al._
(2001); Warren  _et al._ (2006); Tinker  _et al._ (2008); Angulo  _et al._
(2012); Watson  _et al._ (2013); Despali  _et al._ (2016)), as these functions
can in principle be extrapolated to desired halo mass ranges. We found that
the fitting function proposed in Watson et al. 2013 (Watson  _et al._ , 2013,
hereafter W13) proved to be optimal for HM power-spectrum forecasting.
Therefore, we used W13 for our $\mathrm{\Lambda\text{CDM}}$ HMF computations.

For completeness, we now recall the essential steps of Ref. Gupta  _et al._
(2022). Here the target MG HMF is modeled as:

$n(\sigma_{M})_{\text{MG}}=\Delta_{\text{MG}}(\sigma_{M})\cdot
n(\sigma_{M})_{\Lambda\text{CDM}}\,\,,$ (6)

where $\sigma_{M}\equiv\sigma(M)$ is simply the linear mass variance at the
Lagrangian top-hat halo mass scale, $M$.

For $f(R)$ gravity models, the fractional deviation fit is expressed as:

$\Delta_{\text{MG}}\equiv\Delta_{f(R)}=1+a\exp\left[-\frac{(X-b)^{2}}{c^{2}}\right],$
(7)

$X\equiv\ln(\sigma^{-1})$. Here, $(a,b,c)$ are parameters of the fit that were
calibrated using simulations. They depend on the variant of $f(R)$ gravity
model under consideration. See Table 1 for the specific values that we use in
this work.

For nDGP gravity models:

$\Delta_{\text{MG}}\equiv\Delta_{\text{nDGP}}=p+q\arctan{(s\,X+r)}.$ (8)

Here, $X$ is the re-scaled mass density variance,
$X\equiv\ln(\widetilde{\sigma}^{-1})$, $\widetilde{\sigma}=\sigma/\Xi(z)$.
Again, $(p,q,r,s)$ are the parameters of the fit, whose values are determined
by the variant of the nDGP gravity model.

The resolution of our simulations allowed us to probe only intermediate- and
large-mass halos to compute the HMF. In this mass regime, HMF in MG increases
w.r.t. $\mathrm{\Lambda\text{CDM}}$, as small-mass halos accrete matter and
merge faster to form larger structures. However, this enhanced structure
formation at large halo mass-end is happening at the expense of the abundance
of smaller halos used up in this process (see e.g. Hellwing and Juszkiewicz,
2009; Lam and Li, 2012; Falck  _et al._ , 2015; Gupta  _et al._ , 2022). Thus,
we can expect that there should be a simultaneous decrease in the number of
small-mass halos in the MG models when compared to
$\mathrm{\Lambda\text{CDM}}$.

Equation (8) for nDGP allows the possibility of $\Delta_{\text{nDGP}}<1$ for
small mass halos. However, our fit for $\Delta_{f(R)}$ is never below 1. To
admit for low-mass halo deficiency also in the $f(R)$, we impose an artificial
decrease in $f(R)$ HMF for M $<10^{11}$ M${}_{\odot}/h$, when compared to
$\mathrm{\Lambda\text{CDM}}$ results. For low halo masses, we assume that
$\Delta_{f(R)}$ is a linear function of $\ln(\sigma^{-1})$, and is given by:

$\Delta_{f(R)}\rightarrow(m\ln(\sigma^{-1})+n)\times\Delta_{f(R)}$ (9)

We tested for different combinations of the $(m,n)$ parameters values. The
combination $(m,n)=(0.06,0.99)$ turned-out to be optimal for our both $f(R)$
variants. Thus we use these values in this work. A note of caution is in place
here. There is no clear physical justification for our particular choice of
both $m$ and $n$, other than that they are providing optimal HM power spectrum
predictions. An interested reader can play around and search for a different
choice of $(m,n)$. However, the overall impact of the particular $(m,n)$
choice on the resulting HM remains small.

###  III.2 Halo bias

The relation between the clustering amplitude of the underlying DM density
field and halos is quantified in terms of the linear halo bias relation,
$\delta_{h}(M)=b(M)\delta$. In the context of power spectra, it is convenient
to consider the following Fourier-space estimator of the halo bias:

$\hat{b}(k,M)=\frac{P_{hm}(k,M)}{P(k)}$ (10)

Here, $P_{hm}(k,M)$ is the halo-matter cross power spectrum, and $P(k)$ is the
matter power spectrum. One can find an optimal value of the linear bias by
taking a limit, or an average of this estimator at the smallest possible
$k$’s. We consider such a power-spectrum-based bias estimator to use results
from elephant suite for testing and finding optimal analytic bias formula for
HM.

For this purpose, we tested various $b(M)$ fitting functions for
$\mathrm{\Lambda\text{CDM}}$ (Seljak and Warren, 2004; Tinker  _et al._ ,
2010; Comparat  _et al._ , 2017). Sheth et al. 2001 $b(M)$ (Sheth  _et al._ ,
2001, hereafter S01) gave the best match to the simulations. Thus, this will
be our choice for the $b(M)$ computations in this work.

In the Appendix A, we show the performance of S01, both in capturing the ratio
of MG $b(M)$ versus $\mathrm{\Lambda\text{CDM}}$ (Fig. 5), and the absolute
$b(M)$ relation (Fig. 6) . We find that S01 gives reasonable predictions in
both cases. Given that $b(M)$ impacts only the two-halo term, which by
construction matches the $P(k)_{\text{lin}}$ on large scales, the choice of
$b(M)$ does not impact the HM results to a great extent.

###  III.3 Halo density profile: concentration-mass relation

The scale-free nature of structure formation in CDM scenarios results in self-
similar density profiles for individual DM halos, which was first pointed out
by Navarro, Frenk, and White in (Navarro  _et al._ , 1996, hereafter NFW). As
a result, DM density profiles are re-scaled by a characteristic central
density, $\rho_{s}$, and radial scale, $r_{s}$, (or mass $M$ and concentration
$c(M)$, respectively). The $c(M)$ relation is defined as the ratio of the
virial radius, $R_{v}$ of the halo to $r_{s}$, and determines the density
profile of NFW halos.

To obtain relatively unbiased and good-quality NFW fits, the simulated halos
need to be well-resolved. The convergence of the halo density profile depends
on the simulation’s force and mass resolution. Thus $c(M)$ can be reliably
estimated only for a limited halo mass range, usually for halo with masses
corresponding to at least a few $\times$ 103 particles (see e.g. Power  _et
al._ , 2003). The resolution of the elephant suite allow only for probing the
$c(M\geq 10^{13}$M${}_{\odot}/h$). Because of this, we need to resort to the
fitting functions for $c(M)$ here.

We use relations proposed in (Mitchell  _et al._ , 2019) and (Mitchell  _et
al._ , 2021) to compute the $c(M)$ relation in $f(R)$ and nDGP gravity models,
respectively. In these works, direct NFW fitting was used to compute the halo
density profiles, and functional forms were derived for the ratio
$c(M)_{\text{MG}(=f(R),\text{nDGP})}/c(M)_{\Lambda\text{CDM}}$ (refer to Table
1 for explicit expressions). The MG $c(M)$ can be therefore obtained as a
product of this ratio times the concentration-mass relation for
$\mathrm{\Lambda\text{CDM}}$, for which we use the form proposed in (Ludlow
_et al._ , 2016). Considering the ratio instead of absolute MG $c(M)$ would
eliminate the leading-order systematic uncertainties coming from the
background cosmology.

The authors in Mitchell  _et al._ (2019) proposed functional form for
$\log(c(M)_{f(R)}/c(M)_{\Lambda\text{CDM}})$. When expressed as a function of
$M_{500}/10^{p_{2}}$, this ratio is independent of the background scalar field
and $z$. The parameter $p_{2}$ defined in Mitchell  _et al._ (2018),
encapsulates these dependencies, and in turn allows different variants of
$f(R)$ gravity model to be studied in a unified way.

For the case of the nDGP gravity model, in ref. (Mitchell  _et al._ , 2021),
the ratio $c(M)_{\text{nDGP}}/c(M)_{\Lambda\text{CDM}}$ is fitted as a
decreasing function of $M_{200}$. This fitting also captures the $z$
dependence, hence making the ratio only dependent on the nDGP parameter, rcH0.

The halo mass range probed in both Mitchell  _et al._ (2019, 2021) is confined
to $\geq$ 1012 M${}_{\odot}/h$. Therefore, we restrict the use of their
fitting functions to the calibrated mass range, and artificially impose
$c(M)_{\text{$f(R)$, nDGP}}=c(M)_{\Lambda\text{CDM}}$ for $M<10^{12}$
M${}_{\odot}/h$.

##  IV Results

In this section, we combine all the HM components to give an analytical
prediction for matter overdensity PS. As a reference case to gauge our results
against, we always take the PS from elephant simulations. At large scales, the
linear perturbation theory gives accurate and reliable predictions both for
$\mathrm{\Lambda\text{CDM}}$ and MG PS. Hence, we focus here only on the
scales corresponding to mildly and fully non-linear regimes. In practice, we
will be interested in the performance of our models for $k\geq 0.1$
$h/\,\mbox{Mpc}$.

###  IV.1 Halo model predictions for modified gravity

We start by testing the standard set-up for HM, which aims to yield a
theoretical prediction for the PS amplitude in a given cosmology. For
$\mathrm{\Lambda\text{CDM}}$ alone, this approach has at best limited
accuracy, since the classical HM fails to accurately capture PS already in the
mildly non-linear regime, i.e.
$k\stackrel{{\scriptstyle>}}{{{}_{\sim}}}0.2-0.3$ $h/\,\mbox{Mpc}$ Smith  _et
al._ (2003); Tinker  _et al._ (2005); Mead  _et al._ (2015); Mead and Verde
(2021). Thus, we do not expect that it will perform better in MG cosmologies,
which have even richer phenomenology. However, it is still an illustrative
exercise, since we will be using this basic HM set-up to obtain much more
accurate PS predictions for MG.

Using the inputs of the HMF, $b(M)$, and $c(M)$ relation in their MG versions
discussed in the previous section, we compute the resultant power spectra for
a number of redshifts. For this, we employ Eqs. (2)-(4), integrating from
M${}_{\text{min}}$ = 1 M${}_{\odot}/h$ to M${}_{\text{max}}$ = 1016
M${}_{\odot}/h$. We choose a sufficiently broad halo mass range so as to
account for the maximum possible halo masses that still have an impact on the
resulting PS.

For the integral in Eq. (4) to approach unity at large scales, the bias needs
to attain unity when integrated over all the halo masses, i.e.

$\frac{1}{\bar{\rho}}\int_{0}^{\infty}b(M)n(M)dM=1$ (11)

In practice, this integral yields a value below unity, even when the
integration is taken over the maximum possible halo mass range. Changing the
high mass limit for the integration does not impact the results to a great
extent, because on these scales, halos become exponentially rare which makes
their contribution to the total power negligible. On the other hand, we expect
a significant contribution from the low-mass regime. However, owing to
resolution limits, the properties of low mass halos cannot be properly
calibrated using simulations.

Therefore, to add the contribution of the low-mass halos to HM computations,
we use the correction proposed in Schmidt (2016); Mead  _et al._ (2020). This
correction adds the contribution of the missing halos to the two-halo term, in
order to recover $P(k)_{\text{lin}}$ at large scales. The correction term is
simply yielded by:

$A=1-\frac{1}{\bar{\rho}}\int_{M_{\text{min}}}^{M_{\text{max}}}b(M)n(M)dM,$
(12)

and it is used as an additive component in the two-halo term:

$C=\frac{A\tilde{u}(k|\text{M}_{\text{min}})}{\text{M}_{\text{min}}}$ (13)

Here, $\tilde{u}(k|\text{M}_{\text{min}})$ is the normalized Fourier transform
of the density profile for the lowest resolved mass M${}_{\text{min}}$.
Equation (3) is then modified and the resultant two-halo term is given by:

$P(k)_{2h}=P(k)_{\text{lin}}(I_{m}+C)^{2}$ (14)

One could instead replace the $P(k)_{2h}$ term with $P(k)_{\text{lin}}$, as
the former differs from the latter only for $k\geq 1$ $h/\,\mbox{Mpc}$, where
already $P(k)_{1h}$ takes over as the dominant contributor. However, for
completeness, we use the full above expression for the two-halo term.

The results of such direct HM computations for our MG models are illustrated
in Fig. 1, where we compare $P(k)_{\text{HM}}$ (solid lines), as well as
linear theory $P(k)_{\text{lin}}$ (dotted lines), with the elephant
simulations for all our models at $z=0$. The shaded region corresponds to the
uncertainty in the elephant results, which is the inverse of the square root
of the number of statistically independent modes contributing to each $k$-bin,
and the horizontal dashed lines correspond to the 5% accuracy regime. The
performance of HM in these MG models is similar to the
$\mathrm{\Lambda\text{CDM}}$ results and is not much better than the actual
linear theory. With respect to the simulation prediction, $P(k)_{\text{HM}}$
gives better than 5% accuracy for $k\leq 0.2-0.3$ $h/\,\mbox{Mpc}$, and stays
within 10% for $k\leq 0.4-0.5$ $h/\,\mbox{Mpc}$. An interesting exception is
the F5 $f(R)$ variant, where better than 2% accuracy is kept all the way to
$k\sim 0.2$ $h/\,\mbox{Mpc}$.

In all the models, we also encounter an under-prediction w.r.t. the simulation
results for $k\approx 0.5$ $h/\,\mbox{Mpc}$. This is a well-known problem of
the HM formalism in $\mathrm{\Lambda\text{CDM}}$ (Mead  _et al._ , 2015; Mead
and Verde, 2021), and further propagates to the MG scenarios (also seen for
Galileon models in (Barreira  _et al._ , 2014)). Similar behavior is observed
also for other redshifts that our simulations probe, but we do not show them
here for brevity.

These results clearly show that the standard HM cannot be used in MG,
similarly as in $\mathrm{\Lambda\text{CDM}}$, if one is interested in
attaining the accuracy that we expect from future LSS surveys. However,
noticing the above, the positive result here is that HM can be actually
employed to yield predictions for MG power spectra with the same-level
accuracy as for $\mathrm{\Lambda\text{CDM}}$. This is a somewhat surprising
result because the standard HM does not include any room for extra MG physics
(like the fifth-force and screening). Yet it seems that self-consistent
modifications of HMF, $b(M)$, and $c(M)$ are enough to obtain the usual
$\mathrm{\Lambda\text{CDM}}$ HM-level predictions also for different MG
cosmologies. This is very encouraging, and as we show below this can be used
as a strong advantage to build an even better and more accurate PS model for
MG.

Figure 1: Comparison of the power spectrum from halo-model, $P(k)_{\text{HM}}$
(solid lines) and linear theory, $P(k)_{\text{lin}}$ (dotted lines) with
elephant simulations $P(k)_{\text{sim}}$, for $\mathrm{\Lambda\text{CDM}}$ and
the variants of MG models, as in the legend, for $z=0$. The horizontal grey
dashed lines correspond to 5% accuracy regime.

###  IV.2 An improved model for MG power spectrum

Figure 2: The ratio $\Upsilon(k)$ $\equiv$
$P(k)_{\text{MG}}/P(k)_{\Lambda\text{CDM}}$ obtained from linear theory
(dotted lines), elephant simulation (solid lines) and the halo model (dashed
lines), at $z=0$ (top panels) and $z=0.5$ (bottom panels). The left panels
correspond to $f(R)$ gravity variants: F5 and F6, and the right panels
correspond to nDGP gravity variants: N1 and N5. Shaded regions are the
propagated errors from $P(k)_{\text{sim}}$.

Figure 3: Top panel: Matter power spectra obtained from our new approach
($P(k)_{\text{MG}}$ from Eq. (16)) for all the MG variants considered in this
work, at redshifts as indicated in the legends. Bottom panel: Comparison of
$P(k)_{\text{MG}}$, derived with our new method, with $P(k)_{\text{sim}}$. The
shaded region in all the plots corresponds to the uncertainty in the
$P(k)_{\text{sim}}$, and the horizontal dotted lines shows 5% accuracy regime.

In the previous section, we have shown that when HM is applied to model the PS
amplitude, it offers limited accuracy, and is comparable to what can be
achieved for the standard $\mathrm{\Lambda\text{CDM}}$. In this section, we
will demonstrate that we can build a much more accurate PS model for MG. This
can be realized when we apply HM to estimate the fractional departure from the
$\mathrm{\Lambda\text{CDM}}$ baseline, rather than trying to predict the
absolute amplitude of PS alone.

Our starting point will be the ratio of the power spectra:

$\Upsilon(k)\equiv P(k)_{\text{MG}}/P(k)_{\Lambda\text{CDM}}\,.$ (15)

By modeling this ratio, rather than the MG PS itself, we can benefit from a
number of properties, namely: (i) the dependence on the background
cosmological parameters (such as $\Omega_{m}$, $H_{0}$, or $\sigma_{8}$)
should cancel out from the ratio to the leading order; and (ii) the scale of
significant departure from $\mathrm{\Lambda\text{CDM}}$ (i.e. from
$\Upsilon=1$) is naturally determined in terms of the
$\mathrm{\Lambda\text{CDM}}$ baseline, rather than some arbitrary non-linear
amplitude or scale.

In Fig. 2, we compare the ratios $\Upsilon(k)$ estimated from elephant
simulations (solid line), linear theory (dotted line), and HM (dashed line),
for both $f(R)$ (left panels), and nDGP (right panels), at $z=0$ (top panel)
and $z=0.5$ (bottom panel). Naturally, both the simulation and HM results for
$\Upsilon(k)$ are expected to match the linear theory prediction at large
scales for both models. As we approach smaller scales, the departure from
linear predictions increases (namely, linear theory runs away for $f(R)$ and
stays constant for nDGP), and trends peculiar to each model emerge. This is a
well-known result, which highlights the fact that these family of MG models
usually exhibits an increased degree of non-linearity of the density field,
owing to both the fifth force and their respective screening mechanisms
(Koyama  _et al._ , 2009; Sotiriou and Faraoni, 2010; Joyce  _et al._ , 2015;
Puchwein  _et al._ , 2013; Oyaizu  _et al._ , 2008; Schmidt, 2009; Li  _et
al._ , 2013b). For $f(R)$ gravity models, PS approaches
$\mathrm{\Lambda\text{CDM}}$ on the large-scales, and we see a monotonic
increase in the ratio with $k$ (although slower than what the linear theory
would predict). Whereas, for nDGP, $\Upsilon(k)$ enhancement is maximum at the
intermediate scales, and this enhancement decreases for large $k$.

A crucial observation from our study is that HM prediction agrees
qualitatively with the simulations. We note that the agreement is far from
perfect, especially around the peak-like features, but the HM captures the
essential shape and scales of the PS ratios.

As mentioned in the previous section, one perennial problem with the HM has
been the ‘transition’ region, where both two- and one-halo terms have a
similar magnitude, and both contribute equivalently to the predicted signal.
In general, the HM under-predicts the strength of clustering in this region,
with the exact amount depending on redshift and cosmology (Mead  _et al._ ,
2015). We also highlight a similar problem with the HM-based MG predictions in
Fig. 1 at $k\approx 0.5$ $h/\,\mbox{Mpc}$. These scales are also called the
‘quasi-linear’ regime because the evolution of perturbations at these scales
is not exactly governed by linear perturbation theory.

For standard $\mathrm{\Lambda\text{CDM}}$, the inaccuracies of the HM in this
transition regime are addressed by devising empirical fitting functions. One
of the earliest, yet successful examples was halofit Smith  _et al._ (2003),
which is motivated by the principles of HM, and calibrated using N-body
simulations. It was later improved, in particular by (Takahashi  _et al._ ,
2012) who updated its fitting functions from higher resolution simulations and
ameliorated the modeling for dark energy cosmologies. Methods and
prescriptions to predict the non-linear PS in $\mathrm{\Lambda\text{CDM}}$ are
numerous, but in this work we will use halofit as it is sufficiently accurate
for our purposes.

Having seen that the ratio $\Upsilon(k)$ between HM-derived PS for MG and
$\mathrm{\Lambda\text{CDM}}$ correctly captures the simulation trends, we
propose to use it to obtain the fully non-linear PS in MG. This is obtained by
multiplying $\Upsilon_{\mathrm{HM}}$ with an accurate model for the
$\mathrm{\Lambda\text{CDM}}$ baseline $P(k)$. Therefore, we characterize the
beyond-$\mathrm{\Lambda\text{CDM}}$ PS ($P(k)_{\text{MG}}$) as:

$P(k)_{\text{MG}}=\Upsilon(k)_{\text{HM}}\times P(k)_{\Lambda\text{CDM}}$ (16)

In this prescription, $P(k)_{\Lambda\text{CDM}}$ is the non-linear
$\mathrm{\Lambda\text{CDM}}$ PS, for which we take the halofit predictions
using the parameters of a given background cosmology.

The results of applying our proposed methodology are illustrated in Fig. 3,
where we plot $P(k)_{\text{MG}}$ obtained using HM for $\Upsilon(k)$
multiplied by the halofit $\mathrm{\Lambda\text{CDM}}$-baseline. The top
panels present the power spectra directly: $P(k)_{\text{sim}}$ from elephant
(dots), and $P(k)_{\text{MG}}$ derived with Eq. (16) (solid lines). In the
bottom panels, we show departures of thus-obtained $P(k)_{\text{MG}}$ from
$P(k)_{\text{sim}}$ treated as reference. These new results, when compared
with the standard HM predictions from Fig. 1, clearly perform much better. The
standard HM reaches 5% accuracy only up to $k\leq 0.2-0.3$ $h/\,\mbox{Mpc}$.
Now, by using HM only for predicting $\Upsilon(k)$, and combining it with
halofit $\mathrm{\Lambda\text{CDM}}$-baseline, we improve the scale at which
modeling is accurate within 5% by an order of magnitude, reaching up to $k\leq
0.5-2.5$ $h/\,\mbox{Mpc}$ (depending on the model and redshift). We note that
the performance of $P(k)_{\text{MG}}$ generally worsens for higher redshifts,
but still remains significantly improved when compared to the standard HM.

Figure 4: Comparison of our $P(k)_{\text{MG}}$ modeling, with simulation
results from mg-cola Winther  _et al._ (2017); Naidoo  _et al._ (2022), for
two MG variants: N1 (left) and F5 (right). The redshifts are as indicated in
the legends. Dotted grey lines are the 5% accuracy regime. The shaded region
is the simulation error, which is the standard deviation obtained from five
realizations at each redshift.

More generally, the accuracy of $P(k)_{\text{MG}}$ will depend on the user
input of baseline $P(k)_{\Lambda\text{CDM}}$. As already mentioned, other
approaches are being developed to further improve the limited accuracy of
halofit, especially for models departing from the flat Planck-based
$\mathrm{\Lambda\text{CDM}}$. We tested one such alternative way of deriving
the non-linear $\mathrm{\Lambda\text{CDM}}$ PS, that goes into our
$P(k)_{\text{MG}}$ prediction (16): the so-called hmcode (Mead  _et al._ ,
2015, 2021). The results, detailed in Appendix B, indicate that both halofit
and hmcode give similar accuracy, however with different trends at different
scales and redshifts.

Given the fact that we have calibrated our MG HM with a limited-resolution
elephant simulation suite, it is encouraging that this allowed for already an
order-of-magnitude improvement of the scale at which we can obtain accurate PS
predictions. Obtaining accurate MG PS into the fully non-linear regime at
$k\geq 1$ $h/\,\mbox{Mpc}$ with so straightforward modifications to HM opens
up an avenue for even better PS predictions for the MG phenomenology. Even
better PS accuracy could be achieved using possible improvements to HM, better
informed about the clustering and properties of small halo mass regime in MG.

###  IV.3 Testing non-linear $P(k)_{\text{MG}}$ beyond ELEPHANT

In this subsection, we extend our work beyond the elephant simulations to test
the performance of our new approach. For this purpose, we consider different
N-body simulation runs for the F5 and N1 models, performed using mg-cola
Winther  _et al._ (2017), and described in Naidoo  _et al._ (2022). The most
important factors for us are that these runs have a different background
cosmology than elephant (see Section II), and were run using different codes.
Unlike standard $N$-body approach, these simulations employ the COLA method
Tassev  _et al._ (2013), that can straightforwardly trade accuracy at small-
scales in order to gain computational speed without sacrificing accuracy at
large scales. On one hand, this approach is much faster than the standard
$N$-body, but the price to pay is the approximations made, which do not allow
us to use these suite of simulations as the calibration data. Hence, we use
these simulations but only as a test-bed.

Here, the HM ingredients were calculated using the same methodology and setup
as above, described in Section III and summarized in Table 1. The main
difference with respect to Sec. IV.2 was that different background
cosmological parameters were used in the linear power spectra that go into the
particular ingredients of the HM build-up, namely HMF, $b(M)$ and $c(M)$.
Everything else, including the halo mass integration ranges for the HM
components, were the same as before.

Using the HM outputs and $\mathrm{\Lambda\text{CDM}}$ halofit predictions for
the background cosmology of this alternative simulation suite, we computed
$P(k)_{\text{MG}}$. A comparison of our results with the simulation
predictions is in Fig. 4, for both N1 (left plot) and F5 (right plot). Given
the small box size of these simulations ($L$ = 500 $\mbox{Mpc}/h$), we obtain
a discrepancy $>5\%$ with the simulation predictions on large scales, for
$k<0.1-0.2$ $h/\,\mbox{Mpc}$. Now, contrary to the elephant results, for N1
our new PS model performs better than before. However, for F5 the performance
of our approach decreases with increasing redshift. Overall, we see a similar
performance of the new $P(k)_{\text{MG}}$ in both the simulations that we
tested, with 5% accuracy from mildly non-linear to non-linear scales ($k\leq
0.5-2.5$ $h/\,\mbox{Mpc}$).

This test with a different simulation and cosmology reassures us that our new
approach is a valid technique to compute the non-linear PS in these MG
scenarios, and can be successfully extended to simulations and cosmologies
beyond our original data used for calibration and fitting.

##  V Discussion and conclusion

In this work, we combined Halo Model (HM) predictions with an accurate
$\mathrm{\Lambda\text{CDM}}$ baseline for building an analytical framework to
compute the non-linear power spectrum (PS) in modified gravity (MG) scenarios,
where structure formation differs from that in $\mathrm{\Lambda\text{CDM}}$.
For calibration and testing, we used the elephant suite – a set of $N$-body
simulations, which incorporates standard $\mathrm{\Lambda\text{CDM}}$ and two
MG models: Hu-Sawicki $f(R)$ and the normal branch of the Dvali-Gabadadze-
Porrati braneworld (nDGP). HM has been extensively studied for
$\mathrm{\Lambda\text{CDM}}$ (Seljak, 2000; Peacock and Smith, 2000), and we
further extended it to these MG cosmologies. This formalism is advantageous as
a quick and reliable tool to obtain predictions for statistics of density
fields well into the regimes, where linear and perturbation theory fails to
reproduce simulation results.

The HM framework requires the input of three main halo properties: halo mass
function (HMF), which quantifies the number density of halos; linear halo bias
$b(M)$, describing the relation between halos and the underlying DM density
field; and the concentration-mass relation $c(M)$, which describes the
internal distribution of mass in halos. For the HM framework, we needed to
compute these quantities over large range of halo masses, that go much beyond
the range of our simulations. As a result, we relied on fitting functions for
the halo properties in these MG scenarios (Table 1).

Using these three inputs, we obtained the HM-based predictions,
$P(k)_{\text{HM}}$ for these two MG models. We showed that $P(k)_{\text{HM}}$
is within 5-15% of the simulation results across the $k$-ranges, from $k=0.01$
to $k=1$ $h/\,\mbox{Mpc}$. However, MG signatures from these models that
quantify deviations from GR are typically in itself a factor of a few dozen
per cent. Hence, we cannot use HM predictions in its standard form to
complement the expected accuracy from future LSS surveys in order to detect
these MG signals. Additionally, similar to the case of
$\mathrm{\Lambda\text{CDM}}$, HM also faces the consistent problem of under-
prediction of power in the transition regime for both $f(R)$ and nDGP. These
scales correspond to $k\approx 0.5$ $h/\,\mbox{Mpc}$.

To get a better PS model, we further investigated, using HM, the relative
ratio $\Upsilon(k)=P(k)_{\text{MG}}/P(k)_{\Lambda\text{CDM}}$, instead of
employing the absolute PS amplitudes alone. From this, we obtained new
analytical PS by taking a product of $\Upsilon(k)_{\text{HM}}$, with the non-
linear prediction for $\mathrm{\Lambda\text{CDM}}$, $P(k)_{\Lambda\text{CDM}}$
(Eq. 16). For the latter, we used halofit (Smith  _et al._ , 2003; Takahashi
_et al._ , 2012), as it has been a successful approach for
$\mathrm{\Lambda\text{CDM}}$ to circumvent the HM under-prediction of PS in
the intermediate scales, and is widely used to analytically compute non-linear
$\mathrm{\Lambda\text{CDM}}$ PS. One could use other approaches for the input
non-linear $P(k)_{\Lambda\text{CDM}}$ to multiply our ratio
$\Upsilon_{\text{HM}}(k)$ with, and we tested one of them: hmcode (Mead  _et
al._ , 2015, 2021). It gave comparable accuracy as halofit for our cosmology.

Using this approach, we significantly improved the accuracy of PS modeling
compared to standard HM prediction. For HM, we obtained results within 5%
accuracy with the simulation predictions for $k\leq 0.2-0.3$ $h/\,\mbox{Mpc}$.
With our new approach, we now matched simulations within this accuracy for $k$
between $0.5-2.5$ $h/\,\mbox{Mpc}$, with the performance of the method
depending on the MG model and redshift. The $k$-range probed in this work
corresponds to the mildly non-linear and the fully non-linear regime: a range
of scales crucial to constrain modern era cosmological observables. The
sensitivity of these observables to changes in the matter PS will be very
important for making powerful observational cosmological tests of the theory
of gravity, or dark energy.

The main advantage of our approach over using simulations is that it is
computationally inexpensive. The two main inputs, HM and halofit (or e.g.
hmcode) can be flexibly applied to a different background cosmologies,
whereas, in simulations, we need to perform a new run for each new set of
parameters. HM also gives the flexibility of employing different combinations
of HMF, $b(M)$ and $c(M)$ that is best suited to probe a particular cosmology,
scale, halo mass range, or redshift.

To test the limits and accuracy of our approach, we applied Eq. 16 to another
suite of MG N-body simulations, run with mg-cola Winther  _et al._ (2017);
Naidoo  _et al._ (2022). Using the same fitting functions as described for the
elephant simulations, we computed halo properties for the mg-cola cosmology,
and in turn $P(k)_{\text{HM}}$, $\Upsilon(k)_{\text{HM}}$ and
$P(k)_{\text{MG}}$. We compared $P(k)_{\text{MG}}$ with the
$P(k)_{\text{sim}}$ results, and obtained similar accuracy as with the
original data of elephant. For N1, $P(k)_{\text{MG}}$ gives consistency with
simulation results within 5% for $k$ between $0.5-2.5$ $h/\,\mbox{Mpc}$.
However, for the case of F5, overall performance of our approach decreases
with $z$.

We need to appreciate that another promising solution for analytical modeling
of the MG PS is via the fast and reliable emulation techniques (e.g Heitmann
_et al._ , 2006; Habib  _et al._ , 2007; Angulo  _et al._ , 2021).For MG
models, emulators have been proposed in e.g. (Ramachandra  _et al._ , 2021;
Winther  _et al._ , 2019; Arnold  _et al._ , 2022; Harnois-Déraps  _et al._ ,
2022; Ruan  _et al._ , 2023). This approach is sophisticated and promising,
however, is still in its infancy, and has limitations. For instance,
predictions from emulators are confined to the parameter space defined in the
starting base grid of the calibrating simulations. In addition, new extensions
in emulators (e.g. new degrees of freedom, or additional screening mechanisms
in MG models) often requires one to substantially adapt the base grid of
simulations used to build the emulator, which can in itself be computationally
expensive.

On the other hand, HM potentially provides a simple, physically-motivated
semi-analytical picture of the clustering of matter. We showed that HM, in its
standard form, can be qualitatively used to predict estimates for MG
signatures in cosmological observables which relate to matter perturbations.
Furthermore, using HM for modeling the PS ratio $\Upsilon(k)$, and combining
it with a high-quality baseline $\mathrm{\Lambda\text{CDM}}$ predictions
yields significantly better results. This method is advantageous as contrary
to MG scenarios, we have much tighter constraints on
$\mathrm{\Lambda\text{CDM}}$ physics, and the field of modeling
$\mathrm{\Lambda\text{CDM}}$ PS is much more sophisticated and advanced (Mead
_et al._ , 2015; Takahashi  _et al._ , 2012; Smith  _et al._ , 2003; Lawrence
_et al._ , 2010). As a result, more precise $\mathrm{\Lambda\text{CDM}}$
results will provide MG PS with similarly improved performance. Here we
present our results by incorporating the halofit and hmcode predictions for
$\mathrm{\Lambda\text{CDM}}$. These results in themselves give a percent level
of accuracy in both quasi-linear and non-linear regimes.

In order to further improve HM modeling in the MG variants studied here, we
need to probe deeper into the non-linear scales. For this, the behavior of
halo density profiles and HMF in both $f(R)$ and nDGP at low halo masses
requires deeper investigation, as the full effect of the respective screening
mechanisms comes to play in the non-linear regime of gravitational collapse.
As mentioned above, the accuracy of $c(M)$ fitting functions for both $f(R)$
and nDGP has not been tested for M${}_{\text{halo}}<10^{12}$ M${}_{\odot}/h$.
Additionally, we also extrapolated our earlier HMF fits for these MG models
Gupta  _et al._ (2022) to small halo mass scales, which are not resolved by
our $N$-body simulations (the limit being M${}_{\text{halo}}\lesssim 8\times
10^{12}$ M${}_{\odot}/h$). Both $c(M)$ and HMF are important ingredients in
modeling the one-halo term, which is the dominant non-linear contributor in
the HM approach. Such a study will require a completely new set of high-
resolution MG N-body simulations, and we plan it as a future project.

We also note that in this work, we focus on modeling only the dark matter PS.
At our scales of interest ($k$ between $0.1-2.5$ $h/\,\mbox{Mpc}$), PS is not
significantly influenced by baryons, as baryonic suppression in PS is of the
order of a few percent for $k<1-5$ $h/\,\mbox{Mpc}$ Acuto  _et al._ (2021);
Arnold  _et al._ (2019); Chisari  _et al._ (2019). However, Mead  _et al._
(2020, 2021); Acuto  _et al._ (2021) have shown that HM provides the
flexibility, which allows it to add additional parameters that can incorporate
baryonic effects from hydrodynamical simulation. Accounting for such effects
in our MG PS modeling is a significant endeavor, and is well beyond the scope
of this work.

The data used in this work is publicly available from our
website111https://data.cft.edu.pl/UPSILON_PK/UpsilonPk.tar.gz. We provide
$\Upsilon(k)$ for a wide range of $z$, from $z=0$ to $z=2$ for each MG model
considered in this work. A description of the data set is also enclosed in the
directory. Also, the data used to make the figures in this article are
available on request to the authors.

###### Acknowledgements.

We would like to thank Hans A. Winther for kindly providing us with mgcamb
version for specific forms of $\mu(a,k)$ and $\gamma(a,k)$ functions
implementing our $f(R)$ and nDGP models. The authors also thank Krishna Naidoo
for giving his suite of mg-cola simulations that we used as a testbed in our
work. We also acknowledge constructive discussions with Marika Asgari and
Alberto Acuto. Matter power spectrum calculations to compute linear halo bias
have been done using pylians222https://pylians3.readthedocs.io/en/master/.
This work is supported via the research project “VErTIGO” funded by the
National Science Center, Poland, under agreement no. 2018/30/E/ST9/00698. We
also acknowledge the support from the Polish National Science Center within
research projects no. 2018/31/G/ST9/03388, 2020/39/B/ST9/03494 (WAH & MB),
2020/38/E/ST9/00395 (MB), and the Polish Ministry of Science and Higher
Education (MNiSW) through grant DIR/WK/2018/12. This project also benefited
from numerical computations performed at the Interdisciplinary Centre for
Mathematical and Computational Modeling (ICM), the University of Warsaw under
grants no GA67-17 and GB79-7.

##  Appendix A Comparison of theoretical and simulation halo bias

Figure 5: Ratios of halo bias, $b(M)$, between MG and
$\mathrm{\Lambda\text{CDM}}$ for N1 (left plot), and F5 (right plot), across
range of redshifts as indicated in the legends. Solid lines are the analytical
results from Sheth et al. 2001 Sheth  _et al._ (2001), and the respective dots
of the same color are from simulations. Error bars illustrate the propagated
standard deviation from simulations.

Figure 6: Top panels: Linear halo bias, $b(M)$, as a function of halo mass,
$M_{200}$ for N1 (left column) and F5 (right column). The solid lines
correspond to theoretical Sheth et al. 2001 (Sheth  _et al._ , 2001, S01)
predictions, and the respective dots of the same color are the simulation
results obtained using Eq. 10. Error bars correspond to the standard deviation
across five realizations of the simulation box. Bottom panels: Ratio between
S01 and simulation linear halo bias predictions. Grey dotted lines are 5%
accuracy regimes.

Here we compare the simulation results for the linear halo bias, $b(M)$, with
the theoretical predictions from Sheth et al. 2001 (Sheth  _et al._ , 2001,
hereafter S01). The formula proposed by S01 is given by:

$\begin{split}b(M)=\frac{1}{\sqrt{a}\delta_{c}(z)}(\sqrt{a}(a\nu^{2})+\sqrt{a}b(a\nu^{2})^{1-c}\\\
-\frac{(a\nu^{2})^{c}}{(a\nu^{2})^{c}+b(1-c)(1-c/2)},\end{split}$ (17)

with the parameters $a=0.707,b=0.5$ and $c=0.6$.

To apply the S01 expression to our MG variants, we used $P(k)_{\text{lin}}$ to
compute $\sigma(M,z)$, and then $\nu=\delta_{c}(z)/\sigma(M,z)$, specific to
each MG model. For that, we used standard $\mathrm{\Lambda\text{CDM}}$
spherical collapse based $\delta_{c}$ values. We stay with the
$\mathrm{\Lambda\text{CDM}}$ $\delta_{c}$ baseline since we have found that
using slightly different values suggested for either $f(R)$ Kopp  _et al._
(2013), or for the nDGP model Schmidt  _et al._ (2010) impacts the final HM
results by less than a sub-percent.

The results of our substitution are shown in Fig. 5, where we plot the bias
ratios between MG and $\mathrm{\Lambda\text{CDM}}$, as a function of halo
mass, $M_{200}$. Here, we include the two models most departing from
$\mathrm{\Lambda\text{CDM}}$: N1 (left column) and F5 (right column). These
variants illustrate the most extreme behavior in $b(M)$ for the two MG models
we work with. Points illustrate simulation results (with error bars
corresponding to the propagated standard deviation from simulations). For
comparison, ratios of S01 predictions for MG and $\mathrm{\Lambda\text{CDM}}$
are also shown, but they are extended outside of the $M_{200}$ ranges probed
by our simulations, to show the asymptotic behavior at small and large halo
masses. Depending on the redshift and the model, departures in MG $b(M)$ from
$\mathrm{\Lambda\text{CDM}}$ can reach up to $\sim 10\%$. Contrary to the HMF,
MG-induced increase in the strength of gravity lowers the bias, as a result of
enhanced DM clustering. Similar trends have also been reported in (Schmidt
_et al._ , 2009, 2010; Arnold  _et al._ , 2019). The ratios predicted
analytically from the S01 framework do not match the simulation amplitudes
exactly, but they still qualitatively capture the trends.

Further, in Fig. 6, we plot the $b(M)$ in these MG models as a function of
$M_{200}$. The top panels present absolute $b(M)$ values, while the bottom
ones include the ratio between S01 predictions and the simulation-based bias.
Here, we can clearly see that the analytical model matches the simulation
results within $5-10\%$. This affirms our approach in extending S01 to beyond
$\mathrm{\Lambda\text{CDM}}$, for the $f(R)$ and nDGP models we study.

As for the HM build-up, we need a bias prescription for a much wider halo mass
range than what our simulations cover. This overall consistency between the
analytical and simulation results is sufficient for us, and thus we can use
the S01 modeling for $b(M)$ also in our MG variants.

##  Appendix B Comparison of HALOFIT and HMCODE results

Figure 7: Comparison of $P(k)_{\text{MG}}$ obtained from the input of
$P(k)_{\Lambda\text{CDM, HMCODE}}$ (solid lines) and $P(k)_{\Lambda\text{CDM,
HALOFIT}}$ (dashed lines), with the elephant simulation results, for a range
of redshifts as indicated in the legends. The error contours correspond to the
uncertainty in the simulation PS results, and the vertical grey dotted line is
the 5% accuracy regime.

As discussed in Sec. IV.2, for our baseline PS modeling in MG we multiply the
HM-based ratio $\Upsilon(k)$ with halofit PS derived for
$\mathrm{\Lambda\text{CDM}}$. Here we test our approach for the case where the
$\mathrm{\Lambda\text{CDM}}$ PS is obtained from the hmcode Mead  _et al._
(2015) instead. Similarly to halofit, hmcode is also built on the principles
of HM, incorporating however additional corrections in the standard HM build-
up, owing to physical constraints. The parameters of the corrections are based
on high-resolution simulated $\mathrm{\Lambda\text{CDM}}$ power spectra from
the emulator introduced in Ref. Heitmann  _et al._ (2014). Here we use the
latest ‘hmcode-2020’ version333https://github.com/alexander-mead/HMcode Mead
_et al._ (2021).

We compute $P(k)_{\text{MG}}$ by multiplying $\Upsilon(k)_{\text{HM}}$ with
both halofit and hmcode inputs for $P(k)_{\Lambda\text{CDM}}$. Then, in Fig.
7, we compare both predictions with the elephant simulation results. Here we
see a similar performance of both the methods, with some exceptions at small
scales and high-$z$, where hmcode occasionally performs better. Interestingly,
at $z=0$, the halofit framework seems to lead to better results for a range of
$k$-scales. We emphasize however that as our simulations were done for one
particular set of cosmological parameters, these trends between halofit and
hmcode-based predictions could change for other background cosmologies. In any
case, as what we provide is the ratio $\Upsilon(k)$ to be multiplied by the
$\mathrm{\Lambda\text{CDM}}$ PS prediction, one can employ any best-fit
$P(k)_{\Lambda\text{CDM}}$ for the latter to possibly improve the final
accuracy of $P(k)_{\text{MG}}$.

## References

  * Eisenstein  _et al._ (2005) D. J. Eisenstein, I. Zehavi, D. W. Hogg, R. Scoccimarro, M. R. Blanton, and others., ApJ 633, 560 (2005), arXiv:astro-ph/0501171 [astro-ph] .
  * Collaboration (2020) P. Collaboration, A&A 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Cyburt  _et al._ (2016) R. H. Cyburt, B. D. Fields, K. A. Olive, and T.-H. Yeh, Reviews of Modern Physics 88, 015004 (2016), arXiv:1505.01076 [astro-ph.CO] .
  * Alam  _et al._ (2021) S. Alam, M. Aubert, S. Avila, C. Balland, J. E. Bautista, _et al._ , Phys. Rev. D 103, 083533 (2021), arXiv:2007.08991 [astro-ph.CO] .
  * Betoule  _et al._ (2014) M. Betoule, R. Kessler, J. Guy, J. Mosher, D. Hardin, _et al._ , A&A 568, A22 (2014), arXiv:1401.4064 [astro-ph.CO] .
  * Hinshaw  _et al._ (2013) G. Hinshaw, D. Larson, E. Komatsu, D. N. Spergel, C. L. Bennett, _et al._ , ApJS 208, 19 (2013), arXiv:1212.5226 [astro-ph.CO] .
  * Cole  _et al._ (2005) S. Cole, W. J. Percival, J. A. Peacock, P. Norberg, C. M. Baugh, _et al._ , MNRAS 362, 505 (2005), arXiv:astro-ph/0501174 [astro-ph] .
  * Blake  _et al._ (2011) C. Blake, E. A. Kazin, F. Beutler, T. M. Davis, D. Parkinson, S. Brough, M. Colless, C. Contreras, W. Couch, S. Croom, _et al._ , MNRAS 418, 1707 (2011), arXiv:1108.2635 [astro-ph.CO] .
  * Pezzotta  _et al._ (2017) A. Pezzotta, S. de la Torre, J. Bel, B. R. Granett, L. Guzzo, _et al._ , A&A 604, A33 (2017), arXiv:1612.05645 [astro-ph.CO] .
  * Planck Collaboration (2014) Planck Collaboration, A&A 571, A20 (2014), arXiv:1303.5080 [astro-ph.CO] .
  * Bean and Tangmatitham (2010) R. Bean and M. Tangmatitham, Phys. Rev. D 81, 083534 (2010), arXiv:1002.4197 [astro-ph.CO] .
  * Hellwing  _et al._ (2017) W. A. Hellwing, K. Koyama, B. Bose, and G.-B. Zhao, Phys. Rev. D 96, 023515 (2017), arXiv:1703.03395 [astro-ph.CO] .
  * Winther  _et al._ (2017) H. A. Winther, K. Koyama, M. Manera, B. S. Wright, and G.-B. Zhao, J. Cosmology Astropart. Phys 2017, 006 (2017), arXiv:1703.00879 [astro-ph.CO] .
  * Bose  _et al._ (2017) S. Bose, B. Li, A. Barreira, J.-h. He, W. A. Hellwing, _et al._ , J. Cosmology Astropart. Phys 2017, 050 (2017), arXiv:1611.09375 [astro-ph.CO] .
  * Barreira  _et al._ (2015) A. Barreira, S. Bose, and B. Li, J. Cosmology Astropart. Phys 2015, 059 (2015), arXiv:1511.08200 [astro-ph.CO] .
  * Koyama  _et al._ (2009) K. Koyama, A. Taruya, and T. Hiramatsu, Phys. Rev. D 79, 123512 (2009), arXiv:0902.0618 [astro-ph.CO] .
  * DESI Collaboration (2016) DESI Collaboration, arXiv e-prints , arXiv:1611.00036 (2016), arXiv:1611.00036 [astro-ph.IM] .
  * LSST Science Collaboration (2009) LSST Science Collaboration, arXiv e-prints , arXiv:0912.0201 (2009), arXiv:0912.0201 [astro-ph.IM] .
  * Laureijs  _et al._ (2011) R. Laureijs, J. Amiaux, S. Arduini, J. L. Auguères, J. Brinchmann, _et al._ , arXiv e-prints , arXiv:1110.3193 (2011), arXiv:1110.3193 [astro-ph.CO] .
  * Kilbinger  _et al._ (2013) M. Kilbinger, L. Fu, C. Heymans, F. Simpson, J. Benjamin, _et al._ , MNRAS 430, 2200 (2013), arXiv:1212.3338 [astro-ph.CO] .
  * Hikage  _et al._ (2019) C. Hikage, M. Oguri, T. Hamana, S. More, R. Mandelbaum, _et al._ , PASJ 71, 43 (2019), arXiv:1809.09148 [astro-ph.CO] .
  * To  _et al._ (2021) C. To, E. Krause, E. Rozo, H. Wu, D. Gruen, _et al._ , Phys. Rev. Lett. 126, 141301 (2021), arXiv:2010.01138 [astro-ph.CO] .
  * Turner  _et al._ (2023) R. J. Turner, C. Blake, and R. Ruggeri, MNRAS 518, 2436 (2023), arXiv:2207.03707 [astro-ph.CO] .
  * Cooray and Sheth (2002) A. Cooray and R. Sheth, Phys. Rep. 372, 1 (2002), arXiv:astro-ph/0206508 [astro-ph] .
  * Smith  _et al._ (2003) R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, _et al._ , MNRAS 341, 1311 (2003), arXiv:astro-ph/0207664 [astro-ph] .
  * Mead  _et al._ (2015) A. J. Mead, J. A. Peacock, C. Heymans, S. Joudaki, and A. F. Heavens, MNRAS 454, 1958 (2015), arXiv:1505.07833 [astro-ph.CO] .
  * Lawrence  _et al._ (2010) E. Lawrence, K. Heitmann, M. White, D. Higdon, C. Wagner, _et al._ , ApJ 713, 1322 (2010), arXiv:0912.4490 [astro-ph.CO] .
  * Takahashi  _et al._ (2012) R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, ApJ 761, 152 (2012), arXiv:1208.2701 [astro-ph.CO] .
  * Smith and Angulo (2019) R. E. Smith and R. E. Angulo, MNRAS 486, 1448 (2019), arXiv:1807.00040 [astro-ph.CO] .
  * Angulo  _et al._ (2021) R. E. Angulo, M. Zennaro, S. Contreras, G. Aricò, M. Pellejero-Ibañez, and J. Stücker, MNRAS 507, 5869 (2021), arXiv:2004.06245 [astro-ph.CO] .
  * Oyaizu  _et al._ (2008) H. Oyaizu, M. Lima, and W. Hu, Phys. Rev. D 78, 123524 (2008), arXiv:0807.2462 [astro-ph] .
  * Alam  _et al._ (2021) S. Alam, A. Aviles, R. Bean, Y.-C. Cai, M. Cautun, _et al._ , J. Cosmology Astropart. Phys 2021, 050 (2021), arXiv:2011.05771 [astro-ph.CO] .
  * Hu and Sawicki (2007a) W. Hu and I. Sawicki, Phys. Rev. D 76, 104043 (2007a), arXiv:0708.1190 [astro-ph] .
  * Cataneo  _et al._ (2019) M. Cataneo, L. Lombriser, C. Heymans, A. J. Mead, A. Barreira, S. Bose, and B. Li, MNRAS 488, 2121 (2019), arXiv:1812.05594 [astro-ph.CO] .
  * Ruan  _et al._ (2020) C.-Z. Ruan, T.-J. Zhang, and B. Hu, MNRAS 492, 4235 (2020), arXiv:2001.09229 [astro-ph.CO] .
  * Baldi  _et al._ (2014) M. Baldi, F. Villaescusa-Navarro, M. Viel, E. Puchwein, V. Springel, and L. Moscardini, MNRAS 440, 75 (2014), arXiv:1311.2588 [astro-ph.CO] .
  * Winther  _et al._ (2015) H. A. Winther, F. Schmidt, A. Barreira, C. Arnold, S. Bose, _et al._ , MNRAS 454, 4208 (2015), arXiv:1506.06384 [astro-ph.CO] .
  * Puchwein  _et al._ (2013) E. Puchwein, M. Baldi, and V. Springel, MNRAS 436, 348 (2013), arXiv:1305.2418 [astro-ph.CO] .
  * Clifton  _et al._ (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rep. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Baker  _et al._ (2015) T. Baker, D. Psaltis, and C. Skordis, ApJ 802, 63 (2015), arXiv:1412.3455 [astro-ph.CO] .
  * Koyama (2016) K. Koyama, Reports on Progress in Physics 79, 046902 (2016), arXiv:1504.04623 [astro-ph.CO] .
  * Nojiri  _et al._ (2017) S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Phys. Rept. 692, 1 (2017), arXiv:1705.11098 [gr-qc] .
  * Ferreira (2019) P. G. Ferreira, ARA&A 57, 335 (2019), arXiv:1902.10503 [astro-ph.CO] .
  * Sotiriou and Faraoni (2010) T. P. Sotiriou and V. Faraoni, Reviews of Modern Physics 82, 451 (2010), arXiv:0805.1726 [gr-qc] .
  * Baker  _et al._ (2021) T. Baker, A. Barreira, H. Desmond, P. Ferreira, B. Jain, _et al._ , Reviews of Modern Physics 93, 015003 (2021).
  * Koyama (2018) K. Koyama, International Journal of Modern Physics D 27, 1848001 (2018).
  * Joyce  _et al._ (2016) A. Joyce, L. Lombriser, and F. Schmidt, Annual Review of Nuclear and Particle Science 66, 95 (2016), arXiv:1601.06133 [astro-ph.CO] .
  * Bull  _et al._ (2016) P. Bull, Y. Akrami, J. Adamek, T. Baker, E. Bellini, _et al._ , Physics of the Dark Universe 12, 56 (2016), arXiv:1512.05356 [astro-ph.CO] .
  * Seljak (2000) U. Seljak, MNRAS 318, 203 (2000), arXiv:astro-ph/0001493 [astro-ph] .
  * Will (2014) C. M. Will, Living Reviews in Relativity 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Taylor and Weisberg (1982) J. H. Taylor and J. M. Weisberg, ApJ 253, 908 (1982).
  * Abbott  _et al._ (2017) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, _et al._ , Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2016) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, _et al._ , Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Gravity Collaboration  _et al._ (2020) Gravity Collaboration, R. Abuter, A. Amorim, M. Bauböck, J. P. Berger, H. Bonnet, _et al._ , A&A 636, L5 (2020), arXiv:2004.07187 [astro-ph.GA] .
  * Shapiro (1964) I. I. Shapiro, Phys. Rev. Lett. 13, 789 (1964).
  * Joyce  _et al._ (2015) A. Joyce, B. Jain, J. Khoury, and M. Trodden, Phys. Rep. 568, 1 (2015), arXiv:1407.0059 [astro-ph.CO] .
  * Brax (2013) P. Brax, Class. Quant. Grav. 30, 214005 (2013).
  * Sakstein (2018) J. Sakstein, International Journal of Modern Physics D 27, 1848008-228 (2018).
  * Dvali  _et al._ (2000) G. Dvali, G. Gabadadze, and M. Porrati, Physics Letters B 485, 208 (2000), arXiv:hep-th/0005016 [hep-th] .
  * Mead  _et al._ (2016) A. J. Mead, C. Heymans, L. Lombriser, J. A. Peacock, O. I. Steele, and H. A. Winther, MNRAS 459, 1468 (2016), arXiv:1602.02154 [astro-ph.CO] .
  * Zhao (2014) G.-B. Zhao, ApJS 211, 23 (2014), arXiv:1312.1291 [astro-ph.CO] .
  * Tessore  _et al._ (2015) N. Tessore, H. A. Winther, R. B. Metcalf, P. G. Ferreira, and C. Giocoli, J. Cosmology Astropart. Phys 2015, 036 (2015), arXiv:1508.04011 [astro-ph.CO] .
  * Achitouv  _et al._ (2016) I. Achitouv, M. Baldi, E. Puchwein, and J. Weller, Phys. Rev. D 93, 103522 (2016), arXiv:1511.01494 [astro-ph.CO] .
  * Lombriser  _et al._ (2014) L. Lombriser, K. Koyama, and B. Li, J. Cosmology Astropart. Phys 2014, 021 (2014), arXiv:1312.1292 [astro-ph.CO] .
  * Barreira  _et al._ (2014) A. Barreira, B. Li, W. A. Hellwing, L. Lombriser, C. M. Baugh, and S. Pascoli, JCAP 04, 029 (2014), arXiv:1401.1497 [astro-ph.CO] .
  * Barreira  _et al._ (2014) A. Barreira, B. Li, W. A. Hellwing, C. M. Baugh, and S. Pascoli, J. Cosmology Astropart. Phys 2014, 031 (2014), arXiv:1408.1084 [astro-ph.CO] .
  * Schmidt  _et al._ (2010) F. Schmidt, W. Hu, and M. Lima, Phys. Rev. D 81, 063005 (2010), arXiv:0911.5178 [astro-ph.CO] .
  * Hu and Sawicki (2007b) W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007b), arXiv:0705.1158 [astro-ph] .
  * Khoury and Weltman (2004) J. Khoury and A. Weltman, Phys. Rev. D 69, 044026 (2004), arXiv:astro-ph/0309411 [astro-ph] .
  * Deffayet  _et al._ (2002) C. Deffayet, G. Dvali, G. Gabadadze, and A. Vainshtein, Phys. Rev. D 65, 044026 (2002).
  * Vainshtein (1972) A. I. Vainshtein, Physics Letters B 39, 393 (1972).
  * Li  _et al._ (2012) B. Li, G.-B. Zhao, R. Teyssier, and K. Koyama, J. Cosmology Astropart. Phys 2012, 051 (2012), arXiv:1110.1379 [astro-ph.CO] .
  * Li  _et al._ (2013a) B. Li, G.-B. Zhao, and K. Koyama, J. Cosmology Astropart. Phys 2013, 023 (2013a), arXiv:1303.0008 [astro-ph.CO] .
  * Naidoo  _et al._ (2022) K. Naidoo, W. Hellwing, M. Bilicki, N. Libeskind, S. Pfeifer, and Y. Hoffman, arXiv e-prints , arXiv:2209.14386 (2022), arXiv:2209.14386 [astro-ph.CO] .
  * Lewis and Challinor (2011) A. Lewis and A. Challinor, “CAMB: Code for Anisotropies in the Microwave Background,” (2011), ascl:1102.026 .
  * Colombi and Novikov (2011) S. Colombi and D. Novikov, “POWMES: Measuring the Power Spectrum in an N-body Simulation,” Astrophysics Source Code Library, record ascl:1110.017 (2011), ascl:1110.017 .
  * Peacock and Smith (2000) J. A. Peacock and R. E. Smith, MNRAS 318, 1144 (2000), arXiv:astro-ph/0005010 [astro-ph] .
  * Cooray and Hu (2001) A. Cooray and W. Hu, ApJ 554, 56 (2001), arXiv:astro-ph/0012087 [astro-ph] .
  * Watson  _et al._ (2013) W. A. Watson, I. T. Iliev, A. D’Aloisio, A. Knebe, P. R. Shapiro, and G. Yepes, MNRAS 433, 1230 (2013), arXiv:1212.0095 [astro-ph.CO] .
  * Gupta  _et al._ (2022) S. Gupta, W. A. Hellwing, M. Bilicki, and J. E. García-Farieta, Phys. Rev. D 105, 043538 (2022), arXiv:2112.03699 [astro-ph.CO] .
  * Koyama and Silva (2007) K. Koyama and F. P. Silva, Phys. Rev. D 75, 084040 (2007), arXiv:hep-th/0702169 [hep-th] .
  * Sheth  _et al._ (2001) R. K. Sheth, H. J. Mo, and G. Tormen, MNRAS 323, 1 (2001), arXiv:astro-ph/9907024 [astro-ph] .
  * Ludlow  _et al._ (2016) A. D. Ludlow, S. Bose, R. E. Angulo, L. Wang, W. A. Hellwing, _et al._ , MNRAS 460, 1214 (2016), arXiv:1601.02624 [astro-ph.CO] .
  * Mitchell  _et al._ (2019) M. A. Mitchell, C. Arnold, J.-h. He, and B. Li, MNRAS 487, 1410 (2019), arXiv:1901.06392 [astro-ph.CO] .
  * Mitchell  _et al._ (2018) M. A. Mitchell, J.-h. He, C. Arnold, and B. Li, MNRAS 477, 1133 (2018), arXiv:1802.02165 [astro-ph.CO] .
  * Mitchell  _et al._ (2021) M. A. Mitchell, C. Hernández-Aguayo, C. Arnold, and B. Li, MNRAS 508, 4140 (2021), arXiv:2106.13815 [astro-ph.CO] .
  * Press and Schechter (1974) W. H. Press and P. Schechter, ApJ 187, 425 (1974).
  * Bond  _et al._ (1991) J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser, ApJ 379, 440 (1991).
  * Jenkins  _et al._ (2001) A. Jenkins, C. S. Frenk, S. D. M. White, J. M. Colberg, S. Cole, _et al._ , MNRAS 321, 372 (2001), arXiv:astro-ph/0005260 [astro-ph] .
  * Warren  _et al._ (2006) M. S. Warren, K. Abazajian, D. E. Holz, and L. Teodoro, ApJ 646, 881 (2006), arXiv:astro-ph/0506395 [astro-ph] .
  * Despali  _et al._ (2016) G. Despali, C. Giocoli, R. E. Angulo, G. Tormen, R. K. Sheth, G. Baso, and L. Moscardini, MNRAS 456, 2486 (2016), arXiv:1507.05627 [astro-ph.CO] .
  * Tinker  _et al._ (2008) J. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. Warren, _et al._ , ApJ 688, 709 (2008), arXiv:0803.2706 [astro-ph] .
  * Angulo  _et al._ (2012) R. E. Angulo, V. Springel, S. D. M. White, A. Jenkins, C. M. Baugh, and C. S. Frenk, MNRAS 426, 2046 (2012), arXiv:1203.3216 [astro-ph.CO] .
  * Hellwing and Juszkiewicz (2009) W. A. Hellwing and R. Juszkiewicz, Phys. Rev. D 80, 083522 (2009), arXiv:0809.1976 [astro-ph] .
  * Lam and Li (2012) T. Y. Lam and B. Li, MNRAS 426, 3260 (2012), arXiv:1205.0059 [astro-ph.CO] .
  * Falck  _et al._ (2015) B. Falck, K. Koyama, and G.-B. Zhao, J. Cosmology Astropart. Phys 2015, 049 (2015), arXiv:1503.06673 [astro-ph.CO] .
  * Seljak and Warren (2004) U. Seljak and M. S. Warren, MNRAS 355, 129 (2004), arXiv:astro-ph/0403698 [astro-ph] .
  * Tinker  _et al._ (2010) J. L. Tinker, B. E. Robertson, A. V. Kravtsov, A. Klypin, M. S. Warren, _et al._ , ApJ 724, 878 (2010), arXiv:1001.3162 [astro-ph.CO] .
  * Comparat  _et al._ (2017) J. Comparat, F. Prada, G. Yepes, and A. Klypin, MNRAS 469, 4157 (2017), arXiv:1702.01628 [astro-ph.CO] .
  * Navarro  _et al._ (1996) J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 462, 563 (1996), arXiv:astro-ph/9508025 [astro-ph] .
  * Power  _et al._ (2003) C. Power, J. F. Navarro, A. Jenkins, C. S. Frenk, S. D. M. White, _et al._ , MNRAS 338, 14 (2003), arXiv:astro-ph/0201544 [astro-ph] .
  * Tinker  _et al._ (2005) J. L. Tinker, D. H. Weinberg, Z. Zheng, and I. Zehavi, ApJ 631, 41 (2005), arXiv:astro-ph/0411777 [astro-ph] .
  * Mead and Verde (2021) A. J. Mead and L. Verde, MNRAS 503, 3095 (2021), arXiv:2011.08858 [astro-ph.CO] .
  * Schmidt (2016) F. Schmidt, Phys. Rev. D 93, 063512 (2016), arXiv:1511.02231 [astro-ph.CO] .
  * Mead  _et al._ (2020) A. J. Mead, T. Tröster, C. Heymans, L. Van Waerbeke, and I. G. McCarthy, A&A 641, A130 (2020), arXiv:2005.00009 [astro-ph.CO] .
  * Schmidt (2009) F. Schmidt, Phys. Rev. D 80, 123003 (2009), arXiv:0910.0235 [astro-ph.CO] .
  * Li  _et al._ (2013b) B. Li, W. A. Hellwing, K. Koyama, G.-B. Zhao, E. Jennings, and C. M. Baugh, MNRAS 428, 743 (2013b), arXiv:1206.4317 [astro-ph.CO] .
  * Mead  _et al._ (2021) A. J. Mead, S. Brieden, T. Tröster, and C. Heymans, MNRAS 502, 1401 (2021), arXiv:2009.01858 [astro-ph.CO] .
  * Tassev  _et al._ (2013) S. Tassev, M. Zaldarriaga, and D. J. Eisenstein, J. Cosmology Astropart. Phys 2013, 036 (2013), arXiv:1301.0322 [astro-ph.CO] .
  * Heitmann  _et al._ (2006) K. Heitmann, D. Higdon, C. Nakhleh, and S. Habib, ApJ 646, L1 (2006), arXiv:astro-ph/0606154 [astro-ph] .
  * Habib  _et al._ (2007) S. Habib, K. Heitmann, D. Higdon, C. Nakhleh, and B. Williams, Phys. Rev. D 76, 083503 (2007), arXiv:astro-ph/0702348 [astro-ph] .
  * Ramachandra  _et al._ (2021) N. Ramachandra, G. Valogiannis, M. Ishak, K. Heitmann, and LSST Dark Energy Science Collaboration, Phys. Rev. D 103, 123525 (2021), arXiv:2010.00596 [astro-ph.CO] .
  * Winther  _et al._ (2019) H. A. Winther, S. Casas, M. Baldi, K. Koyama, B. Li, _et al._ , Phys. Rev. D 100, 123540 (2019), arXiv:1903.08798 [astro-ph.CO] .
  * Arnold  _et al._ (2022) C. Arnold, B. Li, B. Giblin, J. Harnois-Déraps, and Y.-C. Cai, MNRAS 515, 4161 (2022), arXiv:2109.04984 [astro-ph.CO] .
  * Harnois-Déraps  _et al._ (2022) J. Harnois-Déraps, C. Hernandez-Aguayo, C. Cuesta-Lazaro, C. Arnold, B. Li, _et al._ , arXiv e-prints , arXiv:2211.05779 (2022), arXiv:2211.05779 [astro-ph.CO] .
  * Ruan  _et al._ (2023) C.-Z. Ruan, C. Cuesta-Lazaro, A. Eggemeier, B. Li, C. M. Baugh, _et al._ , arXiv e-prints , arXiv:2301.02970 (2023), arXiv:2301.02970 [astro-ph.CO] .
  * Acuto  _et al._ (2021) A. Acuto, I. G. McCarthy, J. Kwan, J. Salcido, S. G. Stafford, and A. S. Font, MNRAS 508, 3519 (2021), arXiv:2109.11855 [astro-ph.CO] .
  * Arnold  _et al._ (2019) C. Arnold, P. Fosalba, V. Springel, E. Puchwein, and L. Blot, MNRAS 483, 790 (2019), arXiv:1805.09824 [astro-ph.CO] .
  * Chisari  _et al._ (2019) N. E. Chisari, A. J. Mead, S. Joudaki, P. G. Ferreira, A. Schneider, _et al._ , The Open Journal of Astrophysics 2, 4 (2019), arXiv:1905.06082 [astro-ph.CO] .
  * Kopp  _et al._ (2013) M. Kopp, S. A. Appleby, I. Achitouv, and J. Weller, Phys. Rev. D 88, 084015 (2013), arXiv:1306.3233 [astro-ph.CO] .
  * Schmidt  _et al._ (2009) F. Schmidt, M. Lima, H. Oyaizu, and W. Hu, Phys. Rev. D 79, 083518 (2009), arXiv:0812.0545 [astro-ph] .
  * Heitmann  _et al._ (2014) K. Heitmann, E. Lawrence, J. Kwan, S. Habib, and D. Higdon, ApJ 780, 111 (2014), arXiv:1304.7849 [astro-ph.CO] .


# Using the SourceXtractor++ package for data reduction

M. Kümmel,1 A. Álvarez-Ayllón,2 E. Bertin,3,4 P. Dubath,2 R. Gavazzi,4,5 W.
Hartley,2 and M. Schefer,2 on behalf of the Euclid Consortium

###### Abstract

The Euclid satellite is an ESA mission scheduled for launch in September 2023.
To optimally perform critical stages of the data reduction, such as object
detection and morphology determination, a new and modern software package was
required. We have developed SourceXtractor++ as open source software for
detecting and measuring sources in astronomical images. It is a complete
redesign of the original SExtractor, written mainly in C++. The package
follows a modular approach and facilitates the analysis of multiple
overlapping sources over many images with different pixel grids.
SourceXtractor++ is already operational in many areas of the Euclid
processing, and we demonstrate here the capabilities of the current version
v0.19 on the basis of a set of typical use cases, which are available for
download.

1LMU Faculty of Physics, Scheinerstr. 1, 81679 München, Germany;
<EMAIL_ADDRESS>
2Department of Astronomy, University of Geneva, Chemin d’Écogia 16, CH-1290,
Versoix, Switzerland

3Canada-France-Hawaii Telescope, 65-1238 Mamalahoa Hwy, Kamuela, HI 96743, USA

4Sorbonne Université, CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98
bis bd Arago, F-75014 Paris, France

5LAM, Aix-Marseille Université, CNRS, CNES, Marseille, France

##  1 Introduction

Euclid (Laureijs et al. 2011) will cover $15\,000\,{\rm deg^{2}}$ and
$40\,{\rm deg^{2}}$ with its Wide Survey and Deep Survey, respectively and is
expected to detect and analyse $\approx 10^{9}$ individual objects. To
contribute to this task we have developed SourceXtractor++ (Kümmel et al.
2020; Bertin et al. 2020) as a replacement for and an extension of the
original SExtractor (Bertin & Arnouts 1996) software. In the Euclid processing
SourceXtractor++ is being used in the object detection (Kümmel et al. 2015)
and the photometric validation of ground-based data. The software has excelled
in the Euclid Morphology Challenge (Bretonnière & et al. 2022; Merlin & et al.
2022), a recent comparison of several up-to-date model fitting solutions, and
SourceXtractor++ was selected to perform flexible model fitting with Sérsic
and Disk+Sérsic models to all Euclid objects. SourceXtractor++ is
distributed111https://github.com/astrorama/SourceXtractorPlusPlus as an
independent package via source code or various binary distributions (conda,
Fedora/Centos RPMs) to the astronomy community.

##  2 Test Data Set

SourceXtractor++ can be used in a variety of ways to work on different kinds
of input data (single band, multi-band, data cube). To demonstrate this
flexibility we have assembled a series of use cases
222https://cloud.physik.lmu.de/index.php/s/3K4KemBsw5y9yqd. The use cases
contain all necessary input data, the configuration files and, as a reference,
the expected results. The exercises are organized in several sub-directories:

  * •

01_basics: Basic operations such as generating a blank ASCII configuration
file, listing all available object properties or performing a simple object
detection as in SExtractor2 are shown. Different image extensions for
detection and pixel weights are accessed. Sérsic models are fit to all
detected objects. All settings for the model fitting are done in a python
configuration file which is executed at the start of the fitting process.

  * •

02_mband: Disk+Bulge models are simultaneously fit to three measurement images
while using a fourth image for the object detection. The measurement images
are grouped according to the filter information provided in the header, which
results in one independent flux measurement per group. The images coming from
JWST cover different areas and have diverse pixel scales. All relevant
positions and parameters are negotiated with the image WCS. Arbitrarily
defined dependent parameters can be pushed to the output. The usage of python
for the configuration file facilitates pre-processing the data, e.g. band
selection or keyword adjustments;

  * •

03_priors: This example was taken from Quilley & de Lapparent (2022) and
illustrates the refined multi-band fitting of a bright galaxy. With priors the
fit results can be constrained to avoid the influence of neighbour objects
such as bright stars. Priors in connection with dependent parameters allow
large flexibility such as fitting band specific positions while maintaining
meaningful results by keeping the position differences small.

Figure 1.: Fitting a Sérsic+Bar+Ring model to M95 with user defined functions.
The SDSS i-band image is shown in the left, the fitted model in the middle and
the residual image in the right sub-panel.

  * •

04_user_defined: The functional forms, Sérsic, exponential disk, deVaucouleur
Bulge and point-like are provided in SourceXtractor++; the exercise shows how
users can define an ONNX model to define customized functional forms.
SourceXtractor++ uses the ONNX Run time333https://onnxruntime.ai to evaluate
the ONNX models and fit higher level features such as rings or bars. Figure 1
shows the fit of a Sérsic+Bar+Ring to the SDSS i-band data of M95.

  * •

05_cubes: SourceXtractor++ can run on multi-extension FITS images and even on
data cubes. The example uses a broad-band image for detection and then
performs model fitting on several image groups selected from layers in a data
cube from the MUSE/ESO instrument observing the Abell 2744 galaxy
cluster444http://muse-vlt.eu/science/a2744. The layers are selected by
wavelength range in the python configuration file using the WCS of the cube.

  * •

06_extended_assocs: SourceXtractor++ has implemented an extended association
mode which allows the user to funnel external object properties such as
starting positions or shape parameters into the model fitting process. This
mode is a first step towards disentangling the detection and the measurement
stages. Once fully developed it will be possible to run the measurements
solely based on objects properties (positions, segmentation) provided to
SourceXtractor++.

  * •

07_python_control This exercise runs currently only in a special branch of the
v0.20 development version. The sequence of processing steps (detection,
grouping, measurement) is executed directly from Python. Throughout the
processing the user has access to the objects and other processing items and
can apply selections or run their own analysis code on the objects.
SourceXtractor++ can then integrate completely into the more complex
processing pipeline of the user.

Figure 2.: The efficiency of the multi-core processing in SourceXtractor++.

##  3 Status and Outlook

SourceXtractor++ has been available for a few years and has also picked up
users outside of the Euclid community. In addition to the Euclid requirements
we designed SourceXtractor++ as a general processing tool to offer novel
possibilities for reducing data. The efficient parallel processing in
SourceXtractor++ is evidenced by the rather flat curve in Figure 2, which
shows the number of objects processed per hour and core as a function of the
number of processors assigned to the task. SourceXtractor++ contains many of
the features we had in mind when starting the project. Future developments
will concentrate on becoming more Pythonic and separating detection and
measuring, as indicated in the last two exercises. We are certain that
SourceXtractor++ offers new possibilities for reducing data both within and
outside of the Euclid project.

### Acknowledgments

The Euclid Consortium acknowledges the European Space Agency and a number of
agencies and institutes that have supported the development of Euclid, in
particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian
Science Policy, the Canadian Euclid Consortium, the French Centre National
d’Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish
Space Research Institute, the Fundação para a Ciência e a Tecnologia, the
Ministerio de Ciencia e Innovación, the National Aeronautics and Space
Administration, the National Astronomical Observatory of Japan, the
Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the
Romanian Space Agency, the State Secretariat for Education, Research and
Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom
Space Agency. A complete and detailed list is available on the Euclid web site
(http://www.euclid-ec.org).

## References

  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
  * Bertin et al. (2020) Bertin, E., Schefer, M., Apostolakos, N., Álvarez-Ayllón, A., Dubath, P., & Kümmel, M. 2020, in Astronomical Data Analysis Software and Systems XXIX, edited by R. Pizzo, E. R. Deul, J. D. Mol, J. de Plaa, & H. Verkouter, vol. 527 of Astronomical Society of the Pacific Conference Series, 461
  * Bretonnière & et al. (2022) Bretonnière, H., & et al. 2022, arXiv e-prints. 2209.12907
  * Kümmel et al. (2020) Kümmel, M., Bertin, E., Schefer, M., Apostolakos, N., Álvarez-Ayllón, A., & Dubath, P. 2020, in Astronomical Data Analysis Software and Systems XXIX, edited by R. Pizzo, E. R. Deul, J. D. Mol, J. de Plaa, & H. Verkouter, vol. 527 of Astronomical Society of the Pacific Conference Series, 29
  * Kümmel et al. (2015) Kümmel, M., Mohr, J., Fontana, A., Dole, H., Boucaud, A., Cabanac, R., Castellano, M., Petrinca, P., & Pilo, S. 2015, in Astronomical Data Analysis Software an Systems XXIV (ADASS XXIV), edited by A. R. Taylor, & E. Rosolowsky, vol. 495 of Astronomical Society of the Pacific Conference Series, 249
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, et al. 2011, arXiv e-prints. 1110.3193
  * Merlin & et al. (2022) Merlin, E., & et al. 2022, arXiv e-prints. 2209.12906
  * Quilley & de Lapparent (2022) Quilley, L., & de Lapparent, V. 2022, A&A, 666, A170