Realistic simulations of galaxy formation in f(R) modified gravity

TL;DR

This study uses advanced cosmological simulations to explore how f(R) modified gravity influences galaxy formation, dark matter distribution, and observable signals like the hydrogen power spectrum, providing new insights into gravity's role in cosmic structure.

Contribution

First simulation study combining full physics galaxy formation models with f(R) gravity, revealing effects on galaxy properties and large-scale structure observables.

Findings

f(R)-gravity affects the neutral hydrogen power spectrum at high redshift by 20%

More disc galaxies form in F6 compared to GR, fewer in F5

Back reaction between AGN feedback and modified gravity varies between F6 and F5 models

Abstract

We have carried out a set of cosmological hydrodynamical simulations that follow galaxy formation in $f(R)$ modified gravity models. Our simulations employ the Illustris-TNG full physics model and a new modified gravity solver in the AREPO code. For the first time we are able to investigate the degeneracy in the matter power spectrum between the effects of $f(R)$-gravity and feedback from active galactic nuclei (AGN), and the imprint of modified gravity on the properties of galaxies and on the distribution of dark matter, gas and stars in the universe. $f(R)$-gravity has an observable effect on the neutral hydrogen power spectrum at high redshift at a level of 20%. For both the F6 and F5 models, this is significantly larger than the predicted errors for the SKA1-MID survey, making this probe a powerful test of gravity on large scales. A similar effect is present in the power spectrum of the stars at high redshift. We also show that rotationally supported disc galaxies can form in $f(R)$-gravity, even in the partially screened regime. Our simulations indicate that there might be more disc galaxies in F6 compared to GR, and fewer in F5. Finally, we show that the back reaction between AGN feedback and modified gravity in the matter power spectrum is not important in the F6 model but has a sizeable effect in F5.

Equations

$$S=\int {\mathrm{d}}^{4}x\sqrt{-g}\left[\frac{R+f(R)}{16\pi G}+{\mathcal{L}}_m\right],$$

$$G_{\mu \nu }+f_{\mathrm{R}}{R}_{\mu \nu }-\left(\frac{f}{2}-\square f_{\mathrm{R}}\right)g_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu }{f}_{\mathrm{R}}=8\pi G{T}_{\mu \nu },$$

$${\nabla }^2\Phi =\frac{16\pi G}{3}\delta \rho -\frac{1}{6}\delta R,$$

$${\nabla }^2{f}_{\mathrm{R}}=\frac{1}{3}\left(\delta R-8\pi G\delta \rho \right).$$

$$f(R)=-m^2\frac{c_1{\left(\frac{R}{m^2}\right)}^n}{c_2{\left(\frac{R}{m^2}\right)}^n+1},$$

$$\frac{c_1}{c2}=6\frac{{\Omega }_{\Lambda }}{{\Omega }_m}$$

$$f_{\mathrm{R}}\equiv \frac{\mathrm{d}f(R)}{\mathrm{d}R}\approx -\frac{c_1}{c_2^2}{\left(\frac{m^2}{R}\right)}^2,$$

$${\bar{f}}_{\mathrm{R}}(a)={\bar{f}}_{R0}{\left[\frac{{\bar{R}}_0}{\bar{R}(a)}\right]}^2,$$

$$\bar{R}=3m^2\left[a^{-3}+4\frac{{\Omega }_{\Lambda }}{{\Omega }_m}\right].$$

$${\nabla }^2{\mathrm{e}}^u=-\frac{1}{3c^2{\bar{f}}_R(a)}\left[\bar{R}(a)\left(1-{\mathrm{e}}^{-\frac{u}{2}}\right)+8\pi G\delta \rho \right],$$

$$u\equiv \sqrt{\frac{f_{\mathrm{R}}}{{\bar{f}}_{\mathrm{R}}(a)}},$$

$${\nabla }^2\left(u^2\right)=\frac{1}{3}\left[\frac{\bar{R}(a)}{{\bar{f}}_{\mathrm{R}}(a)}\left(\frac{1}{u}-1\right)-\frac{8\pi G}{{\bar{f}}_{\mathrm{R}}(a)}\delta \rho \right].$$

$${\nabla }^2\left(u_{\mathrm{i},\mathrm{j},\mathrm{k}}^2\right)$$

$$\alpha \equiv \frac{\bar{R}(a)}{18{\bar{f}}_{\mathrm{R}}(a)}$$

$$p$$

$$q$$

$$u_{\mathrm{i},\mathrm{j},\mathrm{k}}^3+p{u}_{\mathrm{i},\mathrm{j},\mathrm{k}}+q=0.$$

$$D\equiv {\left(\frac{q}{2}\right)}^2+{\left(\frac{p}{3}\right)}^3.$$

$${\tilde{u}}_{\mathrm{i},\mathrm{j},\mathrm{k}}$$

$$C$$

$${\tilde{u}}_{\mathrm{i},\mathrm{j},\mathrm{k}}$$

$$\Theta$$

$$e_{\mathrm{i},\mathrm{j},\mathrm{k}}=\frac{r_{\mathrm{i},\mathrm{j},\mathrm{k}}}{\bar{\rho }(a)l^2\beta },$$

$$r_{\mathrm{i},\mathrm{j},\mathrm{k}}={\tilde{u}}_{\mathrm{i},\mathrm{j},\mathrm{k}}^3+p{\tilde{u}}_{\mathrm{i},\mathrm{j},\mathrm{k}}+q,$$

$$\max (e_{\mathrm{i},\mathrm{j},\mathrm{k}})<10^{-2}.$$

$${\mathcal{L}}_l(u_l)=0,$$

$${\mathcal{L}}_l({\tilde{u}}_l)=r_l,$$

$$\mathcal{R}{\tilde{u}}_l\equiv \sum _{i=0}^8\frac{{\tilde{u}}_l^i}{8},$$

$${\mathcal{L}}_L(u_L)={\mathcal{L}}_L(\mathcal{R}{\tilde{u}}_l)-\mathcal{R}{r}_l,$$

$${\tilde{u}}_l^{\mathrm{new}}={\tilde{u}}_l+\mathcal{P}({\tilde{u}}_L-\mathcal{R}{\tilde{u}}_l).$$

$${\nabla }^2\Phi =4\pi G(\delta \rho +\delta {\rho }_{\mathrm{eff}}),$$

$$\delta {\rho }_{\mathrm{eff}}=\frac{\delta \rho }{3}-\frac{\delta R}{24\pi G}.$$

Full text

cosmology: theory – methods: numerical

The physical nature of gravity and dark energy are central unsolved problems in modern physics. Einstein’s General Relativity (GR) is generally accepted as the standard theory of gravity. GR has been empirically verified, to remarkably high precision, on small scales, but on cosmological scales the available tests are still not strongly constraining. Planned astronomical surveys will deliver data capable of distinguishing it from alternative theories of gravity. However, predictions for these observables require accurate simulations of the formation of cosmic structure in GR and alternatives, which must include the interactions between baryons and dark matter. Here, we present a set of cosmological hydrodynamical simulations of an attractive alternative to GR, the Hu & Sawicki $f(R)$-gravity, which is representative of the large class of chameleon-type screening modified gravity models. These simulations are the first to include galaxy formation processes, such as feedback from supernovae and super-massive black holes, which can affect observables at similar levels to the modifications to GR. We find that it is possible to form rotationally supported disk galaxies in $f(R)$-gravity and that back-reaction effects on the matter power spectrum are small. We also find that the stellar and neutral hydrogen power spectra are suppressed in $f(R)$-gravity to a level that is detectable with future neutral hydrogen surveys.

The results presented in this work were derived from the SHYBONE (Simulating HYdrodynamics BeyONd Einstein) simulations, a set of full-physics hydrodynamical simulations employing the Illustris-TNG model in Hu-Sawicki $f(R)$-gravity (buchdahl1970, ; husa2007, ), which is one of the most widely-studied models of modified gravity. The theory introduces an extra scalar degree of freedom which mediates a fifth force between matter particles. This enhances standard gravity by up to $4/3$ in low-density environments, while dense regions are screened from such enhancements via the chameleon screening mechanism (khoury2004, ). This way, the model can pass the stringent constraints on gravity in the solar system (will2014, ) and still leave detectable signatures on large scales, making it an excellent tool to explore how possible deviations from GR would be observable.

$f(R)$-gravity predicts identical propagation speeds of gravitational waves and photons, and so is compatible with the recent gravitational wave detection with optical counterparts (lombriser2016, ; sakstein2017, ; ezquiaga2017, ; lombriser2017, ). With appropriate choices of model parameters (see Methods), it can explain the late-time accelerated expansion of our Universe without an explicit cosmological constant (husa2007, ) (Note that the value of $\Omega_{\Lambda}$ still enters the theory as a parameter). In addition to being theoretically well understood sofa2010 ; joyce2015 , $f(R)$-gravity has also been studied extensively in dark matter only (DM-only) (schmidt2010, ; zhao2011, ; li2011, ; lombriser2013, ; puchwein2013, ; hellwing2013, ; zivick2015, ; cautun2018, ; mitchell2018, ; li2018, ) and non-radiative hydrodynamical cosmological simulations (arnold2014, ; arnold2015, ; hammami2015, ; he2016, ). These works were nevertheless unable to include a suitably calibrated baryonic feedback model at the same time, which is nontrivial and computationally expensive.

Such simulations with feedback have recently become possible thanks to a modified arepo (springel2010, ) code, which employs a new and optimized method to solve the fully nonlinear $f(R)$-gravity equations in the quasi-static limit (see Methods). In this work, this more efficient code is combined with the Illustris-TNG galaxy formation model (pillepich2018b, ; springel2018, ; genel2018, ; marinacci2018, ; nelson2018, ) which incorporates prescriptions for gas-hydrodynamics, magnetic fields, star and black hole formation, feedback from supernovae and AGN, gas heating and cooling processes, as well as galactic winds. This new code made the SHYBONE simulations possible, which are dedicated to in-depth studies of the interplay of baryonic physics and modified gravity for the first time. A detailed description of the simulations used in our analysis can be found in Methods.

I Results

I.1 Galaxies in $f(R)$-gravity

Figure 1 shows a selection of galaxies from our full-physics simulations for the F6 (a and b) and F5 model (c and d) in the (total) mass range $1.5\times 10^{12}<M_{\rm 200c}/M_{\odot}<5\times 10^{12}$. Three of the objects are at $z=0$ (a-c), while the galaxy in Panel d is at $z=0.8$ (only one disk galaxy is identified at $z=0$ for F5). The galaxies were selected to have a rotationally supported stellar disk, i.e., a rotational-to-total kinetic energy parameter $\kappa=\sum(j_{z}/r)^{2}/\sum v^{2}>0.57$ (where the $z$-direction is defined by the total angular momentum of the stars in the galaxy) (ferrero2017, ).

It is evident from the central sub-panels of Figure 1 that both galaxies shown for the F6 model are partially screened. The gas in the inner regions is screened and experiences GR-like forces (dark blue regions; $a_{f(R)}/a_{\rm GR}\approx 1$), while the gravitational force experienced by the gas in the outer regions is enhanced by $4/3$ (bright yellow regions in the gas). The size of the central screened region depends on the total mass of the galaxy and its host halo, and the results imply that objects similar to the Milky Way Galaxy are partially screened for F6. In the F5 model, in contrast, the chameleon screening is less efficient, and the objects shown in the bottom panels show either no (c) or very small (d) screened regions in the central panels. Gas and stars will thus experience fully enhanced gravity across the whole object, implying that the model is unlikely to pass local tests of gravity.

Interestingly, this complicated force morphology in F6 does not hinder galaxies from forming rotationally-supported disks within the gas and the stars. We find that there are more galaxies with $\kappa>0.57$ in F6 compared to GR and significantly fewer in F5. The latter is an interesting observation that might be related to the stability of the disks in the unscreened regime or an enhanced galaxy merger rate. However, we caution that the number of disk galaxies is small in all of our simulations ($7,16,1$ objects at $z=0$ and $34,48,22$ objects at $z=0.8$ for GR, F6 and F5, in the large box, respectively and $56,64$ for GR and F6 at $z=1$ in the small box). The number counts of disk galaxies is therefore to be verified in future simulations.

I.2 The matter power spectrum degeneracy

It is well-known that $f(R)$-gravity can cause enhancements in the matter power spectrum while baryonic effects, predominantly AGN feedback, act in the opposite direction puchwein2013 ; arnold2018 . Previous simulation works were nevertheless unable to study both effects simultaneously and therefore could not model the possible back-reactions between them – a challenge that the full-physics simulations in this paper aim to tackle.

To inspect this degeneracy, in Figure 2 we plot the fractional changes of the total matter power spectrum with respect to results from the GR DM-only simulation. The dashed lines show the enhancement predicted by the $f(R)$ DM-only simulations. Solid lines show the combined effects predicted by the full-physics $f(R)$ simulations (for GR they show the impact of baryons only). The dotted lines are produced by simply adding up the relative changes due to modified gravity, as predicted by $f(R)$ DM-only simulations, and due to AGN feedback, as predicted by the GR full-physics simulation.

As one can see from the plot, baryonic feedback suppresses power on scales $k\gtrsim 2\,h{\rm Mpc}^{-1}$, with the suppression reaching $\sim 20\%$ at scale $k\sim 10\,h{\rm Mpc}^{-1}$. These results are in excellent agreement with previous works employing the same hydrodynamical model (springel2018 ) and from the EAGLE simulation (hellwing2016 ). Similar suppression at $k\gtrsim 2\,h{\rm Mpc}^{-1}$ due to baryonic feedback can be seen for F6 and F5.

It is obvious from the plot that the effect of baryons on the power spectrum for the F6 model is negligible for $k<2\,h{\rm Mpc}^{-1}$. The relative difference of the full-physics $f(R)$-gravity power spectrum is thus dominated by the modified gravity effect which increases power by a few percent at these scales. As for GR, the suppression due to baryons becomes active beyond $k=2\,h{\rm Mpc}^{-1}$. Due to the $f(R)$-gravity effects, the minimum of the relative difference is nevertheless only $12\%$ below the baseline. Interestingly, the blue dotted line in Figure 2 shows that the power spectrum of the full-physics F6 simulation can be predicted almost perfectly by an additive combination of the effects due to modified gravity and baryonic feedback. This suggests that there is negligible back-reaction between these effects for F6, and therefore computationally costly full-physics simulations for this model are not necessary as far as the matter power spectrum is concerned. This simplification does not hold for the F5 model, where the back-reaction effect is strong enough that a full-physics simulation is necessary if percent-accuracy predictions of the matter power spectrum are needed even at $k\lesssim 1\,h{\rm Mpc}^{-1}$. This will have important implications on the tests of this and similar models in surveys like Euclid.

The different strengths of interplay between modified gravity and AGN feedback in F6 and F5 can be explained by the chameleon screening mechanism. In the F6 model, the central regions of halos which are massive enough to carry AGN are screened, so that the flow of matter onto the central black hole (the AGN accretion) is not altered by modified gravity, resulting in similar AGN-feedback effects on the power spectrum as in GR. For F5, on the other hand, many of the AGN-hosting halos have become completely unscreened. The flow of matter towards central black holes therefore takes place in environments with by $4/3$ enhanced gravity, leading to modified accretion and feedback efficiency, and in turn a non-negligible back-reaction between the modified gravity and feedback effects in the matter power spectrum (for details on the different types of AGN feedback see pillepich2018 ).

I.3 Power spectra of individual components

The effect on the total matter power spectrum seen above is a combination of effects on the different matter components. In Figure 3 we display the power spectra of these individual components for $z=2,1$ and [math] (Panel a, b and c, respectively). The absolute spectra, which we find to agree with the original Illustris-TNG simulation (springel2018, ), are displayed in the upper panels, the relative differences relative to GR in the lower panels.

The components react in very different ways to modified gravity. The dark matter distribution is to first order only affected by gravity, its power spectrum is thus enhanced in a similar way as in DM-only simulations (e.g., winther2015, ; arnold2018, ). The dark matter and total matter power spectra are converged at $k\lesssim 20\,h{\rm Mpc}^{-1}$ between the $25\,h^{-1}\textrm{Mpc}$ and $62\,h^{-1}\textrm{Mpc}$ boxes, suggesting that the results in Figure 2 can be trusted in this $k$-range.

The gas power spectrum is only mildly influenced by modified gravity at $z=2$, and the effects are stronger at lower $z$. In particular, the gas is more clustered in F6 and F5 than in GR at intermediate ($k<10\,h{\rm Mpc}^{-1}$) scales. These effects are caused by the interplay of a number of processes that we only attempt to explain to first order here: the enhancement in gas clustering on intermediate scales is mainly a result of higher gas densities in the unscreened parts of galaxies – as Figure 1 shows, this will be primarily the case in low-mass unscreened objects or the outer regions of intermediate-mass objects. Comparing the enhancements in the gas power spectrum for F6 at $z=1$ for the small and the large simulation boxes, we see that the results are consistent at few-percent level for $k<10\,h{\rm Mpc}^{-1}$.

Unlike dark matter and gas, the clustering of stars is more strongly effected by modified gravity at high redshift. At $z=2$, the stellar power spectrum from the $62\,h^{-1}\textrm{Mpc}$ box is suppressed by $\sim 20$-$25\%$ in F6 and $\sim 30$-$40\%$ in F5 on scales $k\sim 1$-$10\,h{\rm Mpc}^{-1}$ (thick magenta lines). The predicted suppression for F6 from the $25\,h^{-1}\textrm{Mpc}$ simulation (thin magenta line) is weaker, which is $\sim 10$-$15\%$ in the same $k$ range, indicating that the clustering of stars has not converged in the big box at this redshift. At $z=1$, the two simulation boxes agree to a few percent for $k<5\,h{\rm Mpc}^{-1}$.

A main reason for the suppression of the stellar power spectrum is the gas density in low-mass objects: as smaller halos are already unscreened at high redshift, the gas within these objects is denser and can cool more effectively, compared to the similar objects in GR. This leads to enhanced star formation inside small halos. We have checked that the fraction of $m=10^{10}M_{\odot}$ halos carrying stars is in fact $60\%$ higher in F6 and $100\%$ higher in F5 compared to GR at $z=2$. These halos form from lower initial density peaks, which are intrinsically less clustered. This, together with the larger numbers of star-forming halos at high redshift in $f(R)$-gravity compared to GR, lead to a suppression of the stellar power spectrum in the former. As this process does not involve stellar dynamics, we expect it to be unaffected by potential self-screening of stars (which is not taken into account in our simulations).

I.4 The neutral hydrogen power spectrum

The enhanced gas cooling efficiency in low-mass objects in $f(R)$ gravity at early times also affects the distribution of neutral hydrogen (HI), a large fraction of which resides in small halos. We quantify this using the HI power spectrum, calculated following the procedure outlined in navarro2018 . The results from the $25\,h^{-1}\textrm{Mpc}$ boxes are plotted in Figure 3(a), which shows that at $z=2$ the HI power spectrum on scales $k\sim 0.02$-$10\,h{\rm Mpc}^{-1}$ is suppressed by $\sim 15\%$ for F6 (thin cyan line). The $62\,h^{-1}\textrm{Mpc}$ simulation boxes have not yet reached convergence as they cannot resolve small halos ($m\lesssim 10^{10}M_{\odot}$) which host a significant fraction of the total HI. The small box, on the other hand, resolves halos down to $m\sim 3\times 10^{8}M_{\odot}$, which host the vast majority of HI (navarro2018, ). For both box sizes we checked that the relative suppression of the HI power spectrum in $f(R)$ gravity is almost identical to the suppression of the power spectrum of HI-hosting halos at $k\lesssim 2\,h{\rm Mpc}^{-1}$, showing that it is more sensitive to the clustering of halos than to baryonic feedback. Comparing our results to the expected uncertainty for an intensity mapping experiment similar to SKA1-MID (a future $21\,{\rm cm}$ intensity mapping survey, (ska2015, )), we find that with 1000 hours of observing time, SKA1-MID should be able to distinguish F6 and GR for $z=2$. This is different from most other cosmological tests of this model, which are usually at lower redshifts.

Also interestingly, at $z=2$ the HI power spectrum for F5 is suppressed less than F6, which is because at this redshift F6 has a larger enhancement of the number of small halos while in F5 such small halos have merged to form larger ones. We checked that the F5 HI power spectrum is suppressed more at even higher redshifts.

While both the HI fraction and star formation rate in low-mass objects are sensitive to (stellar) wind feedback, we expect the impact by $f(R)$-gravity to be relatively robust against changes in the feedback model. This is because such changes cannot be arbitrary as that would otherwise lead to tensions between the simulations and observational data for the stellar observables (see Figure B1 in pillepich2018 for a detailed discussion on how different wind feedback implementations affect the stellar properties of galaxies at $z=0$). It is also known altay2013 that details in wind feedback are important in the strong Damped Lyman-$\alpha$ regime where HI column density is $>10^{21.5}{\rm cm}^{-2}$, while we have checked explicitly that modified gravity effects are non-negligible only at column densities $\lesssim 10^{21.2}{\rm cm}^{-2}$.

I.5 Galaxy properties

The SHYBONE simulations use the same hydrodynamical model for all three gravity theories. To check that a re-tuning of this model is not necessary for $f(R)$ gravity, in Figure 4 we plot a few stellar observables which were used in the tuning of the Illustris-TNG model. Panel a displays the stellar mass fraction (SMF) at the centres of galaxies, as a function of the host halo mass from the $62\,h^{-1}\textrm{Mpc}$ box. We find that $f(R)$-gravity introduces a $\lesssim 25\%$ change with respect to GR for both F6 and F5, which is smaller than the $1\sigma$ errors of the observational data (behroozi2013, ) shown for comparison. The absolute value of the SMF in our simulations does – as expected – not match the observations: The TNG model was tuned for a $25\,h^{-1}\textrm{Mpc}$ box with $2\times 512^{3}$ initial resolution elements, while simulations at lower resolution show lower star formation rates (pillepich2018, , see also the cyan lines in Panel d). For comparison, we also show the SMF from TNG simulations at several resolutions (the cyan lines in Panel a). The mass resolution of our $62\,h^{-1}\textrm{Mpc}$ box is close to that in the $25\,h^{-1}\textrm{Mpc}$ TNG box with $2\times 128^{3}$ resolution elements, and so are the SMF results. Supplementary Material Figure 1 further confirms that the SMF from our $25\,h^{-1}\textrm{Mpc}$ box, which has the same resolution as the highest-resolution TNG box, agrees with observational data at $z=1$.

The galaxy stellar mass functions measured within $30\,\rm{kpc}$ from the halo centres are shown in panel b. Again, the results from the $62\,h^{-1}\textrm{Mpc}$ box fall in between the TNG runs with $2\times 128^{3}$ and $2\times 256^{3}$ resolution elements, due to the reduced star formation rate in low-resolution simulations. The relative differences between $f(R)$-gravity and GR are smaller than the uncertainties in the observational data (grey lines).

Panel c shows the galaxy gas fraction as a function of the host halo mass. Unlike the stellar properties, the gas fraction is not resolution dependent. The results from the simulations presented in this work are therefore consistent with the three $25\,h^{-1}\textrm{Mpc}$ TNG boxes. The differences between the gravity models are again much smaller than the spread in observational data from giodini2009 ; lovisari2015 at $M_{200c}\gtrsim 10^{13}M_{\odot}$.

Finally, the cosmic star formation rate densities (SFRD) as a function of redshift are shown in panel d for the $62\,h^{-1}\textrm{Mpc}$ and $25\,h^{-1}\textrm{Mpc}$ (for $z>1$) boxes. The results from the three TNG boxes confirm that the SFRD is resolution dependent. Our $25\,h^{-1}\textrm{Mpc}$ simulations match the observational data of behroozi2013 , and are comparable to the TNG test-box with the same resolution. The $62\,h^{-1}\textrm{Mpc}$ boxes, as expected, show a lower SFRD which falls between the two TNG test-boxes with $2\times 128^{3}$ and $2\times 256^{3}$ resolution elements. The relative difference between F6 and GR is nevertheless converged.

The effects of $f(R)$-gravity on galactic properties are the consequence of the interplay of different physical processes. The star formation rate depends on gas density and temperature, with an enhanced gravity resulting in higher gas densities within galaxies and consequently a higher SFRD, as shown in Figure 4 d. In particular, at high redshift the star formation rates in $f(R)$-gravity and GR are roughly the same because most objects are screened; towards lower redshifts, galaxies gradually become unscreened and the SFRD in the two $f(R)$ models is enhanced with respect to GR. Due to the less efficient chameleon screening, in F5 this process happens earlier than in F6. At even later times, the relative enhancement of the SFRD in $f(R)$-gravity starts to decrease and becomes negative for F5 at $z\lesssim 0.7$. A reason for this is that the cold dense gas gets consumed faster in the $f(R)$ models, leaving less material for star formation at later times. The increased SFRD causes a higher stellar mass function (Panel b), as well as making a larger number of small halos star forming, explaining the weaker clustering of stars.

In addition to the four stellar and gas properties shown in Figure 4, we have also checked our simulation predictions of the galaxy size and black hole mass, which were used to calibrate the Illustris-TNG baryonic model too. In all these cases, the relative difference induced by $f(R)$-gravity does not lead to tensions with observational data, and so we conclude that a re-tuning of the baryonic model is not necessary.

II Discussion

Studying galaxy formation and its application to test gravity models with full-physics hydrodynamical simulations used to be a formidable task, because of the large computational cost and the complexities in calibrating a realistic baryonic model. There were also doubts regarding whether meaningful cosmological signals can be cleanly extracted given the uncertainty in the subgrid physics model. In this work, we present a new suite of high-resolution simulations of $f(R)$-gravity which incorporate a large number of baryonic processes to reproduce a realistic population of galaxies as implemented in the Illustris-TNG model. While this model has been applied in $f(R)$ simulations without changes, we have checked various galaxy gas and stellar properties, and found our simulation predictions in $f(R)$-gravity to agree with observations and the TNG simulations of similar resolutions, so that a re-tuning is not necessary.

A key question addressed here is the degeneracy between baryonic feedback and modified gravity, which has important implications for physical quantities, such as the matter power spectrum, that are central to future galaxy surveys such as Euclid. We find that in weak models like F6 these two effects can be modelled separately and combined additively, with better-than-percent accuracy in the predicted matter power spectrum at $k\lesssim 1\,h{\rm Mpc}^{-1}$, while for stronger models like F5 there is a back-reaction between them, necessitating costly full-physics simulations if sub-percent accuracy is desired.

Despite the complicated screening profile in them, we find rotationally-supported disk galaxies can form even in partially screened halos. There is an indication that strong $f(R)$ models such as F5 produce fewer disk galaxies than GR, which can be caused by the more frequent galaxy mergers in the former; if this is confirmed by future larger simulations, it may offer a possible test that is insensitive to baryonic physics (which has little impact on the merger rate). More interestingly, we find that the enhanced number of small dark matter halos in $f(R)$-gravity can significantly alter the spatial distributions of cool dense gas, stars and neutral hydrogen, and as a result the latter can be used as tracers for the clustering of small halos, which may offer useful model tests with future observations such as 21cm intensity mapping.

As a final remark, we note that though the resolution of the SHYBONE simulations is unprecedented in the case of $f(R)$-gravity, their relatively small boxes inevitably place a limit on a statistical analysis due to cosmic variance. Further, while the results look promising, in some cases, e.g., neutral hydrogen, research is still scarce even in GR, and in this sense we hope that our results will serve as a reality check and a call for more detailed studies in both GR and alternative models.

III Methods

The results presented in this work were derived from a suite of cosmological hydrodynamical simulations of $f(R)$-gravity which, for the first time, allow to study the interplay between detailed baryonic feedback and modified gravity in the same calculation. The simulations were carried out with the arepo (springel2010, ) code, by combining the Illustris-TNG galaxy formation model (pillepich2018b, ; springel2018, ; genel2018, ; marinacci2018, ; nelson2018, ) and a new, much more efficient, modified gravity solver in the code. Below we give an overview of the TNG model and $f(R)$-gravity, and introduce the new modified gravity solver. Further details about the code verification tests are included in the Supplementary Material.

The SHYBONE suite mainly consist of six simulations carried out for $\Lambda$CDM, F6 and F5 universes using the same initial condition. For each model, we carried out a full-physics hydrodynamical simulation and an associated DM-only run. These simulations initially contain $512^{3}$ dark matter particles and the same number of gas cells (for the hydrodynamical simulations only) in a $62\,h^{-1}\textrm{Mpc}$-a-side box, with mass resolutions $m_{\rm DM}=1.28\times 10^{8}h^{-1}M_{\odot}$ and $m_{\rm gas}\approx 2.5\times 10^{7}h^{-1}M_{\odot}$ (the mass within the gas cells can vary by a factor of $2$). For convergence test, two additional full-physics simulations (for F6 and GR) employing the same number of resolution elements in a $25\,h^{-1}\textrm{Mpc}$ box, with $m_{\rm DM}=8.39\times 10^{6}h^{-1}M_{\odot}$ and $m_{\rm gas}\approx 1.6\times 10^{6}h^{-1}M_{\odot}$, were run until $z=1$. All simulation use the same starting redshift, $z=127$. The softening lengths for DM particles and stars are $0.5$ and $1.25\,h^{-1}{\rm kpc}$ for the small and large boxes, respectively. The total run time is $\sim 2.5$ million core hours, and all simulations use Planck 2016 (planck2016, ) cosmological parameters $\sigma_{8}=0.8159$, $\Omega_{m}=0.3089$, $\Omega_{B}=0.0486$, $\Omega_{\Lambda}=0.6911$, $h=0.6774$ and $n_{\rm s}=0.9667$.

Unless stated otherwise, we measure total masses of halos in terms of $m_{200\,{\rm c}}$ which is the mass within a sphere enclosing $200$ times the critical density of the Universe.

III.1 Baryonic physics and feedback

The Illustris TNG model, a successor of the Illustris galaxy formation model (vogelsberger2014, ), incorporates prescriptions for a variety of physical processes necessary to reproduce a realistic galaxy population in cosmological simulations. The TNG model consists of a magneto-hydrodynamics solver on a moving mesh, a recipe for stellar formation, evolution and feedback, a prescription for black hole growth and AGN feedback (weinberger2017, ), as well as an algorithm for the chemical enrichment and cooling of the gas. It has been tuned to reproduce observational results on the stellar mass function, cosmic star formation rate density (SFRD), the relation between black hole mass and galaxy stellar mass, galaxy gas fraction, stellar mass fraction and the size of galaxies, in $\Lambda$CDM simulations (pillepich2018, ). For a more detailed overview of the TNG model, see pillepich2018b ; vogelsberger2014b .

For the simulations carried out for this project, we applied the same TNG model to all gravity theories considered. We are primarily interested in the question how astrophysical and cosmological observables in baryonic simulations are affected by modifications to gravity given a fixed model of galaxy formation. We explicitly checked that the changes to the observational gas and stellar properties against which the TNG model has been tuned are small compared to the uncertainties in the observations, so that a re-tuning is not necessary.

III.2 $f(R)$-gravity

$f(R)$-gravity (buchdahl1970, ) is an extension of Einstein’s general relativity (GR) including an additional, scalar degree of freedom which induces an enhancement to the gravitational force in low-density environments. In dense regions like the Solar system, the model employs the chameleon screening mechanism (khoury2004, ) to screen the modifications to gravity and recover GR-like behaviour.

The theory is constructed by adding a scalar function $f(R)$ of the Ricci scalar $R$ to the Einstein-Hilbert action of GR:

[TABLE]

where $G$ is the gravitational constant, $g$ is the determinant of the metric $g_{\mu\nu}$ and $\mathcal{L}_{m}$ is the matter Lagrangian density. Varying this action with respect to the metric leads to the Modified Einstein equations:

[TABLE]

where $G_{\mu\nu}$ and $R_{\mu\nu}$ denote the Einstein and Ricci tensor, respectively, $T_{\mu\nu}$ the energy momentum tensor and $\Box\equiv\nabla_{\nu}\nabla^{\nu}$, where Einstein summation is used. The derivative of the scalar function, $f_{\rm R}\equiv{\rm d}f(R)/{\rm d}R$, becomes a new degree of freedom whose dynamics is determined by the trace of the modified Einstein equation.

In cosmological simulations one often works in the weak-field and quasi-static limit (see sawicki2015 for a discussion of its validity for $f(R)$-gravity), where Eq. (2) simplifies considerably: similar to the Newtonian limit of GR, the gravitational potential $\Phi$ is given by a modified Poisson equation,

[TABLE]

where $\delta\rho$ is the matter density perturbation, and $\delta R=R-\bar{R}$ is the perturbation to the Ricci scalar $R(f_{R})$, with $f_{R}$ satisfying the following differential equation,

[TABLE]

Throughout this work we adopt the model proposed by husa2007 :

[TABLE]

$c$ denotes the speed of light here and $m$ is a mass scale of the model given by $m^{2}\equiv\Omega_{m}H_{0}^{2}$, with $H_{0}$ being the Hubble constant, and $c_{1}$, $c_{2}$ are model parameters. We choose $n=1$ in this work (This is somewhat arbitrary but the most convenient choice for the computations and widely used in simulations of $f(R)$-gravity). In addition to its relatively simple functional form, this model has two further advantages. First, it allows the so-called chameleon screening mechanism to suppress the modifications to GR in high-density environments, which is necessary to pass the very tight constraints on gravity within the Solar system (will2014, ). Second, the model features a cosmic expansion history which is very close to that of a $\Lambda$CDM universe if one chooses (husa2007, )

[TABLE]

Given the functional form in Eq. (5), under the conditions of Eq. (6) the scalar field $f_{R}$ can be approximated as

[TABLE]

and consequently $R$ can be written as a function of $f_{R}$. Instead of $m,c_{1},c_{2}$, the theory can be more conveniently described by the parameters $\Omega_{m},\Omega_{\Lambda}$ and $\bar{f}_{\rm R0}$, with $\bar{f}_{\rm R0}$ having the physical meaning of the present-day value of the background scalar field. With these parameters specified, the scalar field at arbitrary time can be expressed as

[TABLE]

where $R_{0}$ is the value of the Ricci scalar today and $R(a)$ is its value at scale factor $a$,

[TABLE]

The threshold for the onset of chameleon screening is set by $\bar{f}_{\rm R0}$. Current astrophysical constraints on the theory limit the background scalar field to values smaller than $|\bar{f}_{\rm R0}|=10^{-6}$ (F6); cosmological constraints are weaker (terukina2014, ) (he2018 find that the motions of galaxies in F6 predicted by DM-only simulations and sub-halo abundance matching (SHAM) are in tension with SDSS data. The authors of this paper nevertheless assume, that the same SHAM model can be applied to GR and f(R)-gravity if the mass of halos in the DM-only simulation, which is used as a parameter of the SHAM model in GR, is replaced by the dynamical mass for $f(R)$-gravity. The dynamical or effective mass is the mass a massive test particle would feel under the assumption of standard gravity. This mass measure accounts for the by $4/3$ enhanced gravitational forces in unscreened objects and is the same as the true mass in screened objects. This assumption is to be validated against full-physics MG simulations to asses the robustness of the constraint). We also consider a model with $|\bar{f}_{\rm R0}|=10^{-5}$ (F5) within this work as a toy model which features more prominent effects due to modified gravity in the simulations but is in tension with observations (terukina2014, ).

III.3 A new modified gravity solver in arepo

arepo (springel2010, ) is a massively-parallel, highly-optimized, cosmological simulation code. It employs a second-order Riemann solver using a moving Voronoi mesh to solve the Euler equations of hydrodynamics. The code also features state-of-the-art prescriptions of star formation, feedback from supernovae and active galactic nuclei, and magneto-hydrodynamics (vogelsberger2014, ; weinberger2017, ; pillepich2018, ).

The original gravity solver of arepo is based on the solver in p-gadget3, which is itself based on gadget2 (springel2005c, ). The solver employs an oct-tree algorithm to obtain the short-range gravitational forces while the long-range forces are calculated using a particle-mesh algorithm to reduce the communication load when running in parallel. arepo’s gravity solver has recently been improved and extended in several ways but for the implementation of the modified gravity solver we (currently) stick to the original version.

Because of the similarity between the gravity algorithms in p-gadget3 and arepo, we use the modified gravity solver in mg-gadget (puchwein2013 , also based on p-gadget3) as the starting point for the implementation of the $f(R)$ gravity solver in arepo. As the first step, we ported the solver from gadget to arepo, carrying out the necessary changes to adapt it to the code structure of arepo and optimizing the data structures of the solver to reduce its memory footprint. We will refer to this implementation as the ‘old’ method in the following. In order to improve the convergence behaviour of the multigrid solver, we implemented a ‘new’ method to solve for the scalar field based on the $f(R)$ gravity solver in ecosmog (li2012, ) developed by bose2017 .

In order to ensure that the new $f(R)$ solver implemented in the code gives the correct result, we performed extensive comparison tests between the old and the new methods, as well as convergence tests for the MG-timestep criterion and the residual threshold on the individual grid levels. A selection of these tests are presented in the Supplementary Material.

III.3.1 AMR-grid construction and mass assignment

In order to solve for the gravitational forces one has to solve Eq. (4) first before the scalar field can be used to calculate the gravitational potential through Eq. (3). Due to the chameleon screening mechanism, which causes the scalar field value to drop by many orders of magnitude in dense environments with respect to the background value $\bar{f}_{\rm R}(a)$, Eq. (4) is highly nonlinear. The commonly-adopted method in cosmological simulation codes is therefore to solve the equation iteratively on a grid (li2012, ; puchwein2013, ; llinares2014, ). All modified gravity codes optimized for cosmological simulations employ an adaptively refining mesh (AMR-grid) for this purpose to allow for high spatial resolution in high-density regions while keeping the total computational cost down to a realistic level (winther2015, ).

In arepo we use an AMR-grid to solve for $f_{\rm R}$ as well. This is constructed from the gravitational oct-tree of the standard gravity solver in the same way as in mg-gadget. Up to a maximum level, each tree node represents a mesh cell. In order to avoid a mesh that is too sparse on the finer AMR levels, all eight daughter tree nodes are created whenever a father node contains more than one simulation particle (this differs from the standard tree construction in arepo; for more details see puchwein2013 ). Because the tree is constructed from scratch for every time step in arepo, there is no need to take care of tree-updates for this implementation. It is nevertheless necessary to construct the tree and solve for $f_{\rm R}$ before computing the long range gravity forces on the PM-grid as opposed to the standard gravity implementation in the code, where the tree is stored during the short-range force calculation only (springel2010, ).

As in mg-gadget, the particle masses are assigned to the mesh using a cloud-in-cell (CIC) assignment algorithm allowing to calculate the density contrast field for each mesh cell. The following step, to actually solve for the scalar field $f_{\rm R}$, is where the old and new methods differ.

III.3.2 Scalar field solver

In the old method, Eq. (4) is rewritten in the following form by defining a new variable $u=\ln\left(f_{\rm R}/\bar{f}_{\rm R}(a)\right)$,

[TABLE]

which is solved using the Newton-Raphson relaxation scheme with red black sweeps and multigrid acceleration (see (puchwein2013, ) for more detail). This way unphysical positive values of $f_{\rm R}$ which might occur due to numerical errors during the iterations can be avoided.

The new method makes use of a different re-parametrisation of $f_{\rm R}$. Following bose2017 , we define

[TABLE]

which prevents positive values for $f_{\rm R}$ as well. This leads to a different form of the scalar field equation in terms of $u$,

[TABLE]

This equation is then discretised on the AMR-grid constructed before. On a mesh of cell side-length $l$, the Laplace operator becomes

[TABLE]

where we defined the operator $L$. Further defining

[TABLE]

and

[TABLE]

the discretised version of Eq. (12) becomes bose2017

[TABLE]

This cubic equation has a few advantages over the discretised equation in the old method. First, the result of the cubic equation can be found analytically for each mesh cell in each iteration step. The Newton-Raphson numerical update rule of the old method is thus avoided, leading to improved convergence speed and numerical stability of the solver. Second, the large number of expensive exponential and logarithm operations for each iteration step and mesh cell associated with the definition of $u$ are avoided. Finally, the new method does not introduce extra nonlinearities – through the exponentials in Eq. (10) – in screened regions which slowed down the convergence of the solver in the old method, particularly at high redshift.

Applying Cardano’s method in order to solve Eq. (17), one has to distinguish the three branches of solutions depending on the values of $p$ and $q$. In order to be consistent with (bose2017, ), we define $\Delta_{0}\equiv-3p$ and $\Delta_{1}\equiv 27q$. Independent of the density within a given cell and the field value of its neighbours, $q<0$ and thus $\Delta_{1}<0$ in all cases. The behaviour of the solution depends on the sign of the discriminant $D$:

[TABLE]

A: For the case $D>0$, Eq. (17) takes only one real solution which is given by

[TABLE]

We have used $\tilde{u}_{\rm i,j,k}$ to stress that the solution obtained is an approximation to the true solution $u_{\rm i,j,k}$ because the quantity $L_{\rm i,j,k}$ in $p$, cf. Eq. (15), depends on the values $u_{\rm i\pm 1,j\pm 1,k\pm 1}$ from the the neighbouring cells. Note that this branch of the solution includes the case $p=0$, for which Eq. (19) simplifies to $\tilde{u}_{\rm i,j,k}=\sqrt[3]{-q}$.

B: For $D<1$, Eq. (17) admits three real solutions:

[TABLE]

and $j\in\{0,1,2\}$. The cases $j=0$ and $j=2$ lead to $\tilde{u}_{\rm i,j,k}<0$ which is in contradiction with our definition of $u$. The only physical solution thus is $j=1$.

After the cubic equation is solved for a cell, the value of $u$ therein is replaced by the above solution. In order to allow for efficient parallelisation on distributed-memory compute clusters and to allow for a faster convergence, we carry out the cell updates using a red-black-sweep pattern. We start from an initial guess of $u=1$ for the first time step of the simulation and use the result of the previous time step as an initial guess for all subsequent steps.

At the end of each iteration, an approximate error for $u_{\rm i,j,k}^{2}$ ($\propto f_{\rm R}$) is calculated as

[TABLE]

where

[TABLE]

is the residual for a given approximate solution $\tilde{u}_{\rm i,j,k}$. We consider the solution on a given AMR-level to be converged and stop the iterations once

[TABLE]

This criterion translates to an (approximate) maximum error of $1\%$ on $f_{\rm R}/\bar{f}_{\rm R}(a)$ on a given grid level. We note that the resulting fifth force errors are smaller than the random standard gravity errors originating from the tree walk.

III.3.3 Multigrid acceleration

The relaxation method used above successively reduces the error of the initial guess by iterations. It is well known that the decay of long wave-length modes of the error can be dramatically sped up by multigrid methods. For both the old and the new method we therefore use multigrid acceleration arranged in terms of V-cycles to speed up the convergence of the scalar field solver. Our V-cycle implementation for the old method is identical to mg-gadget, and so we refer the reader to puchwein2013 for further details.

For the new method, we use the standard nonlinear full approximation storage (FAS) algorithm numerical . Eq. (17) is rewritten in terms of a non-linear operator

[TABLE]

The subscript l denotes quantities and operators defined on the ‘fine’ grid level with cells of side-length $l$ here, we will use L for quantities and operators defined on the ‘coarse’ level with cells of side-length $L=2l$. Eq. (23) can be written as

[TABLE]

for some approximate solution $\tilde{u}_{l}$. After a few pre-smoothing iterations to ensure that the short-wavelength modes in $\tilde{u}_{l}$ have been reduced sufficiently, both sides of Eq. (26) are mapped to the coarse level using the restriction operator

[TABLE]

where the summation is over the 8 daughter cells of the coarse cell. This leads to the following equation on the coarse level

[TABLE]

which can be solved either through iterations or further mapping to even coarser levels, to obtain an approximate coarse-level solution $\tilde{u}_{L}$. Using the prolongation operator $\mathcal{P}\tilde{u}_{L}=\tilde{u}_{L}$ the coarse-level solution can then be used to correct the fine-level solution as

[TABLE]

This correction is applied in our multigrid implementation as follows. If no initial guess for $u$ is available, the solver starts to solve Eq. (17) iteratively on level 3 (the top-level node covering the whole simulation volume is level 1), and maps down this solution as the initial guess for the next finer level. This process is continued recursively until the finest level that covers the whole simulation volume (max-level-full-tree, MLFT) is reached. If an initial guess is available, the solver starts on MLFT directly.

All the finer levels (including MLFT) are then solved by V-cycles using corrections from the two next coarser grid levels. The solver would, for example, repeat the cycles 7-6-5-6-7 on level 7 until convergence is reached and then continue on level 8 with 8-7-6-7-8 V-cycles using the solution of level 7 as the initial guess on level 8. During the V-cycles, the solver carries out 4 pre- and 4 post-smoothing iterations before mapping to and after mapping back from the coarser level, to ensure that the short-wavelength modes of the error of $\tilde{u}$ have sufficiently decayed. If the solver determines that the V-cycles do not lead to convergence on a given level, that level will be solved by iterations only. On the finest level which the mesh refinements extend to, the solution is obtained by iterations only for stability reasons.

III.3.4 Force calculation

Both the old and new methods adopt the effective mass approach to calculate the total gravity force as in mg-gadget, making use of the fact that the modified Poisson equation, (3), can be rewritten in terms of an effective density

[TABLE]

where

[TABLE]

The advantage of this approach is that it allows to calculate the gravitational accelerations from Eq.(30) using the standard oct-tree – particle-mesh gravity solver readily implemented in arepo.

III.3.5 Local time stepping

To keep the number of computationally expensive modified gravity force calculation low, we adopt the local time-stepping scheme that was first implemented in mg-gadget (for more details see arnold2016 ). This scheme takes advantage of the fact that the regions which require the smallest time steps, i.e., the regions with the highest accelerations, lie primarily within massive halos and are thus screened by the chameleon mechanism. For reasonable values of $f_{\rm R}$, the maximum modified gravity acceleration in the simulation box at a given time is orders of magnitude smaller than the standard gravity acceleration. Decoupling the modified gravity time step from the standard gravity and hydro time steps can therefore significantly reduce the computational cost of the simulations without loss of accuracy.

In practice, we decouple these time steps using a local time-stepping scheme already implemented in arepo. The modified gravity time step is coupled to the long-range gravity PM time step global, while the short range gravitational forces and hydrodynamics are calculated on smaller time steps for the particles individually. In order to avoid very large time steps for the modified gravity solver, the global modified gravity time step is calculated according to

[TABLE]

where $\Delta t_{\rm MG\,i}$ is the individual modified gravity time step for particle $i$ assessed using an acceleration criterion. The global PM+modified gravity time step is finally determined by

[TABLE]

III.4 Power spectrum calculation

The power spectra of the total density field and the density field of the matter components are calculated using a module implemented within the arepo code. The resolution elements are assigned to a Cartesian mesh using a CIC mass assignment in order to perform the required Fourier transforms (see springel2018, , for a detailed description).

Corresponding author

Please address any correspondence related to this work to Christian Arnold (christian.arnold@durham.ac.uk).

Acknowledgements

We would like to thank the Illustris-TNG collaboration for allowing us to use their baryonic model to carry out the simulations presented in this work. We are particularly grateful to Volker Springel and Rainer Weinberger for their help with the arepo code and for useful discussions on the results, and to Alejandro Benitez-Llambay for making Py-SPHViewer (llambay2015, ) available. Special thanks to Carlos Frenk and Jianhua He for their insightful comments on the results.

This work described in this paper is supported by the European Research Council through an ERC Starting Grant (ERC-StG-716532-PUNCA). BL is additionally supported by STFC Consolidated Grants ST/P000541/1 and ST/L00075X/1.

The cosmological simulations described in this work were run on the DiRAC Data Centric System at Durham University, United Kingdom, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.

The simulation code arepo (springel2010, ) is currently not publicly available. The analysis scripts used to analyze the simulation output can be made available to the reader on request.

Competing interests

The authors declare no competing interests.

Author contributions

C.A. and B.L. planned the project. C.A. developed, implemented and optimized (together with B.L.) the modified gravity solver AREPO, ran the simulations and performed main part of the analysis. M.L. performed the analysis for the HI power spectrum. C.A., B.L. and M.L. interpreted the results. C.A. wrote the manuscript with contributions from M.L. and B.L.

Supplementary Material

Supplementary material and the published version of this preprint are available from this website: https://www.nature.com/articles/s41550-019-0823-y