# On the impact of $f(Q)$ gravity on the Large Scale Structure

Oleksii Sokoliuk1, Simran Arora2, Subhrat Praharaj2, Alexander Baransky3, and
P.K. Sahoo2
1Astronomical Observatory of the National Academy of Sciences of Ukraine (MAO
NASU), Kyiv, 03143, Ukraine
2Department of Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India
3Astronomical Observatory, Taras Shevchenko National University of Kyiv, 3
Observatorna St., 04053 Kyiv, Ukraine
oleksii.sokoliuk<EMAIL_ADDRESS>sp7437473@gmail.comabaransky<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We investigate the exponential $f(Q)$ symmetric teleparallel gravitation,
namely $f(Q)=Q+\alpha Q_{0}(1-e^{-\beta\sqrt{Q/Q_{0}}})$ using ME-GADGET code
to probe the structure formation with box sizes $L_{\mathrm{box}}=10/100$
Mpc$/h$ and middle resolution $N_{p}^{1/3}=512$. To reproduce viable cosmology
within the aforementioned modified gravity theory, we first perform Markov
Chain Monte Carlo (MCMC) sampling on OHD/BAO/Pantheon datasets and constrain a
parameter space. Furthermore, we also derive theoretical values for
deceleration parameter $q(z)$, statefinder pair $\\{r,s\\}$ and effective
gravitational constant $G_{\mathrm{eff}}$, perform $Om(z)$ diagnostics. While
carrying out N-body+SPH simulations, we derive CDM+baryons over
density/temperature/mean molecular weight fields, matter power spectrum (both
2/3D, with/without redshift space distortions), bispectrum, two-point
correlation function and halo mass function. Results for small and big
simulation box sizes are therefore properly compared, halo mass function is
related to the Seth-Tormen theoretical prediction and matter power spectrum to
the standard CAMB output.

###### keywords:

Dark energy – Observations – Large-scale structure of Universe

††pubyear: 2015††pagerange: On the impact of $f(Q)$ gravity on the Large Scale
Structure–On the impact of $f(Q)$ gravity on the Large Scale Structure

##  1 Introduction

Numerous independent cosmological observable predicts that universe undergoes
the accelerated expansion phase at the present time (Riess et al., 1998;
Perlmutter et al., 1999; Suzuki et al., 2012; Hinshaw et al., 2013). It is
well-known that General Theory of Relativity (GR) is a quite successful theory
on various cosmological scales, and it is able to describe the recent
accelerated expansion of the universe by introducing the so-called
cosmological constant (or $\Lambda$ term) in the Einstein-Hilbert action
integral. But, such term gives a rise to various issues, that were found in
the papers of (SAHNI & STAROBINSKY, 2000; Padmanabhan, 2003). Since the
gravitational Lagrangian is not practically restricted by only the linear
Ricci scalar term, one can introduce additional terms to emulate effective
dark energy and reproduce different universe evolutionary phases, such as
cosmological inflation or late time accelerated expansion in order to overcome
the aforementioned problems.

There are a lot of different ways to modify general relativity - for example,
by introducing some matter fields (canonical scalar field, vector and gauge
boson fields, Dirac spinors, etc.). Another way is to present an entirely
different notion of Lorentzian 4-manifold curvature by adjusting the metric-
affine connection (De Felice & Tsujikawa, 2010; Capozziello & De Laurentis,
2011). For example, one could use the so-called torsion or non-metricity,
which are constructed based on Weitzenböck and metric incompatible affine
connections, respectively. Consequently, there are two analoges to GR, namely
the Teleparallel Equivalent of GR (TEGR), introduced in (Hayashi & Shirafuji,
1979; Abedi et al., 2017) and Symmetric Teleparallel Equivalent of GR (STEGR)
(Nester & Yo, 1999; Hohmann, 2021). In the current work, we will focus on the
arbitrary parameterization of STEGR ($f(Q)$ gravitation) in particular.

A key aspect of $f(Q)$ theory is the usage of a flat connection pertaining to
the existence of affine coordinates in which all of its components vanish,
converting covariant derivatives into partial derivatives (Hohmann, 2021;
Dimakis et al., 2022; Zhao, 2022). So, it is possible to distinguish gravity
from inertial effects in $f(Q)$ theory. For many modified gravity theories,
the development of the $f(Q)$ theory provides a fresh starting point.
Additionally, it offers a straightforward formulation in which self-
accelerating solutions spontaneously appear in both the early and late
universe. When compared to other geometric extensions of GR, both $f(T)$ and
$f(Q)$ theories have a substantial benefit in that the background field
equations are always of second order, which means that Ostrogradsky’s theorem
(Motohashi & Suyama, 2015) related instability issues are avoided.

Up to this moment, $f(Q)$ gravity has been incorporated in dozens of studies
and is a very promising theory, that can reproduce the behaviour of both early
(De et al., 2022) and late universe, satisfy constraints from Cosmic Microwave
Background (CMB), SuperNovae (SN) distance modulus, Baryon Acoustic
Oscillations (BAO), Observational Hubble Dataset (OHD) and primordial scalar
index $n_{s}$, standard sirens from LIGO/VIRGO/ET (D’Agostino & Nunes, 2022;
Ferreira et al., 2022). For instance, exponential gravitation were constrained
in the study by (Anagnostopoulos et al., 2021; Atayde & Frusciante, 2021), and
the authors found out that such a theory can challenge concordance
$\Lambda$CDM theory. There was also some studies carried out on the matter of
$f(Q)$ cosmography (Mandal et al., 2020b) and energy conditions (Mandal et
al., 2020a; Koussour et al., 2022). Additionally, observational constraints on
the $f(Q)$ gravity have been established for a number of parameterizations of
the $f(Q)$ function using various observational probes (Lazkoz et al., 2019;
Solanki et al., 2023; Jiménez et al., 2020b). In the context of f(Q)
cosmology, a Hamiltonian formulation has been designed to carry out a
canonical quantization procedure (Dimakis et al., 2021). Aside from these
findings, $f(Q)$ gravity has been the focus of several investigations in
varied applications (Hu et al., 2022; Albuquerque & Frusciante, 2022; Wang et
al., 2022; Esposito et al., 2022; Bajardi et al., 2020; Arora & Sahoo, 2022;
Harko et al., 2018). In the current study, we are going to investigate the
logarithmic $f(Q)$ gravity in terms of observational constraints using Markov
Chain Monte Carlo (MCMC) methodologies and high resolution N-body simulations,
which will be discussed in the following subsections.

###  1.1 N-body simulation as a probe of modified gravity

To probe the validity of a particular modified theory of gravitation, one
needs to incorporate various cosmological observables, ranging from cosmic
expansion rate to clustering and structure formation history. The latter could
be most effectively studied with the use of the so-called N-body simulations,
that are well-known to be the best theoretical probe of the large scale
structure of the universe, that provide information on the matter power
spectrum/bispectrum, $N$-point correlation functions and halo mass function,
void size function etc. Over the last few years, such approach has attracted
some interest in the field of modified gravity (see the work of (Hassani &
Lombriser, 2020)). Authors of the paper (Wilson & Bean, 2022) developed a
pipeline to differentiate modified gravity theories from fiducial $\Lambda$CDM
model and constrain those theories properly using voids in the N-body
simulations, that are known to be less affected by the non-linear and baryonic
physics in relation to the dark matter halos. Besides, in addition to voids,
intrinsic shape alignments of massive halos and galaxy/halo angular clustering
could be used to discriminate MOG theories from $\Lambda$CDM in the presence
of massive neutrino (see (Lee et al., 2022) and (Drozda et al., 2022)
respectively). Aforementioned Halo Mass Function (further-HMF) were examined
for $f(Q)$ and Dvali-Gabadadze-Porrati (DGP) gravities in (Gupta et al.,
2022). Widely used code MG-Gadget, introduced and developed in (Puchwein et
al., 2013) were employed to study $f(R)$ Hu-Sawicki theory (Arnold et al.,
2016, 2015; Giocoli et al., 2018), conformally coupled gravity (Ruan et al.,
2022). In turn, we are going to use ME-GADGET code (for documentation, check
(Zhang et al., 2018a)) to study $f(Q)$ gravity behaviour. Such code was
applied to the case of $f(T)$ teleparallel theory (Huang et al., 2022),
interacting dark energy (Zhang et al., 2018a; Zhao et al., 2022; Liu et al.,
2022) and cubic vector gallileon (Su et al., 2022).

Our paper is organised as follows: in the first Section (1) we provide a
little introduction into the topic of modified theories of gravity and N-body
simulations. Consequently, in the Section (2) we present the foundations of
symmetric teleparallel gravity an it’s arbitrary parameterization, in the
third section we adopt FLRW isotropic line element, derive field equations for
our logarithmic choice of $f(Q)$ function. In the Section (4) we therefore
introduce each observational dataset of our consideration and perform MCMC
analysis, in Section (5) we analyze the provided constraints deriving
theoretical predictions for deceleration parameter, statefinder pair and
$Om(z)$. In the following section we set up the ME-GADGET suite and study the
N-body output for small simulation box size, in (7) we therefore compare
aforementioned results with the ones, obtained for large $L_{\mathrm{box}}$.
Finally, in the last section we present the concluding remarks on the key
topics of our study.

##  2 Modified symmetric teleparallel gravitation

Firstly we are going to start by introducing the fundamentals of the symmetric
teleparallel theories of gravitation. In such theories, it is generally
assumed that the scalar curvature of the manifold does vanishes (and
therefore, $R=0$) as well as torsion, however non-metricity is non-zero (and
describes gravitational interactions). Within the symmetric teleparallel and
related theories, affine connection is metric-incompatible such that
$\nabla_{\mu}g_{\alpha\beta}\neq 0$. In order to present the formalism of
symmetric teleparallel theory, one must firstly define the generalized metric
affine connection (Lin & Zhai, 2021):

$\Gamma^{\alpha}_{\,\,\,\mu\nu}=\widetilde{\Gamma}^{\alpha}_{\,\,\,\mu\nu}+K^{\alpha}_{\,\,\,\mu\nu}+L^{\alpha}_{\,\,\,\mu\nu}.$
(1)

In the equation above $\widetilde{\Gamma}^{\alpha}_{\,\,\,\mu\nu}$ is the
usual Levi-Cevita metric-affine connection, that is widely used within the
General Theory of Relativity:

$\widetilde{\Gamma}^{\alpha}_{\,\,\,\mu\nu}=\frac{1}{2}g^{\alpha\beta}(\partial_{\mu}g_{\beta\nu}+\partial_{\nu}g_{\beta\mu}-\partial_{\beta}g_{\mu\nu}).$
(2)

While other two terms in (1) are namely contortion and deformation tensors and
they could be written below as follows:

$\displaystyle K^{\alpha}_{\,\,\,\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}g^{\alpha\beta}\left(T_{\mu\beta\nu}+T_{\nu\beta\mu}+T_{\beta\mu\nu}\right),$
(3) $\displaystyle L^{\alpha}_{\,\,\,\mu\nu}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}g^{\alpha\beta}\left(Q_{\mu\beta\nu}+Q_{\nu\beta\mu}-Q_{\beta\mu\nu}\right).$
(4)

Here
$T^{\alpha}_{\,\,\,\mu\nu}=\Gamma^{\alpha}_{\,\,\,\mu\nu}-\Gamma^{\alpha}_{\,\,\,\nu\mu}$
is the torsion tensor and (Capozziello & D’Agostino, 2022)

$Q_{\alpha\mu\nu}=\nabla_{\alpha}g_{\mu\nu}=\partial_{\alpha}g_{\mu\nu}-\Gamma^{\beta}_{\,\,\,\alpha\mu}g_{\beta\nu}-\Gamma^{\beta}_{\,\,\,\alpha\nu}g_{\mu\beta}.$
(5)

Is obviously the non-metricity tensor. As we already mentioned, here and
further we assume that both Ricci scalar curvature and torsion terms vanish
and therefore we are left with only non-metricity. Therefore, to proceed with
STEGR case one could derive the non-metricity scalar (fundamental quantity)
from non-metricity tensor and its independent traces
$Q_{\alpha}={Q_{\alpha}^{\,\,\,\mu}}_{\mu}$ and
$\tilde{Q}^{\alpha}=Q_{\mu}^{\,\,\,\alpha\mu}$(Lin & Zhai, 2021):

$Q=-g^{\mu\nu}(L^{\alpha}_{\,\,\,\beta\nu}L^{\beta}_{\,\,\,\mu\alpha}-L^{\beta}_{\,\,\,\alpha\beta}L^{\alpha}_{\,\,\,\mu\nu})=-P^{\alpha\beta\gamma}Q_{\alpha\beta\gamma},$
(6)

Where deformation tensor $\mathbf{L}$ was already defined previously and
superpotential could be expressed in the following way:

$P^{\alpha}_{\,\,\,\mu\nu}=\frac{1}{4}\bigg{[}2Q^{\alpha}_{\,\,(\mu\nu)}-Q^{\alpha}_{\,\,\,\mu\nu}+Q^{\alpha}g_{\mu\nu}-\delta^{\alpha}_{(\mu}Q_{\nu)}-\tilde{Q}^{\alpha}g_{\mu\nu}\bigg{]}.$
(7)

Here symmetric and antisymmetric parts of the tensor are:

$\displaystyle F_{(\mu\nu)}=\frac{1}{2}\bigg{(}F_{\mu\nu}+F_{\nu\mu}\bigg{)},$
(8) $\displaystyle
F_{[\mu\nu]}=\frac{1}{2}\bigg{(}F_{\mu\nu}-F_{\nu\mu}\bigg{)}.$ (9)

The condition of symmetric teleparallelism makes the generic affine connection
to be inertial. The most general connection is

$\Gamma^{\alpha}_{\,\,\,\mu\nu}=\frac{\partial
x^{\alpha}}{\partial\xi^{\sigma}}\frac{\partial^{2}\xi^{\sigma}}{\partial
x^{\mu}\partial x^{\nu}},$ (10)

where $\xi^{\sigma}$ is an arbitrary function of spacetime position. We can
always choose a coordinate $x^{\alpha}=\xi^{\sigma}$ bu utilizing a general
coordinate transformation, where the general affine connection
$\Gamma^{\alpha}_{\,\,\,\mu\nu}=0$. We call this coordinate the coincident
gauge (Jiménez et al., 2018a). Thus, in the coincident gauge, we will have
$Q_{\alpha\mu\nu}=\partial_{\alpha}g_{\mu\nu}$, i.e. all the covariant
derivatives are identical to ordinary derivatives.

That was the fundamentals of symmetric teleparallel analogue of the General
Theory of Relativity (GR). Now we are going to present the formalism of
modified symmetric teleparallel cosmology. Einstein-Hilbert action integral of
the aforementioned theory of gravity is therefore could be written down as
follows (Jiménez et al., 2018a):

$\mathcal{S}[g,\Gamma,\Psi_{i}]=\frac{1}{16\pi
G}\int_{\mathcal{M}}-d^{4}xef(Q)+\mathcal{S}_{\mathrm{M}}[g,\Gamma,\Psi_{i}].$
(11)

In the equation above, $\mathcal{M}$ is the four dimensional Lorentzian
manifold that we work on, $g=\det g_{\mu\nu}=\prod_{\mu,\mu}g_{\mu\nu}$ is the
metric tensor determinant, $e=\sqrt{-g}$ and $f(Q)$ is the arbitrary function
of non-metricity scalar, that defines the modified theory of gravitation.
Moreover, $\Gamma$ is the curvature free affine connection and
$\mathcal{S}_{\mathrm{M}}[g,\Gamma,\Psi_{i}]$ defines the contribution of
additional matter fields $\Psi_{i}$ to the total Einstein-Hilbert action
integral. The reason for the above action and specific selection of the non-
metricity scalar is that GR is recreated, up to a Boundary term for the choice
$f=Q$, i.e., for this choice, we recover the allegedly “symmetric teleparallel
equivalent of GR". By varying the action (11) with respect to the metric
tensor inverse $g^{\mu\nu}$ (using least action principle
$\delta\mathcal{S}=0$) we could obtain the corresponding field equations

$\frac{2}{\sqrt{-g}}\nabla_{\alpha}\left(\sqrt{-g}f_{Q}P^{\alpha\mu}_{\,\,\,\,\,\,\,\nu}\right)+\frac{1}{2}\delta^{\mu}_{\,\,\,\nu}f+f_{Q}P^{\mu\nu\beta}Q_{\nu\alpha\beta}=T^{\mu}_{\,\,\,\,\nu},$
(12)

Where $f_{Q}=\frac{\partial{f}}{\partial{Q}}$ and the energy-momentum tensor
could be easily derived from the variation of matter fields Lagrangian
density:

$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}\mathcal{L}_{\mathrm{M}})}{\delta
g^{\mu\nu}}.$ (13)

The connection of equation of motion can be computed by noticing that the
variation of the connection with respect to $\xi^{\alpha}$ is equivalent to
performing a diffeomorphism so that
$\partial_{\xi}\Gamma^{\alpha}_{\,\,\mu\nu}=-\mathcal{L}_{\xi}\Gamma^{\alpha}_{\,\,\mu\nu}=-\nabla_{\mu}\nabla_{\nu}\,\xi^{\alpha}$
(Jiménez et al., 2020b). Besides this, in the absence of hypermomentum, one
can take the variation of equation (11) with respect to connection

$\nabla_{\mu}\nabla_{\nu}\left(\sqrt{-g}f_{Q}P^{\mu\nu}_{\,\,\,\,\,\alpha}\right)=0.$
(14)

For the metric and connection equations, one can notice that
$\mathcal{D}_{\mu}T^{\mu}_{\,\,\,\,\nu}=0$, where $\mathcal{D}_{\mu}$ is the
metric-covariant derivative.

Therefore, since we already defined all of the necessary quantities, we could
proceed further and set up the background spacetime.

##  3 FLRW cosmology

In order to study the evolution of our universe, it will be handful to assume
that background spacetime is isotropic and homogeneous, namely is Friedmann-
Lemaitre-Robertson-Walker (FLRW) spacetime (we assume that lapse function is
unitary):

$ds^{2}=-dt^{2}+\sum_{i,j}a^{2}(t)dx^{i}dx^{j}.$ (15)

Here $a(t)$ is the scale factor of the universe, it is a fundamental quantity
that defines the evolution of the universe from its beginning. Consequently,
with assumption (15), the non-metricity scalar is written as follows (Caruana
et al., 2020)

$Q=6H^{2},$ (16)

Where $H=\dot{a}/a$ is the well-known Hubble parameter and dot over some
quantity signifies the first order temporal derivative. Finally, one could
evaluate the FLRW field equations of the $f(Q)$ theory:

$\displaystyle
3H^{2}=\kappa^{2}\left(\rho_{\text{m}}+\rho_{\text{eff}}\right)\,,$ (17)
$\displaystyle
3H^{2}+2\dot{H}=-\kappa^{2}\left(p_{\text{m}}+p_{\text{eff}}\right)\,.$ (18)

Here $\kappa^{2}=8\pi G$ is the Einstein gravitational constant squared,
$\rho_{\text{m}}$ and $p_{\text{m}}$ are, respectively, matter energy density
and isotropic pressure. Moreover, in the equation above $\rho_{\text{eff}}$
and $p_{\text{eff}}$ are effective energy density and pressure that define the
contribution of $f(Q)$ gravity to the field equations. For modified STEGR,
fields equations with plugged exact forms of effective quantities read:

$3H^{2}=\frac{\kappa^{2}}{2f_{Q}}\bigg{(}\rho_{m}+\frac{f}{2}\bigg{)},$ (19)
$\left(12H^{2}f_{QQ}+f_{Q}\right)\dot{H}=-\frac{k^{2}}{2}\left(\rho_{m}+p_{m}\right).$
(20)

are effective DE energy density and pressure that define the contribution of
the $f(Q)$ cosmology with

$f_{Q}=\frac{\partial f(Q)}{\partial Q},\quad
f_{QQ}=\frac{\partial^{2}f(Q)}{\partial Q^{2}}$ (21)

The energy-momentum tensor of the cosmological fluid which is given by

$T_{\mu\nu}=\left(\rho+p\right)u_{\mu}u_{\nu}+pg_{\mu\nu},$ (22)

which leads to conservation equation as $\dot{\rho}+3H\left(\rho+p\right)=0$.
In symmetric teleparallel gravity and its extensions, the conservation law
$T^{\mu}_{\,\,\,\nu;\mu}=0$ holds for the matter energy-momentum tensor. The
$T^{\mu}_{\,\,\,\nu;\mu}=0$ holds through (14) for the connection (Jiménez et
al., 2018b; Dimakis et al., 2022; Harko et al., 2018).

###  3.1 Exponential $f(Q)$ gravity

This paper is particularly aimed at the investigation of one $f(Q)$ gravity
model - namely modified exponential $f(Q)$ gravity (which is built from the
linear and exponential terms respectively). In $f(Q)$ theory, numerous cosmic
possibilities have been examined using various exponential models, notably
inflationary cosmology, BBN constraints, and dynamic system analysis.(Harko et
al., 2018; Anagnostopoulos et al., 2023; Khyllep et al., 2023). For that kind
of gravity, $f(Q)$ function reads (we adapt the work of (Linder, 2009, 2010)
for modified STEGR):

$f(Q)=Q+\alpha Q_{0}(1-e^{-\beta\sqrt{Q/Q_{0}}})$ (23)

Where $\alpha$, $\beta$ are free MOG parameters, namely additional degrees of
freedom. We can reduce the number of d.o.f by matching first Friedmann
equation (19) at the present time (i.e. assuming that $z=0$):

$\alpha=-\frac{e^{\beta}(-1+\Omega_{m0}+\Omega_{r0})}{-1+e^{\beta}-\beta}$
(24)

Thus, the complexity of this form is just one step more than the standard
$\Lambda$CDM. The exponential modified gravity could satisfy all stability and
validity and not cross the pahntom divide line (Arora & Sahoo, 2022). In order
to solve the field equations and obtain the numerical form of Hubble
parameter, the definition of $\dot{H}$ is

$\dot{H}=aH\frac{dH}{da}.$ (25)

We will solve the aforementioned field equation numerically, as we already
stated, with Mathematica numerical ODE solver NDSolve. Initial conditions at
the vanishing redshift for $\dot{H}$ could be therefore set up (as a
cosmographical quantity)(Mandal et al., 2020b):

$\dot{H}_{0}=-H_{0}^{2}(1+q_{0}).$ (26)

Where $q_{0}$ is the current deceleration parameter, we fix it to
$q_{0}=-0.55$ (Reid et al., 2019). Additionally, for MCMC training, as a
truths we assume that the present value of the Hubble parameter is
$H_{0}=69\,\mathrm{km/s/Mpc}$ and that matter mass fraction at the present
time is $\Omega_{m0}=0.315\pm 0.007$, following the observational constraints
of Planck2018 (Planck Collaboration et al., 2020a).

##  4 MCMC constraints

In this section, we want to constrain our $f(Q)$ gravity model via
observational datasets. To explore the parameter space, we will be using the
Markov Chain Monte Carlo (MCMC) methodology and Python package emcee (Foreman-
Mackey et al., 2013).

###  4.1 Observational Hubble Data

Observational Hubble Data is one of the most popular and plausible tests of
the universe expansion history beyond GR and $\Lambda$CDM. The OHD sample is
mainly obtained from the differential age of galaxies method (or just DAG) (Yu
et al., 2018; Moresco, 2015). In this method, Hubble rate is usually obtained
from the formula below

$H(z)=\frac{-1}{1+z}\frac{dz}{dt}$ (27)

In the current article, we will primarily use the OHD points derived from the
so-called Cosmic Chronometers (CC) i.e., the massive and passively evolving
galaxies. Using Cosmic Chronometers, ratio $dz/dt$ could be derived from the
$\Delta z/\Delta t$, where $\Delta z$ is the redshift separation in the
galaxies sample and could be easily determined through precise and accurate
spectroscopy. On the other hand, derivation of the $\Delta t$ is much more
challenging and requires some standard clocks. For that purpose, we could use
massive, passively evolving, and old stellar populations that are present
across a wide range of redshifts and therefore could be considered as cosmic
chronometers. To determine the priors and likelihood functions (which are
necessary), we used the $H(z)$ dataset. To constrain our modified gravity
model, we introduce the chi-squared function below

$\chi^{2}_{CC}=\sum_{i=1}^{N_{H}}\bigg{[}\frac{H^{\mathrm{th}}_{i}(p_{1},p_{2},...,p_{n},z_{i})-H^{\mathrm{obs}}_{i}(z_{i})}{\sigma_{H(z_{i})}}\bigg{]}^{2}$
(28)

The likelihood function that we will be using for MCMC sampling have its usual
exponential form

$\mathcal{L}=\exp(-\chi^{2}/2)$ (29)

###  4.2 Pantheon SN Ia Sample

We also used the Pantheon dataset to constrain our modified gravity with dark
energy, which consists of data obtained from the 1048 Ia supernovae
(discovered by the PANSTARRS DR1 (PS1) Medium Deep Survey, Low $z$, SNLS, SDSS
and HST (Scolnic et al., 2018a; Chang et al., 2019)). For this dataset,
redshift varies from $z=0.01$ to $z=2.26$. The corresponding chi-squared
function reads:

$\chi^{2}_{SN}(p_{1},p_{2},...,p_{n})=\sum^{N_{SN}}_{i,j=1}\frac{\Delta\mu_{i}}{\sigma_{\mu(z_{i})}}$
(30)

Where

$\Delta\mu_{i}=\mu^{th}(p_{1},p_{2},...,p_{n})-\mu_{i}^{obs}$ (31)

And distance moduli is (Arora & Sahoo, 2020)

$\mu^{th}=5\log_{10}D_{L}(z)+\mu_{0},\quad\mu_{0}=5\log_{10}\frac{H_{0}^{-1}}{\mathrm{Mpc}}+25$
(32)
$D_{L}(z)=\frac{c(1+z)}{H_{0}}S_{K}\bigg{(}H_{0}\int^{z}_{0}\frac{d\overline{z}}{H(\overline{z})}\bigg{)}$
(33)

Here function $S_{K}(x)$ is just

$S_{K}(x)=\begin{cases}\sinh(x\sqrt{\Omega_{K}})/\Omega_{K},\quad\Omega_{K}>0\\\
x,\quad\quad\quad\quad\quad\quad\quad\;\;\;\Omega_{K}=0\\\
\sin(x\sqrt{|\Omega_{K}|})/|\Omega_{K}|,\quad\Omega_{K}<0\end{cases}\,.$ (34)

It is known that out universe is spatially flat, and therefore $\Omega_{K}=0$.
The nuisance parameters in the Tripp formula (Tripp, 1998)
$\mu=m_{B}-M_{B}+\alpha x_{1}-\beta c+\Delta_{M}+\Delta_{B}$ were retrieved
using the novel method known as BEAMS with Bias Correction (BBC) (Kessler &
Scolnic, 2017), and the observed distance modulus is now equal to the
difference between the corrected apparent magnitude $M_{B}$ and the absolute
magnitude $m_{B}$ $\left(\mu=m_{B}-M_{B}\right)$. Additionally, one can define
the chi-squared function in terms of covariance matrix as follows (Deng & Wei,
2018):

$\chi^{2}_{\mathrm{SN}}=\Delta\boldsymbol{\mu}^{T}\mathbf{C}^{-1}\Delta\boldsymbol{\mu}$
(35)

Where covariance matrix consists of systematic and statistical uncertainties
respectively (Conley et al., 2011):

$\mathbf{C}=\mathbf{D}_{\mathrm{stat}}+\mathbf{C}_{\mathrm{sys}}$ (36)

In the current work we assume that diagonal matrix of statistical
uncertainties looks like
$\mathbf{D}_{\mathrm{stat},ii}=\sigma^{2}_{\mu(z_{i})}$. Besides, systematic
uncertainties are derived using the Bias Corrections (BBC) method, introduced
and developed in (Scolnic et al., 2018b):

$\mathbf{C}_{ij,\mathrm{sys}}=\sum^{K}_{k=1}\bigg{(}\frac{\partial\mu^{obs}_{i}}{\partial
S_{k}}\bigg{)}\bigg{(}\frac{\partial\mu^{obs}_{j}}{\partial
S_{k}}\bigg{)}\sigma^{2}_{S_{k}}$ (37)

Indexes $\\{i,j\\}$ denote the redshift bins for distance modulus, $S_{k}$
here denotes the magnitude of systematic error, $\sigma_{S_{k}}$ is
respectively it’s standard deviation uncertainty.

Figure 1: MCMC best fits from OHD, Pantheon and BAO datasets and joint
distribution for exponential $f(Q)$ model

###  4.3 Baryon Acoustic Oscillations

Finally, we use Baryon Acoustic Oscillations (BAOs) to constrain our modified
gravity model. BAOs arise in the early times of universe evolution. In earlier
times, fermions and photons are strongly connected to each other due to the
Thompson scattering. This mixture of baryons and photons behaves like a single
fluid and can not gravitationally collapse. Moreover, this fluid oscillates
because of the huge photonic pressure. These oscillations are called BAOs. The
Characteristic scale of the BAO is defined by the so-called sound horizon
$r_{s}$, which is seen at the photon decoupling epoch with redshift $z_{*}$:

$r_{s}=\frac{c}{\sqrt{3}}\int^{\frac{1}{1+z_{*}}}_{0}\frac{da}{a^{2}H\sqrt{1+(3\Omega_{b0}/4\Omega_{\gamma
0})a}}$ (38)

Here $\Omega_{b0}$ is known as the baryon mass density at present ($z=0$) and
$\Omega_{\gamma 0}$ is, respectively, photon mass density at present. Also, as
it was noticed, angular diameter distance is derived directly from the BAO
sound horizon. In this work, to constrain our MOG model with BAO, we will use
the observational datasets with $d_{A}(z_{*})/D_{V}(z_{BAO})$ data. Here we
consider that $d_{A}(z_{*})$ is the angular diameter distance in the comoving
coordinates:

$d_{A}(z)=\int^{z}_{0}\frac{dz^{\prime}}{H(z^{\prime})}$ (39)

And $D_{V}(z_{BAO})$ is the dilation scale:

$D_{V}(z)=(d_{A}(z)^{2}z/H(z))^{1/3}$ (40)

Finally, we also consider that photon decoupling epoch arise at the redshift
(Planck Collaboration et al., 2016):

$z_{*}=1048[1+0.00124(\Omega_{b}h^{2})^{-0.738}][1+g_{1}(\Omega_{m}h^{2})^{g_{2}}]$
(41)

Where,

$g_{1}=\frac{0.0783(\Omega_{b}h^{2})^{-0.238}}{1+39.5(\Omega_{b}h^{2})^{-0.763}}$
(42) $g_{2}=\frac{0.560}{1+21.1(\Omega_{b}h^{2})^{-1.81}}$ (43)

This dataset was gathered from the works of (Blake et al., 2011; Percival et
al., 2010; Jarosik et al., 2011; Eisenstein et al., 2005; Giostri et al.,
2012).

Consequently, to perform the MCMC sampling, we need to define the chi squared
function for our BAO dataset:

$\chi^{2}_{BAO}=X^{T}C^{-1}X$ (44)

Where $X$ is the matrix of form (Giostri et al., 2012):

$X=\left(\begin{array}[]{c}\frac{d_{A}(z_{*})}{D_{V}(0.106)}-30.95\\\
\frac{d_{A}(z_{*})}{D_{V}(0.2)}-17.55\\\
\frac{d_{A}(z_{*})}{D_{V}(0.35)}-10.11\\\
\frac{d_{A}(z_{*})}{D_{V}(0.44)}-8.44\\\
\frac{d_{A}(z_{*})}{D_{V}(0.6)}-6.69\\\
\frac{d_{A}(z_{*})}{D_{V}(0.73)}-5.45\end{array}\right)\,$ (45)

We also performed the joint analysis form the combined $OHD+SN+BAO$ by
minimizing $\chi^{2}_{OHD}+\chi^{2}_{SN}+\chi^{2}_{BAO}$.

Datasets $H_{0}$ $\Omega_{m0}$ $\beta$ Hubble (OHD) $66.9\pm 3.3$
$0.320^{+0.055}_{-0.070}$ $4.3\pm 1.9$ OHD+SNa $68.9\pm 1.7$
$0.290^{+0.028}_{-0.020}$ $5.3^{+1.8}_{-1.0}$ OHD+SNa+BAO $68.9\pm 1.6$
$0.292\pm 0.016$ $5.6\pm 1.25$ $\rm Models$ $\chi_{\rm min}^{2}$ $\rm AIC$
$\rm BIC$ $\Lambda$CDM 58.700 67.248 76.127 $f(Q)$ 57.616 68.499 79.137

Table 1: Best-fit values of model parameters and statistical analysis

The results are, therefore, numerically derived from MCMC trained on OHD,
Pantheon, BAO and joint datasets. Besides, results are placed on the Table (1)
above for model free parameters $H_{0}$, $\beta$ and $\Omega_{m0}$.
Furthermore, the $1-\sigma$ and $2-\sigma$ likelihood contours for the
possible subsets of parameter space are presented in Fig. 1.

###  4.4 Statistical evaluation

To evaluate the success of our MCMC analysis, one should perform the
statistical evaluation using the so-called Akaike Information Criterion (AIC)
and Bayesian Information Criterion (BIC). The first quantity, namely AIC can
be expressed as follows (Akaike, 1974):

$\mathrm{AIC}=\chi^{2}_{\mathrm{min}}+2d$ (46)

With $d$ being the number of free parameters in a chosen model. To compare our
results with the well-known fiducial $\Lambda$CDM model, we are going to use
the AIC difference between our modified gravity model and fiducial cosmology
$\Delta\mathrm{AIC}=|\mathrm{AIC}_{\Lambda\mathrm{CDM}}-\mathrm{AIC}_{\mathrm{MOG}}|$.
In that case, if $\Delta\mathrm{AIC}<2$, there is a strong evidence in favor
of MOG model, while for $4<\Delta\mathrm{AIC}\leq 7$ there is a little
evidence if favor of MOG model of our consideration. Finally, for the case
with $\Delta\mathrm{AIC}>10$ there is practically no evidence in favor of MOG
(Liddle, 2007). In addition, BIC is defined through the relation, written down
below:

$\mathrm{BIC}=\chi^{2}_{\mathrm{min}}+d\ln N$ (47)

For that case, $N$ is the number of data points being used for MCMC. For BIC,
if $\Delta\mathrm{BIC}<2$, there is no strong evidence against chosen model
that deviate from $\Lambda$CDM, if $2\leq\Delta\mathrm{BIC}<6$ there is
evidence against the MOG model and finally for $\Delta\mathrm{BIC}>6$ there is
strong evidence against MOG model. We therefore store the
$\chi^{2}_{\mathrm{min}}$/AIC/BIC data for modified gravity model of our
consideration in the Table (1). As we see, $\Delta\rm AIC=1.25$ and $\Delta\rm
BIC=3.01$ so that our model can very precisely mimic $\Lambda$CDM one.

##  5 Validity of cosmological constraints

Figure 2: Dimensionless mass density for matter and effective dark energy
within exponential $f(Q)$ gravitation Figure 3: Hubble parameter $H(z)$,
deceleration parameter $q(z)$ and distance modulus for exponential $f(Q)$
gravity with best fit values from MCMC used. For comparison, we as well show
the fiducial $\Lambda$CDM results

In order to check the validity of the aforementioned applied cosmological
constraints (such as Pantheon, BAO, or OHD), we will probe the behavior of
some quantities, such as deceleration parameter $q(z)$ or statefinder pair.
Furthermore, one could define the so-called deceleration parameter:

$q=-\frac{\dot{H}}{H^{2}}-1.$ (48)

From the above equation, one could easily notice that both deceleration and
Hubble parameters are related to each other by higher order temporal
derivatives of scale factor $a$. Consequently, to differentiate our model from
other numerous MOG, DE/DM models, one could present the pair of parameters,
the so-called statefinder pair (Sahni et al., 2003; Alam et al., 2003; Pasqua
et al., 2015; Xu et al., 2018):

$r=\frac{\dddot{a}}{aH^{3}},$ (49) $s=\frac{r-1}{3(q-1/2)}.$ (50)

For the sake of simplicity, we could redefine the statefinder pair $\\{r,s\\}$
fully in terms of deceleration parameter:

$r(z)=q(z)(1+2q(z))+q^{\prime}(z)(1+z),$ (51)
$s(z)=\frac{r(z)-1}{3(q(z)-1/2)}.$ (52)

We are going to construct the phase plane $r(z)-s(z)$, in which different
points correspond to the various universe states, such that:

  * •

$\Lambda$CDM corresponds to $(s=0,r=1)$,

  * •

Chaplygin Gas (CG) corresponds to $(s<0,r>1)$,

  * •

SCDM corresponds to $(r=1,q=0.5)$,

  * •

Quintessence corresponds to $(s>0,r<1)$.

Consequently, we plot both statefinder parameter phase portraits, deceleration
parameter, and additionally $H(z)$ as probes of model validity in the
cosmological sense in Figs. (3) and (4). Statefinder diagnostics and $q(z)$
were performed only for the joint dataset, since other datasets shows the
similar behavior as the joint solution. Remarkably, a transition from
deceleration to acceleration phase on the third plot of the aforementioned
figure is seen. A valid interval for $q_{0}$ is marked as a gray area. As
stated already, one may check the universe evolutionary scenario using
statefinder pairs $\\{r,s\\}$ and $\\{r,q\\}$. From the $r-s$ plane of our
model, one could observe that the initial universe was filled with
quintessence, then passed the $\Lambda$CDM phase and is currently reverting
towards the quintessence scenario. On the other hand, in the $r-q$ plane, it
is evident that the universe once also passed through the $\Lambda$CDM phase.
However, now our space-time is generally filled with quintessential fluid, it
is expected that the future universe will eventually turn to the de-Sitter
state (when $\Lambda$ term will fully dominate). The point on quintessential
fluid also coincides with MCMC observational constraints.

The very last probe of cosmological validity is the well known $Om(z)$
diagnostics, firstly presented in the paper (Sahni et al., 2008), where
$Om(z)$ is defined through equation (Pan et al., 2018; Harko et al., 2022):

$Om(z)=\frac{E^{2}(z)-1}{(1+z)^{3}-1}$ (53)

This parameter was derived to distinguish $\Lambda$CDM from other, more
complicated cosmological models. In the equation (53), it is a handful to
define $E^{2}(z)=H^{2}(z)/H_{0}^{2}$, which is exactly the Hubble flow,
dimensionless quantity normalized by the current Hubble parameter value.
$Om(z)$ parameter have a constant value for the $\Lambda$CDM model, which is
same as the current matter mass density $\Omega_{m0}$.

Figure 4: Statefinder pairs and $Om(z)$ function for exponential $f(Q)$
gravity, $\Lambda$CDM and $\omega$ varying $\Lambda$CDM cosmologies

Consequently, we place numerical solution of $Om(z)$ function for $f(Q)$ model
in the Fig.(4). For the sake of comparison, we as well plot $Om(z)$ solutions
within classical $\Lambda$CDM model and within $\omega$ varying $\Lambda$CDM
cosmologies. As one could easily notice, for our $f(Q)$ model, $Om(z)$ shows
only $Om(z)<\Omega_{m0}$ behavior in the distant past, which could lead to the
presence of phantom fluid (for more information on the subject, see (Mostaghel
et al., 2017)).

Figure 5: N-body simulations snapshot (CDM over density) for $f(Q)$ gravity
with best fit MCMC values on different redshifts Figure 6: SPH simulation
snapshots of $f(Q)$ gravity for gas temperature $T$ and mean molecular weight
$\mu$

However, at $z\approx 2$, our model transits $\Lambda$CDM and has constantly
growing trend, therefore in the near past and present times, quintessential
fluid appears, which converges well with the statefinder diagnostics and MCMC.

Finally, we also analyse both matter and effective dark energy mass densities
for our model of modified gravitation to conclude on its validity.
Corresponding results are plotted on the Figure (2). One can easily notice,
both $\Omega_{m0}\land\Omega_{\Lambda 0}\in[0,1]$ and their sum always
converges to unity, epoch of the equality appears at redshift $z\approx 0.35$,
which is very near to the $\Lambda$CDM estimate.

##  6 N-body simulations of LSS with small $\mathrm{L_{\mathrm{Box}}}$

As we already remarked previously, the main purpose of this paper is to
perform N-body simulations of the comoving box that contain DM+baryonic matter
and dark energy in exponential $f(Q)$ gravitation and compare our results with
the Large Scale Structure of concordance $\Lambda$CDM cosmology. For that aim,
we will use the publicly available code ME-GADGET, a modification of the well-
known hydrodynamical N-body code GADGET2. It was modified for generality so
that one can perform simulations for practically any cosmological model. The
code above was described in the pioneering works of (An et al., 2019; Zhang et
al., 2019), whereas the tests are provided in (Zhang et al., 2018b). This code
as an input needs tables with Hubble flow $H/H_{0}$ and the deviation of
effective gravitational constant from the Newtonian one
$G_{\mathrm{eff}}/G_{N}$ (in some models of modified gravity, namely screened
ones, such deviation exists only up to some scale $k_{\mathrm{screen}}$
because of the so-called fifth force). The effective gravitational constant
exact form was found in the paper (Jiménez et al., 2020a):

$G_{\mathrm{eff}}=\frac{G_{N}}{f_{Q}}$ (54)

Equation above is being numerically solved assuming appropriate best fit
values for free parameters of our model. As one could easily notice, at the
very early time (high-$z$ epochs), $f(Q)$ gravity has Newtonian-like
gravitational constant and then, at approximately $a\approx 0.1$,
$G_{\mathrm{eff}}$ is being separated from $G_{\mathrm{N}}$ for our model.

Since we already defined needed inputs for ME-GADGET code, we could proceed
further to fine-tuning our simulation setup.

###  6.1 Simulation setup

One needs to define various parameters to produce the simulations and initial
conditions, based on the second order Lagrangian Perturbation Theory (namely,
2LPT). We want to obtain the mid-resolution simulations, therefore, particle
number is $N=512^{3}$ and mesh size is respectively $N_{\mathrm{mesh}}=2\times
512^{3}$ as well. The simulation box has periodic vacuum boundary conditions
and sides with length $10\mathrm{Mpc/h}$. Initial conditions were produced
with the Simp2LPTic code (see GitHub repository
https://github.com/liambx/Simp2LPTic), glass files (pre-initial conditions)
were generated with the use of ccvt-preic (check
https://github.com/liaoshong/ccvt-preic). We assumed that glass tile fraction
is unitary. Moreover, cosmological parameters were borrowed from our MCMC
constraints, discussed earlier: $h=100H_{0}=0.689\pm 0.016$ (so-called
"little-h"), $\Omega_{m0}=0.292\pm 0.016$, leading to $\Omega_{\Lambda
0}=0.708$, if one will not take into account radiation and massive neutrino
species. On the other hand, baryon mass density equals to $\Omega_{b}=0.0493$
(relation between total matter density and baryon one decides how much gas
particles are present in the simulation). Moreover, matter power spectrum
amplitude at $k=8\mathrm{Mpc/h}$ is assumed to be $\sigma_{8}=0.811\pm 0.006$
and initial power spectrum is linear, constructed from the Eisenstein & Hu
transfer function (Eisenstein & Hu, 1998) (power spectrum were constructed
using code CAMB, see (Lewis & Challinor, 2011)). Initial conditions are being
generated at the redshift $z=10$ and spectrum index of scalar perturbations is
$n_{s}=0.9649\pm 0.0042$ (Planck Collaboration et al., 2020b).

###  6.2 Results

In the current subsection, we are going to discuss the main results, obtained
from the N-body simulations of the Large Scale Structure of the Universe.

Firstly, we demonstrate the spatial slices of CDM over density
$\delta_{\mathrm{CDM}}=\overline{\rho}_{\mathrm{CDM}}/\rho_{\mathrm{CDM}}$ for
our $f(Q)$ gravity model with different values of redshift $z$ on the Figure
(5).

In addition to the over density measurements, we also show the temperature of
gas $T$ that arise from Smoothed Particle Hydrodynamics (SPH) and mean
molecular weight $\mu=\overline{m}/m_{\mathrm{HI}}$, which defines the
relation between mean particle mass and neutral hydrogen particle mass on the
Figure (7) respectively. As one can easily notice, the DM walls are
represented by smaller value of mean molecular weight. Besides, temperature
maps show the well-known hot "bubbles" within the Inter-Galactic Medium (IGM)
that are formed due to impinging galactic winds.

Now we are going to investigate the matter power spectrum for our model. In
comparison, we are going to use fiducial $\Lambda$CDM cosmology power
spectrum, generated with the use of CAMB code (Lewis & Challinor, 2011; Lewis
& Bridle, 2002; Lewis et al., 2000; Howlett et al., 2012)111Documentation for
this code is stored in camb.readthedocs.io. In order to extract $P(k)$ for
some value of redshift within our N-body framework, we used Python-based code
Pylians3 (Villaescusa-Navarro, 2018)222For installation procedure and full
documentation, refer to the pylians3.readthedocs.io.

We consequently compare the matter power spectrum on the Figure (7)
with/without Redshift-Space Distortions (RSDs) directed along both $X$, $Y$
and $Z$ axes. As we noticed during numerical analysis, up to some $k$ near
$k_{\mathrm{Box}}$ limit for our simulation, $P(k)$ spectrum in Fourier space
does reconstruct non-linear matter power spectrum, given by CAMB, while
Redshift-Space Distorted (RSD) one behave like linear matter power spectrum,
as expected. Also, it is worth to notice that difference between RSD and
regular matter power spectrum is bigger for CDM+Gas case. Finally, effect of
RSDs in our simulations is almost isotropic, so that $\Delta(\mathrm{RSD})$
differs only by few percents with the change of RSD direction axis.

Figure 7: Matter power spectrum with/without RSDs for $f(Q)$ gravity vs. CAMB
linear/non-linear $P(k)$ for $\Lambda$CDM. Dashed N-body $P(k)$ represent the
CDM-only power spectrum, while solid line represent CDM+Gas $P(k)$. Error bars
represent Ly $\alpha$ forest observations on high $z$ Figure 8: Monopole
redshift-space distorted two point correlation function with
$L_{\mathrm{Box}}=10/100h^{-1}$Mpc for $f(Q)$ log-like modified gravity. In
relation we plot Quijote simulation correlation function for Planck fiducial
cosmology with Gpc wide box. Also, for each case we display the scale, at
which BAO bump occurs

###  6.3 Halo mass function

Now we correspondingly derive the well-known halo mass function (further -
just HMF), that obviously define the number of halos at a certain mass.
Firstly, we built the halo catalogue for all of our snapshots with the use of
halo/subhalo structure finder, namely ROCKSTAR (for more information on the
subject, refer to the documentation paper (Behroozi et al., 2013)).
Afterwards, we built the binned halo mass function, which is based on the
values of $M_{\mathrm{200c}}$ (the mass of enclosed halo volume with energy
density 200 times bigger that critical density of the universe
$\rho_{\mathrm{cr}}$). Results are respectively plotted on the Figure (9) with
the added Seth-Tormen theoretical prediction for halo mass function, based on
Planck2018 fiducial cosmology and CAMB power spectrum at the $z=0$. Seth-
Tormen HMF were computed using the python package pynbody (Pontzen et al.,
2013).

Figure 9: Halo mass function for $f(Q)$ gravitation with
$L_{\mathrm{Box}}=10\mathrm{\;Mpc}/h$ and theoretical prediction for HMF by
Seth-Tormen

From the Figure (9) shown above one could easily notice that generally our
prediction for halo mass function from the MOG N-body simulation shows values
of $n$ that approximately coincide with the ones that theoretically predicted
by the Seth-Tormen HMF.

##  7 LSS with large $L_{\mathrm{Box}}$: comparison

As we previously mentioned, now we are going to perform an analysis of N-body
simulation for bigger simulation box size, namely with
$L_{\mathrm{box}}=100h^{-1}$Mpc. In that case, we only differ force resolution
($\epsilon=3.9$kpc), other cosmological parameters are assumed to be the same.
At first, we as usual plot the CDM over density field for vanishing redshift
at the Figure (10).

In addition, we also plot the matter power spectrum for CDM, CDM+baryons on
the Figure (11). As an obvious consequence of a larger box size, one can
notice that maximum wave number $k$ grew to $k_{\mathrm{max}}\approx
20\;h/$Mpc. Even at such big scales, our matter power spectrum, derived from
the corresponding N-body simulation converge with theoretical prediction from
CAMB code with up to sub-percent accuracy. As we noticed previously for small
simulation box, axis of redshift-space distorsions had a vary small impact on
the matter power spectrum. This statement holds for large $L_{\mathrm{box}}$
as well.

In the previous section we discussed the halo mass function for the case with
smaller simulation box. Now we can discuss the same matter but for the larger
$L_{\rm Box}$. As it appears, HMF extracted from the simulation replicates the
Seth-Tormen HMF almost perfectly up to $M\approx 10^{14}M_{\odot}$ (see Figure
(9)). However, at bigger halo masses, our simulation HMF slightly differs from
the theoretical prediction, which is usually observed in N-body simulations.

Figure 10: Present day snapshot of CDM over density for
$L_{\mathrm{box}}=100h^{-1}$Mpc run Figure 11: Matter power spectrum
with/without RSDs for $f(Q)$ gravity vs. CAMB linear/non-linear $P(k)$ for
$\Lambda$CDM. Dashed N-body $P(k)$ represent the CDM-only power spectrum,
while solid line represent CDM+Gas $P(k)$. Error bars represent Ly $\alpha$
forest observations on high $z$. For this case, we have assumed large
simulation box size of 100$h^{-1}$Mpc Figure 12: Relation of reduced matter
power spectrum $Q(k_{1},k_{2},k_{3})$ for both small and large simulation
volumes with $1\sigma$ deviation, varying $k_{1}=2k_{2}$ Figure 13: Two-
dimensional matter power spectrum for small and large simulations within
$f(Q)$ gravitation with/without RSDs

###  7.1 Reduced bispectrum from 3PCF

Finally, we also are going to introduce the so-called reduced bispectrum,
which is derived from the regular bispectrum and matter power spectrum via the
following relation, written down below:

$Q=\frac{B}{P_{1}P_{2}+P_{2}P_{3}+P_{1}P_{3}}$ (55)

Where $P_{i}=P_{m}(k_{i})$. We consequently plot the relation between
bispectrum of large and small cosmological volumes on the Figure (12). It is
easy to notice that for smaller wave number ($k_{1}=5$), relation between
$Q(k_{1},k_{2},k_{3})$ for both cases has a mean value of $\approx 1.9$ for
all bins of angle $\theta$ (where it’s maximum value is $\theta=\pi$, which is
the angle between two sides of triangle $k_{1}$ and $k_{2}$). On the other
hand, for relatively big $k_{1}$ (in our case, it’s $k_{1}=6h/\rm Mpc$),
deviation of large $L_{\mathrm{box}}$ reduced bispectrum from the small one is
smaller because of the the fact that wave number span is shifted towards
bigger values, while in the first case $k_{1}=5h/\rm Mpc$ were on the box size
limit for smaller simulation, which distorted the results and caused the
deviation to grow.

###  7.2 2PCF for $f(Q)$ gravitation

In addition to the matter power spectrum/reduced bispectrum and halo mass
function, we as well derive the two point correlation function (further - just
2PCF) in a real space for CDM halos. Generally, 2PCF is defined as follows:

$\xi(|\mathbf{x}_{1}-\mathbf{x}_{2}|)=\langle\delta_{m}(\mathbf{x}_{1})\delta_{m}(\mathbf{x}_{2})\rangle$
(56)

Where $\mathbf{x}_{i}$ is three-dimensional position of an $i$-th CDM halo and
$\delta_{m}$ is CDM over density parameter. We show the monopole redshift-
space distorted two point correlation functions for both large and small
simulations boxes on the Figure (8), where we added Quijote simulations
(Villaescusa-Navarro et al., 2020) 2PCF, that admits Planck fiducial
cosmology. For the sake of completeness, we additionally marked BAO bumps for
each case with color-coded dotted line. It is obvious, that in the case of
small simulation box size, permitted range of $R$ is very small (up to
$R\approx 2\times 10^{0}h^{-1}$Mpc) and because of the small box size,
correlation function is under sampled and does not coincide with Quijote one.
On the other hand, for $L_{\rm Box}=100h^{-1}$Mpc simulation, correlation
function corresponds to the Quijote one with sub-percent accuracy for range
$R\in[2\times 10^{0},10^{1}]$. Now, we can proceed to the latest topic of our
consideration, namely two-dimensional power spectra.

###  7.3 2D matter power spectrum

We plot the two-dimensional matter power spectrum for small and large boxes on
the Figure (13) with/without the presence of redshift-space distortions. As
one can easily notice, on the plots with RSDs, the so-called "Finger of God"
effect is observed, which arise because of the large scatter of galaxies
recessional velocities at the small scales. Also, it is worth to inform that
2D matter power spectra for both box sizes are very alike. Now, since we
discussed all of the topics for both simulation volumes within the modified
theory of gravitation, we can proceed to the concluding remarks on the key
findings within our study.

##  8 Conclusions

One can describe gravity using several geometric bases. The STGR, which
attributes gravity to the nonmetricity tensor, has recently drawn much
attention. A fascinating method for studying modified gravity is $f(Q)$
gravity, an extension of STGR. This study examined large scale structure
formation observables using N-body simulations of $f(Q)$ gravitation for the
first time to assess the theory’s validity to cosmological context.
Simulations were run with the use of ME-GADGET code, modification of the
widely known GADGET-2 code for two simulation boxes, namely $L_{\rm
Box}=10h^{-1}$Mpc and $L_{\rm Box}=100h^{-1}$Mpc to decide on an optimal box
size and compare the results for both simulation volumes.

We first performed Markov Chain Monte Carlo (MCMC) analysis for our
exponential $f(Q)$ model to obtain best-fit values of MOG free parameters in
Section (4). To test the fits provided by MCMC, we obtained theoretical
predictions for the dimensionless mass densities $\Omega_{m0}$ and
$\Omega_{\Lambda 0}$, Hubble parameter $H(z)$, deceleration parameter $q(z)$
and statefinder pair $\\{r,s\\}$, $Om(z)$ parameter, placing graphical results
on the Figures (2), (3) and (4) respectively. As we noticed, Hubble parameter
respected low redshift observations as well as deceleration parameter provided
correct values of $q_{0}$ and transitional redshift within the constrained
range. Moreover, statefinder diagnostics predict that initially universe was
in Quintessence phase, passed the $\Lambda$CDM state up to Quintessence again.
Finally, $Om(z)$ demonstrated that at the high-$z$ range, our universe was
filled with phantom fluid, passed $\Lambda$CDM EoS at $z\approx 2$ and now
again has a phantom equation of state. After theoretical predictions, we
started working on the N-body simulations whose primary findings corresponding
to the quantities of interest are as follows:

  * •

Three-dimensional matter power spectrum monopole $P_{k}$: this was the first
probe of a large-scale structure we studied in the present work. We plotted
non-linear matter power spectra (with/without RSDs) for both small and large
simulation volumes on the Figures (7) and (11) respectively, where we plotted
CAMB linear/non-linear fiducial power spectra and observational data from
Ly-$\alpha$ forest for the sake of comparison. One can notice that within the
permitted range of wave number $k$ (limited by mean inter particle separation
and simulation box length), non-linear matter power-spectra from small/large
N-body simulations coincide with the CAMB one. However, for $L_{\rm
Box}=10h^{-1}$Mpc, non-linear $P_{k}$ coincide with linear CAMB prediction too
early because of the small box size.

  * •

Halo Mass Function: the second significant quantity that can solely conclude
the validity of a simulation. We place Seth-Tormen’s theoretical HMF and the
ones extracted from our N-body simulations in Figure (9). As we found, our
small box size can cannot provide sufficient halo masses and reproduce viable
halo mass function at all mass ranges up to the resolution limit, while large
simulation follows Seth-Tormen prediction very precisely within the large span
of halo masses $\log_{10}M/M_{\odot}\in[10,14]$, but gets slightly smaller
than theoretical prediction for higher masses.

  * •

Two-Point Correlation Function monopole $\xi_{0}(r)$: we as well investigate
the redshift-space distorted correlation function monopoles in the Figure (8),
where Quijote simulations correlation function is plotted to compare our
results to fiducial cosmology. As remarked during numerical analysis, small
box simulation fails to predict correct CDM halo correlations. On the other
hand, in the range, $R\in[2\times 10^{0},10^{1}]h^{-1}$Mpc large box
simulation precisely reconstructs Quijote data.

  * •

Reduced bispectrum $Q(k_{1},k_{2},k_{3})$: for the reduced bispectra case, we
plotted the relation $(Q_{\rm Large}-Q_{\rm Small})/Q_{\rm Small}$ with
different $k_{1}$ values (which acts as a triangle side length) on the Figure
(12). We observed that this relation is generally close to $\approx 1.5$
across all bins of the angle between $k_{1}$ and $k_{2}$ (namely $\theta$) if
one will assume value of $k_{1}$ that is not on the resolution limit for both
cases (while it is worth to notice that we only adopted the case, where
$k_{2}=2k_{1}$).

  * •

Two-dimensional matter power spectrum $P_{m}(k_{\parallel},k_{\perp})$: This
is the last quantity extracted from our N-body simulations. We plotted 2D
power spectra for both sims in Figure (13) with/without redshift-space
distortions. From the plots, we noticed the so-called "Finger of God" effect
that is present in the RSD case because of the elongated positions of CDM
halos.

In conclusion, considering all the above points, the small simulation volume
experiment failed to recreate matter power spectrum and correlation function
correctly. However, second one, namely more extensive N-body simulation
provided both viable 3D/2D matter power spectrum, correlation function, and
halo mass functions and, therefore, we can consider exponential $f(Q)$ model
to be a viable substitution of fiducial $\Lambda$CDM cosmology, since it not
only satisfy many large scale structure constraints mentioned above, but also
provide correct distance modulus up to high redshift values, where
$\Lambda$CDM fails.

In the following papers, it will be interesting to investigate this modified
gravity model using hybrid N-body and SPH simulations that incorporate
supernova/AGN feedback, star and galaxy formation, jets, etc. using code GIZMO
that allows the use of tabulated Hubble parameter and effective gravitational
constant. It will, however, require a lot more computational resources (on the
scale of millions of CPU hours).

## Data Availability Statement

There are no new data associated with this article.

## Acknowledgements

Sokoliuk O. performed the work in frame of the "Mathematical modeling in
interdisciplinary research of processes and systems based on intelligent
supercomputer, grid and cloud technologies" program of the NAS of Ukraine. The
authors gratefully acknowledge the computing time provided on the high
performance computing facility, Sharanga, at the Birla Institute of Technology
and Science - Pilani, Hyderabad Campus. PKS acknowledges the Science and
Engineering Research Board, Department of Science and Technology, Government
of India for financial support to carry out the Research project No.:
CRG/2022/001847. SA acknowledges BITS-Pilani, Hyderabad campus for the
Institute fellowship. We are very much grateful to the honorable referee and
to the editor for the illuminating suggestions that have significantly
improved our work in terms of research quality, and presentation.

## References

  * Abedi et al. (2017) Abedi H., Wright M., Abbassi A. M., 2017, Physical Review D, 95, 064020
  * Akaike (1974) Akaike H., 1974, IEEE Transactions on Automatic Control, 19, 716
  * Alam et al. (2003) Alam U., Sahni V., Deep Saini T., Starobinsky A., 2003, Monthly Notices of the Royal Astronomical Society, 344, 1057
  * Albuquerque & Frusciante (2022) Albuquerque I. S., Frusciante N., 2022, Physics of the Dark Universe, 35, 100980
  * An et al. (2019) An R., Costa A. A., Xiao L., Zhang J., Wang B., 2019, MNRAS, 489, 297
  * Anagnostopoulos et al. (2021) Anagnostopoulos F. K., Basilakos S., Saridakis E. N., 2021, Physics Letters B, 822, 136634
  * Anagnostopoulos et al. (2023) Anagnostopoulos F. K., Gakis V., Saridakis E. N., Basilakos S., 2023, The European Physical Journal C, 83, 58
  * Arnold et al. (2015) Arnold C., Puchwein E., Springel V., 2015, Mon. Not. Roy. Astron. Soc., 448, 2275
  * Arnold et al. (2016) Arnold C., Springel V., Puchwein E., 2016, Mon. Not. Roy. Astron. Soc., 462, 1530
  * Arora & Sahoo (2020) Arora S., Sahoo P. K., 2020, Physica Scripta, 95, 095003
  * Arora & Sahoo (2022) Arora S., Sahoo P. K., 2022, Annalen der Physik, 534, 2200233
  * Atayde & Frusciante (2021) Atayde L., Frusciante N., 2021, Physical Review D, 104, 064052
  * Bajardi et al. (2020) Bajardi F., Vernieri D., Capozziello S., 2020, The European Physical Journal Plus, 135, 1
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, Astrophys. J., 762, 109
  * Blake et al. (2011) Blake C., et al., 2011, Monthly Notices of the Royal Astronomical Society, 418, 1707
  * Capozziello & D’Agostino (2022) Capozziello S., D’Agostino R., 2022
  * Capozziello & De Laurentis (2011) Capozziello S., De Laurentis M., 2011, Physics Reports, 509, 167
  * Caruana et al. (2020) Caruana M., Farrugia G., Levi Said J., 2020, Eur. Phys. J. C, 80, 640
  * Chang et al. (2019) Chang Z., Zhao D., Zhou Y., 2019, Chinese Physics C, 43, 125102
  * Conley et al. (2011) Conley A., et al., 2011, ApJS, 192, 1
  * D’Agostino & Nunes (2022) D’Agostino R., Nunes R. C., 2022, Phys. Rev. D, 106, 124053
  * De Felice & Tsujikawa (2010) De Felice A., Tsujikawa S., 2010, Living Reviews in Relativity, 13, 1
  * De et al. (2022) De A., Saha D., Subramaniam G., Sanyal A. K., 2022
  * Deng & Wei (2018) Deng H.-K., Wei H., 2018, Eur. Phys. J. C, 78, 755
  * Dimakis et al. (2021) Dimakis N., Paliathanasis A., Christodoulakis T., 2021, Classical and Quantum Gravity, 38, 225003
  * Dimakis et al. (2022) Dimakis N., Paliathanasis A., Roumeliotis M., Christodoulakis T., 2022, arXiv preprint arXiv:2205.04680
  * Drozda et al. (2022) Drozda P., Hellwing W. A., Bilicki M., 2022, Phys. Rev. D, 106, 043513
  * Eisenstein & Hu (1998) Eisenstein D. J., Hu W., 1998, ApJ, 496, 605
  * Eisenstein et al. (2005) Eisenstein D. J., et al., 2005, The Astrophysical Journal, 633, 560
  * Esposito et al. (2022) Esposito F., Carloni S., Cianci R., Vignolo S., 2022, Physical Review D, 105, 084061
  * Ferreira et al. (2022) Ferreira J., Barreiro T., Mimoso J., Nunes N. J., 2022, Phys. Rev. D, 105, 123531
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Giocoli et al. (2018) Giocoli C., Baldi M., Moscardini L., 2018, Mon. Not. Roy. Astron. Soc., 481, 2813
  * Giostri et al. (2012) Giostri R., dos Santos M. V., Waga I., Reis R., Calvão M., Lago B. L., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 027
  * Gupta et al. (2022) Gupta S., Hellwing W. A., Bilicki M., García-Farieta J. E., 2022, Phys. Rev. D, 105, 043538
  * Harko et al. (2018) Harko T., Koivisto T. S., Lobo F. S., Olmo G. J., Rubiera-Garcia D., 2018, Physical Review D, 98, 084043
  * Harko et al. (2022) Harko T., Asadi K., Moshafi H., Sheikhahmadi H., 2022, ] 10.48550/ARXIV.2203.08907
  * Hassani & Lombriser (2020) Hassani F., Lombriser L., 2020, MNRAS, 497, 1885
  * Hayashi & Shirafuji (1979) Hayashi K., Shirafuji T., 1979, Physical Review D, 19, 3524
  * Hinshaw et al. (2013) Hinshaw G., et al., 2013, The Astrophysical Journal Supplement Series, 208, 19
  * Hohmann (2021) Hohmann M., 2021, Physical Review D, 104, 124077
  * Howlett et al. (2012) Howlett C., Lewis A., Hall A., Challinor A., 2012, J. Cosmology Astropart. Phys., 1204, 027
  * Hu et al. (2022) Hu K., Katsuragawa T., Qiu T., 2022, arXiv preprint arXiv:2204.12826
  * Huang et al. (2022) Huang Y., Zhang J., Ren X., Saridakis E. N., Cai Y.-F., 2022, Phys. Rev. D, 106, 064047
  * Jarosik et al. (2011) Jarosik N., et al., 2011, The Astrophysical Journal Supplement Series, 192, 14
  * Jiménez et al. (2018a) Jiménez J. B., Heisenberg L., Koivisto T., 2018a, Physical Review D, 98, 044048
  * Jiménez et al. (2018b) Jiménez J. B., Heisenberg L., Koivisto T. S., 2018b, Journal of Cosmology and Astroparticle Physics, 2018, 039
  * Jiménez et al. (2020a) Jiménez J. B., Heisenberg L., Koivisto T., Pekar S., 2020a, Phys. Rev. D, 101, 103507
  * Jiménez et al. (2020b) Jiménez J. B., Heisenberg L., Koivisto T., Pekar S., 2020b, Physical Review D, 101, 103507
  * Kessler & Scolnic (2017) Kessler R., Scolnic D., 2017, The Astrophysical Journal, 836, 56
  * Khyllep et al. (2023) Khyllep W., Dutta J., Saridakis E. N., Yesmakhanova K., 2023, Physical Review D, 107, 044022
  * Koussour et al. (2022) Koussour M., Shekh S. H., Bennai M., 2022, JHEAp, 35, 43
  * Lazkoz et al. (2019) Lazkoz R., Lobo F. S. N., Ortiz-Baños M., Salzano V., 2019, Phys. Rev. D, 100, 104027
  * Lee et al. (2022) Lee J., Ryu S., Baldi M., 2022
  * Lewis & Bridle (2002) Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511
  * Lewis & Challinor (2011) Lewis A., Challinor A., 2011, CAMB: Code for Anisotropies in the Microwave Background, Astrophysics Source Code Library, record ascl:1102.026 (ascl:1102.026)
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Liddle (2007) Liddle A. R., 2007, Monthly Notices of the Royal Astronomical Society: Letters, 377, L74
  * Lin & Zhai (2021) Lin R.-H., Zhai X.-H., 2021, Phys. Rev. D, 103, 124001
  * Linder (2009) Linder E. V., 2009, Physical Review D, 80, 123528
  * Linder (2010) Linder E. V., 2010, Phys. Rev. D, 81, 127301
  * Liu et al. (2022) Liu Y., Liao S., Liu X., Zhang J., An R., Fan Z., 2022, Mon. Not. Roy. Astron. Soc., 511, 3076
  * Mandal et al. (2020a) Mandal S., Sahoo P. K., Santos J. R. L., 2020a, Phys. Rev. D, 102, 024057
  * Mandal et al. (2020b) Mandal S., Wang D., Sahoo P., 2020b, Physical Review D, 102, 124029
  * Moresco (2015) Moresco M., 2015, Monthly Notices of the Royal Astronomical Society: Letters, 450, L16
  * Mostaghel et al. (2017) Mostaghel B., Moshafi H., Movahed S. M. S., 2017, Eur. Phys. J. C, 77, 541
  * Motohashi & Suyama (2015) Motohashi H., Suyama T., 2015, Physical Review D, 91, 085009
  * Nester & Yo (1999) Nester J. M., Yo H.-J., 1999, Chin. J. Phys., 37, 113
  * Padmanabhan (2003) Padmanabhan T., 2003, Physics Reports, 380, 235
  * Pan et al. (2018) Pan S., Mukherjee A., Banerjee N., 2018, Monthly Notices of the Royal Astronomical Society, 477, 1189
  * Pasqua et al. (2015) Pasqua A., Chattopadhyay S., Momeni D., Raza M., Myrzakulov R., 2015, Journal of Cosmology and Astroparticle Physics, 2017
  * Percival et al. (2010) Percival W. J., et al., 2010, Monthly Notices of the Royal Astronomical Society, 401, 2148
  * Perlmutter et al. (1999) Perlmutter S., et al., 1999, The Astrophysical Journal, 517, 565
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Planck Collaboration et al. (2020a) Planck Collaboration et al., 2020a, A&A, 641, A6
  * Planck Collaboration et al. (2020b) Planck Collaboration et al., 2020b, A&A, 641, A10
  * Pontzen et al. (2013) Pontzen A., Roškar R., Stinson G., Woods R., 2013, pynbody: N-Body/SPH analysis for python, Astrophysics Source Code Library, record ascl:1305.002 (ascl:1305.002)
  * Puchwein et al. (2013) Puchwein E., Baldi M., Springel V., 2013, Monthly Notices of the Royal Astronomical Society, 436, 348
  * Reid et al. (2019) Reid M. J., Pesce D. W., Riess A. G., 2019, Astrophys. J. Lett., 886, L27
  * Riess et al. (1998) Riess A. G., et al., 1998, The Astronomical Journal, 116, 1009
  * Ruan et al. (2022) Ruan C.-Z., Hernández-Aguayo C., Li B., Arnold C., Baugh C. M., Klypin A., Prada F., 2022, JCAP, 05, 018
  * SAHNI & STAROBINSKY (2000) SAHNI V., STAROBINSKY A., 2000, International Journal of Modern Physics D, 09, 373
  * Sahni et al. (2003) Sahni V., Saini T. D., Starobinsky A. A., Alam U., 2003, Journal of Experimental and Theoretical Physics Letters, 77, 201
  * Sahni et al. (2008) Sahni V., Shafieloo A., Starobinsky A. A., 2008, Phys. Rev. D, 78, 103502
  * Scolnic et al. (2018a) Scolnic D. M., et al., 2018a, The Astrophysical Journal, 859, 101
  * Scolnic et al. (2018b) Scolnic D. M., et al., 2018b, Astrophys. J., 859, 101
  * Solanki et al. (2023) Solanki R., Arora S., Sahoo P. K., Moraes P. H. R. S., 2023, Universe, 9, 12
  * Su et al. (2022) Su C., et al., 2022
  * Suzuki et al. (2012) Suzuki N., et al., 2012, The Astrophysical Journal, 746, 85
  * Tripp (1998) Tripp R., 1998, Astronomy and Astrophysics, 331, 815
  * Villaescusa-Navarro (2018) Villaescusa-Navarro F., 2018, Pylians: Python libraries for the analysis of numerical simulations, Astrophysics Source Code Library, record ascl:1811.008 (ascl:1811.008)
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro F., et al., 2020, The Astrophysical Journal Supplement Series, 250, 2
  * Wang et al. (2022) Wang W., Chen H., Katsuragawa T., 2022, Physical Review D, 105, 024060
  * Wilson & Bean (2022) Wilson C., Bean R., 2022
  * Xu et al. (2018) Xu T., Cao S., Qi J., Biesiada M., Zheng X., Zhu Z.-H., 2018, JCAP, 06, 042
  * Yu et al. (2018) Yu H., Ratra B., Wang F.-Y., 2018, The Astrophysical Journal, 856, 3
  * Zhang et al. (2018a) Zhang J., An R., Liao S., Luo W., Li Z., Wang B., 2018a, Phys. Rev. D, 98, 103530
  * Zhang et al. (2018b) Zhang J., An R., Liao S., Luo W., Li Z., Wang B., 2018b, Phys. Rev. D, 98, 103530
  * Zhang et al. (2019) Zhang J., An R., Luo W., Li Z., Liao S., Wang B., 2019, ApJ, 875, L11
  * Zhao (2022) Zhao D., 2022, The European Physical Journal C, 82, 1
  * Zhao et al. (2022) Zhao Y., Liu Y., Liao S., Zhang J., Liu X., Du W., 2022