High-redshift test of gravity using enhanced growth of small structures probed by the neutral hydrogen distribution

TL;DR: In this paper, the effects of modified gravity on neutral hydrogen abundance and power spectra were investigated using SHYBONE simulations, which model galaxy formation in fðRÞ modified gravity using the IllustrisTNG==================¯¯¯¯¯¯model, and they showed that this suppression can be detectable by SKA1-MID with approximately 1000 hours of exposure time.

Abstract: Future 21-cm intensity mapping surveys such as SKA can provide precise information on the spatial distribution of the neutral hydrogen (HI) in the postreionization epoch. This information will allow us to test the standard Λ cold dark matter paradigm and with that the nature of gravity. In this work, we employ the SHYBONE simulations, which model galaxy formation in fðRÞ modified gravity using the IllustrisTNG model, to study the effects of modified gravity on HI abundance and power spectra. We find that the enhanced growth low-mass dark matter halos experience in fðRÞ gravity at high redshifts alters the HI power spectrum and can be observable through 21-cm intensity mapping. Our results suggest that the HI power spectrum is suppressed by ∼13% on scales k ≲ 2h Mpc−1 at z ¼ 2 for F6, a fðRÞ model which passes most observational constraints. We show that this suppression can be detectable by SKA1-MID with 1000 hours of exposure time, making HI clustering a novel test of gravity at high redshift.

Summary

Introduction

• Durham Research Online Deposited in DRO: 24 September 2019 Version of attached le: Published Version Peer-review status of attached le: Peer-reviewed Citation for published item: Leo, Matteo and Arnold, Christian and Li, Baojiu (2019) 'A high-redshift test of gravity using enhanced growth of small structures probed by the neutral hydrogen distribution.', Physical review D., 100 (6). 064044.
• The model employs the so-called chameleon screening mechanism [16,17] to ensure that the modifications to standard gravity are suppressed and GR-like behavior is recovered in high-density regions like the Solar System.
• 21-cm intensity mapping can be used to trace the underlying distribution of matter [30–33] and with that the low-mass halos in the Universe (as suggested in [34]).

II. THEORETICAL MODELS AND SIMULATIONS

• With an appropriate choice of the functional form and parameters of fðRÞ, the theory can mimic the late time expansion history of a ΛCDM universe without explicitly having a cosmological constant Λ (the accelerated expansion in these theories is achieved via some form of quintessence/cosmological constant and is not due to the modification of gravity itself [40–42]).
• The theory is therefore fully specified by Ωm and the present-day value of the background scalar field, f̄R0.
• The chameleon screening has been described in great detail in the literature and thus the authors will not discuss it further here, but instead simply mention that it becomes effective when fR becomes close to zero, such that δR ≈ 8πGδρ according to Eq. (3), and then Eq. (2) reduces to the standard Poisson equation in Newtonian gravity.
• The screening is more likely to take place at earlier times when matter density is high and the background value of the scalar field, jf̄Rj, is small.
• As the authors will see in the next sections, 21-cm intensity mapping is sensitive to the abundance of halos down to 109 M⊙, making it a very promising probe of differences at the low-mass end of the halo mass function, without the need to resolve individual halos.

B. Full-physics simulations in MG

• In order to quantify how modifications to gravity affect the 21-cm signal, the authors analyze the SHYBONE simulations [46], a set of high-resolution full-physics hydrodynamical simulations of HS fðRÞ gravity, carried out with the moving mesh simulation code AREPO [47].
• The full-physics simulations use the IllustrisTNG hydrodynamical model [49–57], incorporating a prescription of star and black hole formation and feedback, gas cooling, galactic winds, and magnetohydrodynamics on a moving Voronoi mesh [53,57].
• The equations for fðRÞ gravity are solved to full nonlinearity in the Newtonian limit by the modified gravity solver in the code [46], fully capturing the effects of the chameleon screening.
• The postprocessing gives the total fraction of hydrogen that is nonionized: atomic (HI) and molecular hydrogen (H2).

A. Overall neutral hydrogen density

• The authors follow the common definition for the overall HI abundance, ΩHIðzÞ ¼ ρ̄HIðzÞ=ρc0, where ρ̄HIðzÞ is the mean HI density in their 064044-3 simulations at a given redshift z and ρc0 is the present-day critical density as defined above.
• First, the authors note that in Fig. 1 the HI abundance (for each model) predicted by S25 is higher than that measured from the low-resolution S62.
• A similar effect was found in [34] comparing the low- and high-resolution TNG simulations.
• This behavior can be understood as follows.
• At high redshifts (z≳ 4–5), modified gravity effects on the matter and halo distribution are screened for the models considered here, and thus F5 and F6 both behave similarly to GR.

B. HI mass in halos

• The authors present and discuss the halo HI mass function, i.e., the average HI mass enclosed in halos as a function of the halo mass.
• In Fig. 2, the authors also compare their best-fit curves with the GR results taken from [34] (cyan dotted lines).
• As the authors can see, the fitting results in the two works agree very well at z ¼ 2 for the entire halo mass range.
• It is nevertheless possible to correct the total HI abundance in their simulations for the missing contribution from high-mass halos, ΩcorrHI .

C. HI clustering

• In the previous subsections, the authors have focused on the overall HI abundance and the halo HI mass function.
• Similar trends can be found for the realspace power spectrum.
• In the case of S62, the authors show the HMFs for both the full-physics (dashed lines) and DMO simulations, while for S25 they only show the fullphysics ones (solid lines).

D. Additional tests

• To further check the above result about the different behavior of PHI in MG and GR, the authors have carried out two additional tests.
• To further show the importance of simulation resolution for accurate predictions of HI, Fig. 5 displays the number density of halos with HI mass MHI ≥ 106 M⊙ (i.e., the halos that contribute significantly to PHI).
• For masses lower 064044-7 than Mpeak200 , the number of HI-rich halos decreases, though the halo mass function keeps increasing.
• This can explain why the authors find a very similar degree of suppression in the power spectra of these two models in the previous subsection.
• Because of its less efficient chameleon screening, one may naively expect that F5 is able to turn more initial density peaks into halos (above a certain mass) than F6, and so this model always has more HI-rich halos than F6 at a given time.

E. Observational forecast

• The authors compute the expected 1σ errors for 1000 observing hours, as shown by the shaded areas in the central lower subpanels of Fig.
• In the calculation, the authors have used the monopole of the HI redshift space power spectrum for GR measured from the S25 simulation.
• The integrated signal-tonoise (S/N) ratios for distinguishing MG to GR are shown in Table II for the two redshifts considered in this analysis.

F. Uncertainties in subgrid physics

• Current state-of-the-art simulations are still unable to fully resolve the formation and evolution of stars and galaxies from first principles, largely due to the huge dynamical range between the scales at which star formation and black hole accretion take place and the scales which must be covered in order to realistically reproduce the large-scale structure of the Universe.
• An important point to understand is to what extent the uncertainties in active galactic nuclei (AGN) and stellar feedback can affect the results above.
• Much stronger variations to the feedback mechanism may alter the HI content within halos of a given mass range and consequently the HI power spectrum; such strong changes would nevertheless lead to tensions with the low redshift observables used to tune the IllustrisTNG model [53].
• From the discussion above, the authors expect the main conclusions of this work to be relatively robust against changes of a subgrid physics parameter in the simulations.

IV. SUMMARY AND CONCLUSIONS

• The screening mechanism of chameleon-type MG models, such as fðRÞ gravity, is particularly efficient at high redshift and for massive objects.
• The authors propose that this enhancement should be observable through 21-cm intensity mapping and use the SHYBONE simulations, a set of state-of-the-art hydrodynamical simulations employing the IllustrisTNG model which were carried out for two different fðRÞ gravity models (F5 and F6), to analyze the viability of this approach.
• The authors fitting curves for GR are in agreement with those found in [34].

Figures

FIG. 2. HI mass contained in each halo at z ¼ 3 and z ¼ 2 measured from the S25 simulations for GR, F6, and F5, respectively (as labeled). The green solid line in each panel shows the best-fit curve obtained from Eq. (9). The best-fit values of the free parameters are displayed in Table I. The cyan dotted line indicates the best-fit curve of Eq. (9) found in [34] for GR. The blue vertical line indicates the mass of halos with around 50 DM (simulation) particles in S25.

TABLE II. The integrated S/N ratios for distinguishing a MG model from GR using redshift space PHIðkÞ with kmax ¼ 2h Mpc−1, with ðS=NÞ2 ≡Pkmaxkmin ½PMGHI ðkÞ − PGRHI ðkÞ 2=σ2ðkÞ; here σðkÞ is the expected 1σ error from SKA1-MID for 1000 hours, while kmin is set by the value of the box length L ¼ 25h−1 Mpc of our high-resolution simulations, kmin ¼ 2π=L. The above results are calculated using the redshift space HI power spectra measured from S25 only.

FIG. 3. Neutral hydrogen power spectra and halo mass functions at z ¼ 2 (upper panels) and z ¼ 3 (lower panels) for GR (black), F6 (blue), and F5 (green). Left panels: Real-space HI power spectra PHI. Dashed lines show the results for S62, while solid lines show those for S25. The real-space power spectra are taken from Ref. [46]. Central panels: Same as the left panels but for the redshift space (monopole) HI power spectra. The cyan shaded areas in the lower subpanels represent the expected errors for SKA1-MID measurements for GR, assuming 1000 observing hours. The error bars are calculated using the GR redshift space power spectrum measured from S25. Right panels: Differential halo mass functions as shown in [69], for comparison. Solid (dashed) lines show the results from full-physics S25 (S62) simulations. Symbols indicate the results from DMO S62 simulations. The lower subpanels show the relative differences to GR.

FIG. 4. The ratios of the F5 (green) and F6 (blue) HI power spectra with respect to GR for the actual real-space HI power spectrum (symbols) and the halo HI power spectrum (lines), PHI;halo, as defined in the text, at z ¼ 2.0 (upper panel) and z ¼ 3.0 (lower panel). Solid lines and full circles show the results from S25, while dashed lines and empty circles show the results from S62.

FIG. 5. Number density of halos with HI mass ≥ M⋆HI ¼ 106 M⊙, for GR (black), F6 (blue), and F5 (green) at z ¼ 2.0 (left panel) and 3.0 (right panel). Solid lines show the results for S25, while dashed lines represent those from S62. The vertical grey dashed lines indicate halos with ∼50 particles in S62. We do not show the corresponding limit for S25 because it is at around ∼6 × 108 M⊙, so it is outside the halo mass range shown in the figure.

FIG. 1. Top panel: Overall HI abundance,ΩHIðzÞ ¼ ρ̄HIðzÞ=ρc0, where ρ̄HIðzÞ is the mean HI density, from GR (black), F6 (blue), and F5 (green), compared with observationally measured values (symbols). Solid lines refer to S25 simulations, while dashed lines refer to S62 simulations. Bottom panel: The relative differences of the simulation predictions from F6 (blue) and F5 (green) with respect to GR.

Full text

Durham Research Online
Deposited in DRO:
24 September 2019
Version of attached le:
Published Version
Peer-review status of attached le:
Peer-reviewed
Citation for published item:
Leo, Matteo and Arnold, Christian and Li, Baojiu (2019) 'A high-redshift test of gravity using enhanced
growth of small structures probed by the neutral hydrogen distribution.', Physical review D., 100 (6). 064044.
Further information on publisher's website:
https://link.aps.org/doi/10.1103/PhysRevD.100.064044
Publisher's copyright statement:
Reprinted with permission from the American Physical Society: Leo, Matteo, Arnold, Christian Li, Baojiu (2019). A
high-redshift test of gravity using enhanced growth of small structures probed by the neutral hydrogen distribution.
Physical Review D 100: 064044.
c

2019 by the American Physical Society. Readers may view, browse, and/or download
material for temporary copying purposes only, provided these uses are for noncommercial personal purposes. Except as
provided by law, this material may not be further reproduced, distributed, transmitted, modied, adapted, performed,
displayed, published, or sold in whole or part, without prior written permission from the American Physical Society.
Additional information:
Use policy
The full-text may be used and/or reproduced, and given to third parties in any format or medium, without prior permission or charge, for
personal research or study, educational, or not-for-prot purposes provided that:
•
a full bibliographic reference is made to the original source
•
a link is made to the metadata record in DRO
•
the full-text is not changed in any way
The full-text must not be sold in any format or medium without the formal permission of the copyright holders.
Please consult the full DRO policy for further details.
Durham University Library, Stockton Road, Durham DH1 3LY, United Kingdom
Tel : +44 (0)191 334 3042 | Fax : +44 (0)191 334 2971
https://dro.dur.ac.uk

High-redshift test of gravity using enhanced growth of small
structures probed by the neutral hydrogen distribution
Matteo Leo ,
1,2,*
Christian Arnold,
2
and Baojiu Li
2
1
Institute for Particle Physics Phenomenology, Department of Physics, Durham University,
Durham DH1 3LE, United Kingdom
2
Institute for Computational Cosmology, Department of Physics, Durham University,
Durham DH1 3LE, United Kingdom
(Received 16 January 2019; published 23 September 2019)
Future 21-cm intensity mapping surveys such as SKA can provide precise information on the spatial
distribution of the neutral hydrogen (HI) in the postreionization epoch. This information will allow us to
test the standard Λ cold dark matter paradigm and with that the nature of gravity. In this work, we employ
the
SHYBONE
simulations, which model galaxy formation in fðRÞ modified gravity using the IllustrisTNG
model, to study the effects of modified gravity on HI abundance and power spectra. We find that the
enhanced growth low-mass dark matter halos experience in fðRÞ gravity at high redshifts alters the HI
power spectrum and can be observable through 21-cm intensity mapping. Our results suggest that the HI
power spectrum is suppressed by ∼ 13% on scales k ≲ 2h Mpc
−1
at z ¼ 2 for F6, a fðRÞ model which
passes most observational constraints. We show that this suppression can be detectable by SKA1-MID with
1000 hours of exposure time, making HI clustering a novel test of gravity at high redshift.
DOI: 10.1103/PhysRevD.100.064044
I. INTRODUCTION
Our standard model of cosmology—the Λ cold dark
matter model (ΛCDM)—has proved very successful in
describing almost all currently available observational data
of the Universe. Its underlying theory of gravity, Einstein’s
general relativity (GR), has been tested to remarkably high
precision on small scales [1]. In recent years, in the wake
of high-precision astronomical observations, tests of GR on
cosmological scales have become possible and common-
place as well [2], although until now these tests have
primarily focused on comparatively high-mass objects and
low redshifts (e.g., [3–12]). Because of the screening
mechanisms which many alternatives to GR employ to
pass the stringent Solar System tests, these objects are less
suited to distinguish GR from alternative, modified gravity
(MG) theories with screening mechanisms.
As a representative example, we consider a particular
one of these MG models in this paper—fðRÞ gravity [13],
though we expect our conclusions to hold at least quali-
tatively for general thin-shell screening [14] models. fðRÞ
gravity is a generalization of GR which alters cosmic
structure formation through a factor-of-4=3 enhanced
gravitational force. We adopt the popular variant proposed
in Ref. [15], which, with certain choices of the free
parameters of the model, can produce a cosmic expan-
sion history very close to that of the ΛCDM paradigm.
The model employs the so-called chameleon screening
mechanism [16,17] to ensure that the modifications to
standard gravity are suppressed and GR-like behavior is
recovered in high-density regions like the Solar System.
The model considered here has been widely studied using
numerical simulations, e.g., [9,18–24]. In this work we
consider two instances of it [15]: F6 and F5, with its model
parameter f
R0
equal to −10
−5
and −10
−6
, respectively (see
Sec. II A for more details).
Although the constraints from our local environment are
very tight, as mentioned above, the previous constraints
from cosmological scales are much weaker since the
objects used in these tests are generally more massive
and well screened. One way to overcome this limitation is
to study low-mass objects which are less likely to be
screened and hence experience larger deviations from GR.
However, a major challenge to this approach is the
difficulty to accurately detect, resolve, and trace such small
objects in observations, even at low redshifts.
In this paper, we propose a novel test of gravity at
intermediate scales and high redshifts (z ≥ 2), using the
distribution of neutral hydrogen (HI) in our Universe,
which is observable in 21-cm experiments (some current
and future instruments of this kind include SKA [25],
MeerKAT [26], LOFAR [27], CHIME [28], and BINGO
[29]). 21-cm intensity mapping can be used to trace the
underlying distribution of matter [30–33] and with that the
low-mass halos in the Universe (as suggested in [34]).
In order to determine how possible deviations from GR
*
matteo.leo@durham.ac.uk
PHYSICAL REVIEW D 100, 064044 (2019)
2470-0010=2019=100(6)=064044(12) 064044-1 © 2019 American Physical Society

would affect the HI distribution, we employ the
SHYBONE
simulations, a set of full-physics hydrodynamical simula-
tions of fðRÞ modified gravity. Comparing the fðRÞ
simulations to their ΛCDM counterpart allows us to
quantify the size of the MG effects on HI observables
(such as the overall neutral hydrogen abundance and the HI
power spectrum), and to assess if these effects can be
observable with future 21-cm intensity mapping experi-
ments. HI clustering has been proposed as a probe for a
number of nonstandard cosmological models, e.g., massive
neutrinos, warm DM, dark energy, and modified gravity
[35–39], but this study reveals new features, thanks to the
high resolution of our simulations.
This paper is structured as follows. In Sec. II we briefly
introduce the two instances of Hu-Sawicki (HS) fðRÞ
gravity [15] used in our investigation and the suite of
hydrodynamical simulations employed to quantify the
abundance and clustering of HI. In Sec. III, we show
and discuss the main results of this paper, including the
overall neutral hydrogen density (III A), the HI abundance
in halos (III B), and the HI power spectra in both real and
redshift space ( III C). Additional tests to explain the physics
behind our results are performed in Sec. III D and obser-
vational forecasts for a future 21-cm intensity mapping
experiment are discussed in Sec. III E. In Sec. III F,we
comment on the dependence of our results on the galaxy
formation model employed in our simulations. Finally, we
conclude our findings in Sec. IV.
II. THEORETICAL MODELS
AND SIMULATIONS
A. f(R) gravity
fðRÞ gravity is a popular class of MG models that is
obtained by adding a scalar function fðRÞ to the Ricci
scalar R in the standard Einstein-Hilbert action [13] of
general relativity. With an appropriate choice of the func-
tional form and parameters of fðRÞ, the theory can mimic
the late time expansion history of a ΛCDM universe
without explicitly having a cosmological constant Λ (the
accelerated expansion in these theories is achieved via
some form of quintessence/cosmological constant and is
not due to the modification of gravity itself [40–42]).
The action for the f
ðRÞ gravity can be written as
S ¼
Z
d
4
x
ffiffiffiffiffiffi
−g
p

R þ fðRÞ
16πG
þ L
m

; ð1Þ
where G is the gravitational constant, g is the determinant
of the metric, g
μν
, and L
m
is the standard matter/radiation
Lagrangian density. The simulations considered here
employ the weak-field and quasistatic limit (see [43] for
more details on the validity of these approximations), so
that the equations of motion obtained by varying the action
in Eq. (1) can be simplified to a (modified) Poisson
equation plus an equation for the scalar degree of freedom,
f
R
≡ dfðRÞ=dR (the so-called scalar field),
∇
2
Φ ¼
16πG
3
δρ −
1
6
δR; ð2Þ
∇
2
f
R
¼
1
3
ðδR − 8πGδρÞ; ð3Þ
where δρ ≡ ρ −
¯
ρ and δR ≡ R −
¯
R are the matter density
perturbation and the Ricci scalar perturbation, respectively
(and
¯
ρ and
¯
R are their background values).
The HS variant of the theory [15] uses
fðRÞ¼−m
2
c
1
ð
R
m
2
Þ
n
c
2
ð
R
m
2
Þ
n
þ 1
; ð4Þ
where m is a new mass scale of the model, m
2
≡ Ω
m
H
2
0
,
H
0
is the Hubble constant, Ω
m
is the total nonrelativistic
matter energy density at present time in units of the pre-
sent-day critical energy density of the Universe, ρ
c0
≡
3H
2
0
=8πG, and c
1
, c
2
, n are model parameters. We choose
n ¼ 1 hereafter for simplicity. Furthermore, if we tune the
parameters c
1
and c
2
such that
c
1
c
2
¼ 6
Ω
Λ
Ω
m
and
c
2
jRj
m
2
≫ 1; ð5Þ
the model leads to a cosmic expansion history which is very
close to that of a ΛCDM universe [15]. Ω
Λ
in the above
equation represents the cosmological constant energy
density in units of ρ
c0
for the ΛCDM universe; in the case
of fðRÞ gravity, Ω
Λ
still enters in the theory as a parameter.
If one further assumes that the Universe is spatially flat,
then Ω
Λ
is simply given by Ω
Λ
¼ 1 − Ω
m
. This is our
default assumption in this work.
The scalar field f
R
in this model can be approximated as
f
R
≡
dfðRÞ
dR
≈−
c
1
c
2
2

m
2
R

2
; ð6Þ
and its background evolution can be expressed in terms of
the background Ricci scalar
¯
R,
¯
f
R
ðaÞ¼
¯
f
R0

¯
R
0
¯
RðaÞ

2
; ð7Þ
where
¯
R
0
is the value of the Ricci scalar today and
¯
RðaÞ¼3m
2

a
−3
þ 4
Ω
Λ
Ω
m

: ð8Þ
The theory is therefore fully specified by Ω
m
and the
present-day value of the background scalar field,
¯
f
R0
.
In order to satisfy the stringent constraints on the
possible deviations from standard gravity in our local
MATTEO LEO, CHRISTIAN ARNOLD, and BAOJIU LI PHYS. REV. D 100, 064044 (2019)
064044-2

environment [1], the theory employs the above-mentioned
chameleon screening mechanism [16,17] to suppress
modifications to gravity and restore GR in high-density
regions. The chameleon screening has been described in
great detail in the literature and thus we will not discuss
it further here, but instead simply mention that it becomes
effective when f
R
becomes close to zero, such that
δR ≈ 8πGδρ according to Eq. (3), and then Eq. (2) reduces
to the standard Poisson equation in Newtonian gravity.
The screening is more likely to take place at earlier times
when matter density is high and the background value of
the scalar field, j
¯
f
R
j, is small. At a given time, this
mechanism screens regions where their density is high
and therefore the Newtonian potential is deep. The tran-
sition between screened and unscreened regimes depends
on the choice of
¯
f
R0
.InfðRÞ gravity, the speed of the
gravitational wave is equal to the speed of light and the
model passes recent constraints from gravitational wave
observations [44], making it one of the most promising
alternatives to GR.
As mentioned in the Introduction, in this work we will
focus on the F5 and F6 instances of HS fðRÞ gravity, for
which
¯
f
R0
is equal to −10
−5
and −10
−6
, respectively. In F5,
the effects of MG are stronger than in F6 and this model is
now ruled out by observational constraints [12] (see [45]
for a review about the recent constraints on chameleon
gravity). However, it is used here as a toy model to assess
the effects of a stronger deviation from GR on the HI
distribution. The F6 model is consistent with most
cosmological observations.
An important characteristic of the chameleon screening
is that this mechanism becomes inefficient for small
structures at high redshift, while more massive objects
and denser environments become unscreened at later
times (lower redshift). At high redshift, low-mass halos
are already unscreened and affected by modified gravity.
As a consequence, by observing such small structures we
can, in principle, place constrains on the fðRÞ deviations
from GR. As we will see in the next sections, 21-cm
intensity mapping is sensitive to the abundance of halos
down to 10
9
M
⊙
, making it a very promising probe of
differences at the low-mass end of the halo mass function,
without the need to resolve individual halos.
B. Full-physics simulations in MG
In order to quantify how modifications to gravity affect
the 21-cm signal, we analyze the
SHYBONE
simulations
[46], a set of high-resolution full-physics hydrodynamical
simulations of HS fðRÞ gravity, carried out with the
moving mesh simulation code
AREPO
[47]. The suite
includes two subsets of simulations: a large-box set with
a box size of L ¼ 62h
−1
Mpc (S62 hereafter) and a small-
box set for which L ¼ 25h
−1
Mpc (S25 hereafter), both
with roughly 2 × 512
3
resolution elements [h is the
dimensionless Hubble constant, given by h ≡ H
0
=
ð100 km s
−1
Mpc
−1
Þ]. The S62 simulations have a mass
resolution of m
DM
¼ 1.3 × 10
8
h
−1
M
⊙
for DM particles
and roughly m
gas
¼ 2.5 × 10
7
h
−1
M
⊙
for gas cells, and
they have been run for GR, F6, and F5 up to z ¼ 0.
The S25 simulations have a mass resolution of m
DM
¼
8.4 × 10
6
h
−1
M
⊙
and m
gas
¼ 1.6 × 10
6
h
−1
M
⊙
and have
been run for the same three models, up to z ¼ 0 for GR and
F6 and up to z ¼ 1 for F5 (the enhanced gravitational
interactions in the F5 model considerably increase the
computational cost of the simulations compared to their
GR counterpart). S25 features a higher resolution, but its
smaller box means that we inevitably lose some informa-
tion of large-scale modes and massive halos (see the
discussion in the next section). The S62 suite features also
DM-only (DMO hereafter) counterparts for all the runs,
which are used to compare the halo mass function from
full-physics and DMO simulations below. All simulations
adopt the Planck 2016 [48] cosmology with Ω
m
¼ 0.3089,
Ω
B
¼ 0.0486, Ω
Λ
¼ 0.6911, h ¼ 0.6774, σ
8
¼ 0.8159,
and n
s
¼ 0.9667, where Ω
B
is the present-day baryon
density parameter, σ
8
is the root-mean-squared matter
density fluctuation over spherical regions with radius
8h
−1
Mpc at z ¼ 0, and n
s
is the index of the primordial
power spectrum.
The full-physics simulations use the IllustrisTNG hydro-
dynamical model [49–57], incorporating a prescription of
star and black hole formation and feedback, gas cooling,
galactic winds, and magnetohydrodynamics on a moving
Voronoi mesh [53,57]. The equations for fðRÞ gravity are
solved to full nonlinearity in the Newtonian limit by the
modified gravity solver in the code [46], fully capturing the
effects of the chameleon screening.
To calculate the neutral hydrogen fraction in each
Voronoi cell, we follow the prescription in Sec. n 2.2 in
[34]. For non-star-forming gas, we use the neutral hydrogen
fraction calculated on the fly in the simulations, while for
star-forming gas we postprocess the outputs, recalculating
the neutral hydrogen fraction in each cell assuming a
temperature of T ¼ 10
4
K and following the approach in
[58] to take into account self-shielding corrections. The
postprocessing gives the total fraction of hydrogen that is
nonionized: atomic (HI) and molecular hydrogen (H
2
).
Because we are solely interested in HI, we calculate and
subtract the fraction of H
2
for each Voronoi cell as in [34].
III. RESULTS
A. Overall neutral hydrogen density
In Fig. 1 we show the overall neutral hydrogen density
measured from the S62 (in the range 0 ≤ z ≤ 3 for all
models) and S25 (in the range 0 ≤ z ≤ 5 for GR and
F6 and 1 ≤ z ≤ 5 for F5) simulations. We follow the
common definition for the overall HI abundance, Ω
HI
ðzÞ¼
¯
ρ
HI
ðzÞ=ρ
c0
, where
¯
ρ
HI
ðzÞ is the mean HI density in our
HIGH-REDSHIFT TEST OF GRAVITY USING ENHANCED … PHYS. REV. D 100, 064044 (2019)
064044-3

simulations at a given redshift z and ρ
c0
is the present-day
critical density as defined above. We also show a selection
of observational data for the HI abundance at different
redshifts from [59–63].
First, we note that in Fig. 1 the HI abundance (for each
model) predicted by S25 is higher than that measured from
the low-resolution S62. A similar effect was found in [34]
comparing the low- and high-resolution TNG simulations.
This discrepancy between simulations at different resolu-
tion can be understood as follows. The neutral hydrogen in
the postreionization epoch is concentrated in halos, where
shielding effects screen them from ionization. It was shown
(see, e.g., [34]; but see also the next section) that there is a
significant amount of HI in halos with masses as low as
10
9
M
⊙
at z ≤ 5. This implies that resolving halos of this
mass in simulations is essential to measure the HI abun-
dance accurately. However, because of its lower resolution,
S62 does not fully resolve halos with masses <10
10
M
⊙
,
and therefore predicts a lower value for the HI abundance at
all redshifts considered in this analysis. On the other hand,
S25 can resolve halos down to 6 × 10
8
M
⊙
, and therefore
produces more reliable results of the HI abundance.
However, we note that because of their small-box size,
we do not have a statistically robust sample of halos with
masses ≳10
12
M
⊙
in S25. This explains why our high-
resolution simulations predict slightly lower HI abundance
at z ≤ 5 compared to that measured in [34] for GR using
the TNG-100 simulation (performed in a box of comoving
length L ¼ 75h
−1
Mpc at the same resolution as S25).
Considering the differences between fðRÞ gravity and
GR, in S25 the ratios of Ω
HI
ðzÞ with respect to GR (solid
lines in the lower panel of Fig. 1) show a similar trend in F6
and F5. Indeed, in both models, Ω
HI
is similar to GR at high
redshifts, larger than GR at intermediate redshifts (with an
enhanced peak at z ¼ 1 and z ¼ 2 for F6 and F5,
respectively) and falls below the GR values for z<1.
Overall, HI is ∼5% (18%) more abundant in F6 (F5) than in
GR at z ¼ 3.Atz ¼ 2 there is ∼12% (22%) more HI in F6
(F5) than GR, while at z ¼ 1 we find more HI in F6 than in
F5 and GR. This behavior can be understood as follows. At
high redshifts (z ≳ 4–5), modified gravity effects on the
matter and halo distribution are screened for the models
considered here, and thus F5 and F6 both behave similarly
to GR. At intermediate redshifts (z ¼ 2 ∼ 3), low-mass
halos in F5 and F6 become unscreened and experience
enhanced growth, leading to increased abundance of these
low-mass objects compared to GR. Since neutral hydrogen
can survive only in self-shielding halos, this implies that in
F5 and F6 there are more HI-hosting halos than in GR and,
consequently, these models are characterized by larger
overall HI abundance. At low redshifts, baryonic effects
become important. At these redshifts, we suspect that
processes of gas heating can be more efficient in MG than
in GR, reducing the overall HI abundance in F5 and F6
compared to GR. A closer inspection of Fig. 1 suggests
that F6 behaves as a “retarded” version of F5, with the
maximum enhancement with respect to GR shifted to lower
redshifts. This is expected as the screening is more efficient
in F6 and thus leads to a later onset of the MG force
enhancement compared to F5.
B. HI mass in halos
In this subsection, we present and discuss the halo HI
mass function, i.e., the average HI mass enclosed in halos
as a function of the halo mass. In the postreionization
epoch, the majority of HI resides in halos, and therefore an
accurate knowledge of the halo HI mass function can be
used to predict the HI power spectrum in real and redshift
space without requiring the full hydrodynamical simulation
apparatus, but instead by painting the HI on top of dark
matter halos from DMO simulations [64]. This approxi-
mation has been applied to extract information on the HI
power spectrum for nonstandard cosmological models of
dark matter and dark energy [35–37].
Here, for each halo in our simulations we measure its
enclosed HI mass. The halos are identified using the
SUBFIND
[65] algorithm implemented in
AREPO
. The halo
mass M
halo
is defined as M
200
, the mass contained in a
sphere of radius r
200
, within which the average density is
200 times the critical density at the specified redshift.
In Fig. 2, we plot the HI mass as a function of the host halo
mass for halos in the S25 simulations. As the figure shows,
the HI mass increases monotonically with the halo mass on
average in all models considered in our analysis, which is in
agreement with what was found in [34]. The plot also
shows that HI is present within very low-mass halos,
FIG. 1. Top panel: Overall HI abundance, Ω
HI
ðzÞ¼
¯
ρ
HI
ðzÞ=ρ
c0
,
where
¯
ρ
HI
ðzÞ is the mean HI density, from GR (black), F6
(blue), and F5 (green), compared with observationally measured
values (symbols). Solid lines refer to S25 simulations, while
dashed lines refer to S62 simulations. Bottom panel: The relative
differences of the simulation predictions from F6 (blue) and F5
(green) with respect to GR.
MATTEO LEO, CHRISTIAN ARNOLD, and BAOJIU LI PHYS. REV. D 100, 064044 (2019)
064044-4

Frequently asked questions

Q1. What can be done to break the potential degeneracy between modified expansion history and structure growth?

As a side remark, the authors note that 21-cm intensity mapping can also give accurate information of the expansion history [76–78], which may be used to break potential degeneracy between modified expansion history and structure growth. 

Q2. What was the HI distribution used in the IllustrisTNG model?

The IllustrisTNGmodel was calibrated against the stellar mass function, stellar mass fraction, galaxy central black hole masses and gas fraction, galaxy sizes at redshift z ¼ 0, and the cosmic star formation rate history. 

Q3. What is the effect of the HI power spectrum in redshift space?

For the monopole of the redshift space power spectrum, the authors find that at z ¼ 3 the F6 (F5) PHI;redshiftðkÞ is suppressed by 8% (14%) for k ∼ 2h Mpc−1, while the effect is even stronger at higher wave numbers. 

Q4. What is the main characteristic of the chameleon screening?

An important characteristic of the chameleon screening is that this mechanism becomes inefficient for small structures at high redshift, while more massive objects and denser environments become unscreened at later times (lower redshift). 

Q5. How is the fR0 a good probe of differences at the low mass end?

As the authors will see in the next sections, 21-cm intensity mapping is sensitive to the abundance of halos down to 109 M⊙, making it a very promising probe of differences at the low-mass end of the halo mass function, without the need to resolve individual halos. 

Q6. What is the halo mass function in the S62 suite?

The S62 suite features also DM-only (DMO hereafter) counterparts for all the runs, which are used to compare the halo mass function from full-physics and DMO simulations below. 

Q7. Why do fR models have a working chameleon screening mechanism?

This is because viable fðRÞ models with a working chameleon screening mechanism to restore GR in the Solar System must have practically identical expansion history to ΛCDM, but structures in these models still grow differently at high z.