Spherical accretion of collisional gas in modified gravity I: self-similar solutions and a new cosmological hydrodynamical code

The spherical collapse scenario has great importance in cosmology since it captures several crucial aspects of structure formation. The presence of self-similar solutions in the Einstein-de Sitter (EdS) model greatly simplifies its analysis, making it a powerful tool to gain valuable insights into the real and more complicated physical processes involved in galaxy formation. While there has been a large body of research to incorporate various additional physical processes into spherical collapse, the effect of modified gravity (MG) models, which are popular alternatives to the $\Lambda CDM$ paradigm to explain the cosmic acceleration, is still not well understood in this scenario. In this paper, we study the spherical accretion of collisional gas in a particular MG model, which is a rare case that also admits self-similar solutions. The model displays interesting behaviours caused by the enhanced gravity and a screening mechanism. Despite the strong effects of MG, we find that its self-similar solution agrees well with that of the EdS model. These results are used to assess a new cosmological hydrodynamical code for spherical collapse simulations introduced here, which is based on the hyperbolic partial differential equation engine ExaHyPE 2. Its good agreement with the theoretical predictions confirms the reliability of this code in modelling astrophysical processes in spherical collapse. We will use this code to study the evolution of gas in more realistic MG models in future work.


INTRODUCTION
Spherical collapse is a widely studied phenomenon in cosmology. It describes the evolution of a spherically symmetric overdense region: how it slows down and decouples from the Hubble flow, turns around, and finally collapses into a singularity or some virialised matter distribution. Despite its simplicity, this scenario is of great importance, as it can describe several crucial aspects of structure formation of different matter components (e.g., collisionless dark matter and collisional baryonic gas), thus providing valuable insights into the real and more complicated cosmological process. Some cosmological hydrodynamical simulation codes also adopt this scenario as a test of reliability (e.g., , Teyssier 2002). The study of spherical collapse has a long history, with some of the early works including Gunn & Gott (1972); Fillmore & Goldreich (1984); Ryden & Gunn (1987); Subramanian et al. (2000). Among them, Bertschinger (1985) revealed an elegant self-similarity in the solution for a matter-dominated, Einstein-de Sitter (EdS), universe, for both collisionless and collisional matter. Using the turnaround radius, ta ( ), in the EdS model, the various quantities in the system of evolution equations can be rescaled, such that all the dependencies on the spherical radius and time are reduced into the dependence on a single variable ≡ / ta ( ). This gives a unique set of solutions of physical quantities, expressed in terms of , which can be used to ★ E-mail: han.zhang3@durham.ac.uk obtain the status of the evolution at arbitrary ( , ). Spherical collapse is therefore one of the few scenarios where a detailed semi-analytical solution is known in cosmology.
In the past decades, a lot of effort has been made to incorporate more physical processes into the spherical collapse model. Based on original radial collapse of matter, there are studies that look into the effects of angular momentum (Ryden 1988;Sikivie et al. 1997;Le Delliou & Henriksen 2003), dynamical friction (Antonuccio-Delogu & Colafrancesco 1994;Popolo 2009) and shears (Del Popolo et al. 2013;Pace et al. 2014). On the thermodynamics side, some research also studies cooling and heating process during the collapse (Abadi et al. 2000;Uchida & Yoshida 2004;McCarthy et al. 2007) for a more realistic thermal history of structure formation.
Nowadays, studies of cosmological structure formation have entered a highly advanced stage, with ever more complicated physical processes added into increasingly sophisticated hydrodynamical simulations (e,g, Schaye et al. 2015;McCarthy et al. 2017;Springel et al. 2018), which can realistically reproduce the observed properties of galaxies and clusters; see, e.g., Borgani & Kravtsov (2011), for a review. Comparatively, therefore, the role of spherical collapse as a stand-alone simulation experiment has declined. However, this scenario can still be served as an useful benchmark test to access the accuracy and reliability of new simulation codes. This explains partially why the original spherical collapse model with self-similarity still attracts attention (e.g., Halle et al. 2019;Alard 2020). It is also worth noting that the self-similarity can still hold under some other circumstances, if the physics added (e.g., the cooling function) follows certain assumptions (Sikivie et al. 1997;Uchida & Yoshida 2004).
This paper concerns spherical collapse in modified gravity (MG) models, which are an alternative solution to avoid several problems of the current concordance ΛCDM cosmological model. The ΛCDM model suggests that the majority of energy density in the Universe is cold dark matter (CDM), a species of non-baryonic and non-realistic particles, and dark energy, an energy component with exotic properties (e.g., negative pressure), in the form of a positive cosmological constant Λ (Amendola & Tsujikawa 2010), which is needed to explain the accelerated Hubble expansion. However, the hypothetical Λ suffers from long-standing theoretical problems (Weinberg 1989). Modified gravity models hope to overcome those issues by extending the standard General Relativity (GR) rather than assuming extra unobserved components (e.g., Sotiriou & Faraoni 2010;Linder 2010). In recent years, there has been growing interest in the MG models, because constraining them in various astrophysical and cosmological observations offers a powerful way to test our theory of gravity.
In this context, we want to mention that there are already various studies which look into the spherical collapse scenarios in different modified gravity models (e.g., Martino et al. 2009;Schmidt et al. 2010;Li & Efstathiou 2012;Lombriser et al. 2014; Barreira et al. 2014;Lopes et al. 2018;Contigiani et al. 2019). However, these studies generally focus on models that no longer uphold the property of self-similarity. This is not surprising, because even within GR there are strict conditions which must be satisfied to have self-similar solutions. For example, the EdS model loses its self-similarity property once a cosmological constant is added.
In this paper, we investigate the spherical collapse scenario for collisional gas in both Einstein-de Sitter universe and a slightly modified version of the Dvali-Gabadadze-Porrati (DGP, Dvali et al. 2000) braneworld model. The latter is a popular class of MG models that has attracted much attention in the last two decades, featuring an enhanced strength of the total gravitational force and the Vainshtein screening mechanism (Vainshtein 1972) which suppresses deviations from GR near massive objects to give the model a chance of passing the stringent Solar System and lab constraints. Despite its complexity, we find that the self-similarity property can still be achieved in this model under certain conditions which are not unnatural. Our self-similar solutions in this model will provide insights into how these mechanisms of modified gravity may affect structure formation in similar, but more realistic models where self-similarity no longer happens.
We then implement the spherical collapse scenario in a new cosmological hydrodynamical code, and use the above derived selfsimilar solutions to assess its ability in handling simulations for different gravity models. Our code is based upon the publicly-available hyperbolic partial differential equation (PDE) engine E H PE 2, which implements a blockstructured adaptive mesh refinement (AMR) (Dubey et al. 2016) Finite Volume (LeVeque 2002) code on spacetrees (Weinzierl 2019), and is parallelised through a combination of MPI, OpenMP BSP parallelism and a task formulation (Li et al. 2022). A comparison of the theoretical solutions and simulation results proves the reliability of our code, and we are going to use it to study more realistic and complicated modified gravity models in future work.
Our paper is organised as follows: in Section 2 we briefly introduce the DGP model ( § 2.1), and review the self-similar solution in an Einstein-de Sitter universe ( § 2.2.1) discovered by Bertschinger (1985). We then derive the self-similar solution for a slightly modified version of the DGP model, following a similar approach ( § 2.2.2), and compare the behaviour of the solutions in this model, for several difference parameter choices, with that of the EdS model ( § 2.3). This analysis reveals some interesting features of the solutions, which will be discussed in detail there. In Section 3 we describe our numerical code and the simulation configuration we use for the spherical collapse scenario, paying particular attention to the implementation details and certain tricky issues in the settings including the initial and boundary conditions ( § 3.4). The simulation results are presented and discussed in Section 4, and finally Section 5 is devoted to discussions and conclusions.

THEORIES
In this section, we introduce the physics we investigated and implemented in the code.
Throughout this paper, we assume that the background cosmology is that of the Einstein-de Sitter universe, i.e., a flat matter-dominated background(for simplicity we assume that this still holds even in the DGP models). It used to be the standard cosmological model before the ΛCDM model replaced it in the face of growing evidence that the cosmic expansion rate has been accelerating at late times, and it still serves as a good approximation for the real Universe between redshifts 300 and 2. The EdS universe assumes a zero cosmological constant and flat spatial curvature, and the equation of state of its non-relativistic matter content is ( ) = 0. With these parameters, the evolution of the scale factor of the universe, , can be derived analytically from the Friedmann equation as ( ) = 2/3 , where is the cosmic time, ≡ −2/3 0 is a constant and 0 is the cosmic time today (when = 1). This is an important assumption we will use to derive the self-similar solution later.

The DGP gravity model
The Dvali-Gabadadze-Porrati (DGP) braneworld model is a modified gravity model in a spacetime with an extra, fifth, dimension. The base assumption of this model is that the universe is a four-dimensional "brane" embedded in a five-dimensional spacetime, which is called a "bulk".
This model provides an explanation as why gravity is much weaker than other fundamental forces: all matter components are assumed to be confined on the brane, while gravitons could propagate through, or leak into, the extra spatial dimension.
The spacetime action of the DGP model is given by where is the Ricci scalar, is the determinant of the metric tensor, is Newton's constant, and the superscript (5) means the corresponding quantities live in the five-dimensional bulk. Others without it are normal four-dimensional quantities.
The modified Einstein equation for the DGP models can be derived from the variation of the gravitational action Eq. (1), which further leads to the following modified Friedmann equation that governs the cosmic expansion history ( ): where 0 = ( = 1) is the Hubble constant today (when the scale factor is = 1), Ω m0 is the present-day density parameter of matter (we have neglected the presence of radiation and massive neutrinos here since they are not relevant for the interest of this work), Ω DE ( ) represents the density parameter of a possible additional dark energy species at time , and Ω rc ≡ 2 /(4 2 0 2 ). Here, is the so-called crossover scale, which is a new free model parameter that indicates the scale above which the gravity begins to deviate from the standard Einsteinian: .
( 3) It is easy to see that, Eq.
(2) goes back to the usual form of the Friedmann equation when 0 → ∞. The ± in Eq.
(2) shows that this model has two branches of solutions. There is a "self-accelerating" branch (sDGP, the "+" branch) that can realise an accelerated Hubble expansion at late times without the need of a cosmological constant or dark energy, i.e., Ω DE ( ) = 0. However, this branch has several unsolved theoretical issues (Koyama 2007). Additionally, its predicted cosmological history is significantly different from that of ΛCDM and the observation also disfavours this model (e.g., Song et al. 2007).
The other branch, the so-called normal branch of DGP (nDGP) gravity, where Eq. (2) takes the "−" sign, can not provide an accelerated Hubble expansion by itself, and thus some additional dark energy component is needed (Ω DE ≠ 0) to explain the observation. This model has attracted much attention in recent years as it serves as a useful testbed of the Vainshtein screening mechanism (e.g., Brax 2013), despite its unappealing property of being still in need for additional dark energy. We will describe Vainshtein screening in more detail below.
The (modified) Poisson equation of DGP gravity and corresponding equations of the scalar field have been derived by Koyama & Silva (2007): and where Φ and are the gravitational potential and the scalar field of the model, respectively. They are also known as the brane-bending mode, which represents the position of the brane in the fifth dimension. ∇ is the spatial gradient (wrt to comoving coordinates), is the speed of light and m = m −¯m is the matter density perturbation (throughout this paper an overbar denotes the background value of a quantity). is a time-dependent function: for the two branches, which for the normal branch can be simplified as While we are interested in the DGP model, our main focus in this paper will be the effect of a fifth force that is mediated by the scalar field , denoted by the second term on the right-hand side of Eq. (4). To gain flexibility and to ensure self-similarity of the resulting model behaviour, we take the liberty to keep the main features of Eq. (5) but allow deviations from the exact behaviour of the sDGP or nDGP models. More explicitly, we will promote to a time-dependent function, and also allow to differ from Eq. (6). We remark that such variations from the original DGP model are not uncommon in other modified gravity models involving the Vainshtein mechanism, notably the cubic Galileon (Nicolis et al. 2009;Deffayet et al. 2009) and the Proca (Heisenberg 2014) theories.

Self-similar behaviour in collapse of collisional gas
In this subsection, we describe the self-similar collapse of collisional and non-radiative gas in some models. We first review the classic result from Bertschinger (1985), which applies to standard gravity in EdS universe. Then we proceed to show that self-similarity can also be achieved in the DGP model with Vainshtein screening. These can be used as a test case to verify our numerical implementation with E H PE 2 for both the standard and modified gravity scenarios, though our implementation of modified gravity is not restricted to the DGP model where self-similarity holds.

Einstein-de Sitter universe
Consider a uniform spherical overdensity region in the matter dominated universe background. Its initial condition could be written as < 0 is the initial cosmic time for this scenario to begin, = /¯ 1 the density contrast at , where¯=¯( ) and are respectively the mean matter density at time and the density perturbation, and is the initial radius of the spherical overdensity region. At the beginning, the Hubble flow is approximately unperturbed as 1. Thus, we have = and = 2/(3 ). As the universe expands, the matter inside starts to decelerate and decouple from the Hubble flow because of the slightly higher density. At some point it stops expanding completely (so-called "turnaround") and turns into a collapse. The turnaround for the mass shell at initially happens at a cosmic time and max radius (Bertschinger 1985 where the subscript ita stands for "initial turnaround". Matter inside the initial overdensity region starts to collapse first and all matter there infalls at the same time. No shell crossing happens. The matter initially in more distant shells (i.e., at initial radii > ) will start to collapse in progressively later times. The radius at which they turn around can be calculated using the Lagrangian picture. For the mass element initially located at , its evolution obeys the Newton's gravity law: Here accounts all mass interior to the shell we are considering. As no shell crossing happens during the evolution, it can be written as where = 1/(6 2 ) is the background density at for an Einsteinde Sitter universe. We then recast equation (10) using the following dimensionless time and radius variables, ≡ / and ≡ / , as d 2 d 2 = − 2 9 (1 + Δ) 1 2 .
Integrating this equation twice and using the assumption Δ 1, the solution can be expressed implicitly as (Bertschinger 1985 where we have defined the variables and for later use. As turnaround happens when reaches its maximum, this yields to ta = (where a subscript ta means "turnaround"). From Eq. (13) this corresponds to a time = (3 /4)Δ −3/2 and ta = ta / = Δ −1 . Combining these two expressions with the relationship between , Δ and given in Eq. (11), it is straightforward to derive the following expression of the turnaround radius: Now we switch to the fluid picture. The motion of a collisional gas in this system is governed by the gravity-driven Euler equations: where = ( , ), = ( , ) and = ( , ) are the density, velocity and pressure of the fluid at radius and time . ≡ (< ) represents the total mass within a given radius , and the adiabatic index. We now use Eq. (15) and define the new radial coordinate: as well as the dimensionless quantities , , and : where H = H ( ) is the critical density at time , which is equal to the mean matter density¯m ( ) in the EdS model. These allow us to cast Eq. (16 -19) as the following new dimensionless fluid equations (Bertschinger 1985): where a prime means the derivative wrt . Those equations only have one variable and thus could be solved directly given proper boundary conditions (see Section 2.3 below). No further time or length scales are involved, which means that the solutions to the system would remain identical throughout the evolution if expressed in terms of the coordinate. This is where self-similarity comes from. Obviously, if any new terms added in Eqs. (16 -19) depend on other scales besides , the solution to that new system will deviate from this self-similar solution.

DGP gravity model
In the spherically symmetric system, the scalar field equation (5) gets simplified significantly (see, e.g., Li et al. 2013): This equation does not contain scale factor as we use physics radius here. We then definê where we have usedˆto distinguish from ( ) introduced in Eq. (10), sinceˆdoes not account for the background matter density. Eq. (29) then can be integrated once to give: solving which gives the radial gradient of scalar field directly as where we have dropped the other branch of solution that is unphysical. The equation could be simplified further by defining the "Vainshtein radius" as follows: 3 = 16 2 9 2 2ˆ( ta ) .
Note that here we have usedˆwithin ta ( ) to define the Vainshtein radius, which differs from the usual definition that only accounts for the mass within the tophat radius -this is for convenience, because in this way we end up with a generic expression that does not depend on the particular size of any tophat. Now the gradient of scalar field reads as: Note that / determines the strength of the fifth force, and one can easily see the following limiting behaviour: If the scale of studied problem is significantly smaller than the Vainshtein radius , the gradient of the scalar field is also much smaller than that of the Newtonian potential, such that the fifth force is negligible compared with the standard Newtonian force. This is the idea behind the Vainshtein screening. Our next step is to try to recast the expression of the fifth force in the self-similar form (which, needless to say, is not always possible) similar to what we get above for the Einstein-de Sitter universe. This means that we hope that the ratio between the fifth force and standard Newtonian gravity, i.e., the coefficient in front of Eq. (34), depends on time and radius only through the combination ta ( ). We again define ≡ / ta ; note that this ta is the same as in Eq. (15)-this is mainly for convenience, but it does mean the ta in this expression is no longer the true turnaround radius in the DGP model. The mass can then be rewritten, using the definition of ( ) give in Eq. (24), aŝ As mentioned above, we have removed the contribution from the background mass as the fifth force only depends on density perturbations. The Vainshtein radius now reads as 3 = 16 2 9 2 2 2 9 3 ta and . (39) To achieve the self-similarity, we need to ensure that Eq. (40) only depends on . The dependence of 2 / 2 2 need to be cancelled out.
However, also appears in Eq. (34) in the overall factor 4/(3 ), and thus should be constant over time to avoid reintroducing an explicit dependency. This then leads to ∝ ∝ 3/2 , with the second proportionality true in an Einstein-de Sitter universe.
Denoting ( ) = 0 ( / 0 ), where 0 is the cosmic time today and 0 is the value of at 0 , and defining the dimensionless constant we get Therefore, the solution can be written as This expression shows that the fifth-force-to-Newtonian-gravity ratio can be written in a form that only depends on , which satisfies the requirement of self-similarity. It is straightforward to show that the coefficient of N in the above equation is always smaller than 2/(3 ), which means that the Vainshtein screening always works (though not necessarily always strong). Let us briefly comment that, according to its definition in Eq. (41), is the ratio between the crossover radius ( ) and . The latter can be considered as some characterisation of the size of the Einstein-de Sitter universe (actually it is 3 ). Therefore, the fact that this ratio is a constant in time implies that the Vainsthein screening mechanism is always effective on scales that correspond to a fixed fraction of the size of the universe, and therefore it should not be surprising that the self-similar properties of the EdS model have been preserved for this particular choice of ( ). Since characterises the length scale beyond which gravity is modified in the DGP model, we expect that for any physically interesting scenario we need to have ∼ O (1). The choice of = 2/3, for example, corresponds to 0 0 / = 1, which leads to a similar Vainshtein screening efficiency to that for a typical parameter choice in studies of the nDGP model for the same value of .
The actual strength of the fifth force is 1 2 , which means that the final expression for the fifth-force-to-Newtonian-gravity ratio is given by Turning to the derivation of the self-similar equations in the DGP model, i.e., the counterparts of Eqs. (25 -28), it is evident that only Eq. (26) needs to be modified. It is the only place where the law of gravity enters the calculation. However, instead of simply multiplying the − 2 9 2 by 1 + ( ), the correct final version of Eq. (26) is slightly more complicated. This is because ( ) is the ratio between the fifth force and N , which itself does not receive any contribution from the background matter density, c.f., Eq. (31). On the other hand, the term − 2 9 2 contains contributions from the background matter. Taking this into account leads to the following DGP version of Eq. (26): or equivalently A similar modification also appears in the DGP counterpart of Eq. (10), which now reads = − 2 9 (1 + Δ) where has been defined in Eq. (44), but is now expressed in terms of the dimensionless radius and time, and . More explicitly: This equation is needed for the exact solution of our equations in the next section. Before concluding this subsection, let us note that one limit of the DGP model arises from → 0, in which Eq. (43) approaches and so the fifth-force-to-Newtonian-gravity ratio approximately becomes 1/(3 ), which is the linear-regime (i.e., no screening) solution. This corresponds to a time-and scale-independent enhancement of Newton's constant by a factor of 1/(3 ) since we are assuming to be a constant here.

Self-similar solutions
Our next step is to find the exact solution to our self-similar equations Eqs. (16 -19): the profile of ( ), ( ), ( ) and ( ).
At the beginning stage, the spherical collapse can be described by a pressureless infall. Outside the radius of the tophat, the inner spherical shells infall at a greater speed than the outer shells, meaning that there is no shell-crossing or squeezing. However, when the infall speed of a given shell increases to a point where it exceeds the sound speed of the fluid, the shell impacts upon the fluid element inside it before there is enough time for the latter to adjust. A discontinuity of fluid properties, such as velocity, pressure and density, then starts to arise there, which is known as a shock. The shock location is our primary quantity of interest when we validate the outcome of our simulation. We assume the shock happens at radius or ≡ / ta (the subscript means shock), where we can apply the Rankine-Hugoniot jumping conditions, written in dimensionless forms: Here, a subscript 1 or 2 is used to denote the preshock and postshock values of a quantity, respectively, and is the dimensionless speed of the shock position itself. Physically, the three jumping conditions represent the continuity of mass, momentum and energy across the shock.
One can analytically calculate the preshock solutions in terms of using Eq. (12) and its solutions, Eqs. (13, 14) for Δ 1: are the values of and at , and, ≡ 1 − 3 2 . Combining Eqs. (50 -56), we get the boundary condition for the other side (post side) of the shock: The entire postshock solution can then be obtained by numerically integrating Eqs. (25 -28) inwards from = , using these boundary conditions. However, since is not known a priori, this is a trial and error process where the value of is updated iteratively until when the corresponding solutions meet the following physical boundary conditions in the centre of the system: This is how Bertschinger (1985) got his self-similar solution and we plot our reproduced result here in Figure 1. While the use of the variable to write the solution to Eq. (12) in the implicit forms of Eqs. (13, 14) is convenient, this is impossible for the DGP model where the corresponding spherical collapse equation takes a more complicated form. However, the introduction of in the EdS model is largely a matter of choice for convenience, and the same physics can be produced using as well. Because this is what we shall use for the DGP model, we decide to also use instead of to obtain the numerical self-similar solutions for the EdS model. This means that we need to express the preshock solutions to , , and at . For the velocity, using its definition we obtain where an overdot denotes the derivative wrt , and = ( ). For , using For , using For , we have 1 = 0 again. The following steps are the same as before. We can use the Rankine-Hugoniot jumping conditions to obtain 2 , 2 , 2 and 2 , and numerically integrate the equations again to find the postshock solutions. This time we need to vary for our trial-and-error process after Δ is specified. and can be calculated numerically from by using Eq. (10) for the EdS model and Eq. (47) for the DGP model. For EdS, we have explicitly checked that using the -based approach to set up the boundary conditions for the postshock solutions gives identical answer as using Eqs. (53 -56), as expected. We summarise our result for self-similar solutions in DGP gravity with in Figure 1. The black curves in the figure are the self-similar solutions to , , and for the EdS model, which we find to be in excellent agreement with literature results (e.g., Bertschinger 1985). The coloured curves show the results for several variants of the 'artificial' DGP model described in Section 2.2.2, with the case = 0 (red) corresponding to a constant enhancement of by 1/(3 ). The cases with = 1, 5 and 10 represent progressively more efficient Vainshtein screening, which explains why they are in between the EdS and = 0 cases. In particular, we see that at = 10 the screening is already very efficient so that the brown curves are very close to EdS. The qualitative trend also agrees with what one should expect for a model with enhanced gravity: the infall becomes faster such that the preshock solution of becomes more negative and the shock happens at larger radius; the density and pressure are also higher due to the stronger structure formation, and the latter explains why the enclosed mass within a given radius is larger.
The DGP results in Figure 1 are obtained with the parameter Δ = 0.001. The Δ dependence of the solution, which we have explicitly checked, can already be seen at the equation level, cf. Eqs. (62-64), and also in Fig. 2, where we show the two sets of self-similar solutions for two models, EdS (dashed lines) and DGP with = 1.0, = 1 (solid lines). Different colours indicates the value of Δ for each curve, as shown by the legends. We can see that the DGP model has a very strong Δ dependence. A similar dependence is also present in the EdS case, but is much weaker there -indeed, it is known that in EdS there is approximately no Δ dependence in the limit Δ 1 (Bertschinger 1985). This Δ dependence comes from our choice of using the turnaround radius ta ( ), given in Eq. (15) dimensionless coordinate , where ta itself depends on Δ. For the EdS model, rescaling using this turnaround radius helps to cancel out the Δ dependence from Eq. (12), because this ta is calculated from the same dynamical equation and has the physical meaning of where the shell start to collapse. But such a cancellation should not be expected to happen when we use the same Eq. (15) to define for the DGP (and generally other gravity) models, since it does not represent the true turnaround radius anymore.
Using the same ta to define in all models above certainly has its advantages. One of these is that Eq. (15) is an analytical function with a power-law dependence on , which is convenient when deriving the dimensionless equations governing the self-similar evolution. It also allows these equations to take the similar form between the DGP and EdS models. For example, Eq. (25-28) remain almost the same for the DGP model, with only some slight changes of Eq. (26) to Eq. (46). In addition, Figure 1 clearly shows the effect of modified gravity law on the collapse of collisional gas and on the formation of shock: this also benefits from the fact that we have used the same 'turnaround' radius, ta ( ), to define the rescaled quantities in all models, so that the differences in the rescaled quantities reflect the differences in the same quantities pre-rescaling. Nevertheless, for theoretical interest, we also want to see the results when we actually define using the true turnaround radius of each model. Because there are no analytical expressions for ta for the DGP model, this has to be done in a "post-processing" way: after getting the profile ( ) by following the above steps, we can obtain the real turnaround radius in the preshock ( ) solution, by looking for the value of ta where ( ta ) crosses 0; we then get the correct turnaround radius as: and use ta to rescale our solutions for the other quantities, which is equivalent to performing the following 're-rescaling': The new result is summarised in Figure 3. The Δ dependence in EdS case is negligible as the solutions assume the limit Δ 1 (Bertschinger 1985). On the other hand, the results for DGP model shows a clear dependence on Δ, which we explain in the text. the results obtained by using Δ = 0.001 here, we find that using Δ = 0.01, 0.1, 0.2 give very similar results. One notable property is that the new rescaled profiles are very close to that in Einstein-de Sitter universe, i.e., the DGP model behaves similarly to standard gravity if expressed in terms of the coordinate which is defined using the true turnaround radius of the model. As the real physical evolutions of these models are very different, this similarity is quite interesting, since it suggests that self-similarity works (at least to a good approximation) in more general models than just EdS.
As we shall see below, this "re-rescaling" idea using the true turnaround radius can also be applied to the numerical simulation result from E H PE 2, and help to check its reliability on handling this scenario.

NUMERICAL SIMULATIONS WITH E H PE 2
In this section, we first introduce the numerical code we implement on E H PE 2, then describe how we configure the spherical collapse scenario with it.
Our simulations are based upon an adaptive Cartesian mesh host-ing a Finite Volume discretisation with an explicit Euler. The code is realised through E H PE, which is a publicly available engine designed for generic hyperbolic PDEs that arise in different branches of sciences and engineering. We rely on the second-generation E -H PE 2 code which is a rewrite that has been used for astrophysical challenges before (e.g., Reinarz et al. 2020).

Spatial and temporal discretisation
E H PE 2 constructs the spatial discretisation from a spacetree formalism (Weinzierl 2019) combined with block-structured adaptive mesh refinement (AMR; Dubey et al. 2016): The computational domain is embedded in a cube and split into three equal parts along each coordinate axis. This yields 3 3 = 27 smaller cubes. We continue recursively, i.e., decide for each cube whether to cut it into 27 subcubes again. The refinement decision or criterion is subject of discussion below. The process yields an adaptive refined Cartesian mesh. Starting from an initial adaptive mesh, dynamic adaptivity could be realised by applying the splitting in between time steps to yield a finer mesh.
Each cube hosts a × × Cartesian mesh. We call these Cartesian meshes patches and make them carry the actual solution representation: each mesh element in the patch holds a piecewise constant solution of the governing equations, i.e., defines one "finite volume". Every patch thus consist of 3 volumes. The patch of volumes is augmented with a "halo 1 layer" of width one around it. The patches hence yield a non-overlapping domain decomposition of the computational domain, while the haloes introduce an overlap between them. Let the vector ì : R × R + ↦ → R 5 denote the unknowns of interest as they evolve over time, where the symbol ì highlights that this is a data (rather than space) vector that in our case has a dimensionality of 5. We approximate the time derivatives with forward finite differences, i.e., d ì d ≈ ì new − ì old with a given time step size and ì old , ì new representing the values of ì at the start and end of the time step. Our equations are a set of generic first-order hyperbolic PDEs d ì d where ( ì ) and ( ì ) are the flux and source term, respectively. Here we have used bold symbols to denote space vectors to distinguish them from the notation for data vectors introduced above. The generic first-order hyperbolic PDEs can be written in a weak formulation for one timestep as  (LeVeque 2002) to solve the Riemann problem that arises once we assume that the solution remains constant within every timestep and every volume , and set all test function as characteristic function of one finite volume, i.e., they are ( , ) = 1 within and vanish anywhere else. The integration of Eq. (68) over time gives us where d is the (oriented) area element of the surface of the volume , . Here the closed-surface integration is decomposed into the summation of multiple faces that have constant normal vector respectively. At the same time, we assume the solution vector ì to be piece-wise constant, so we apply the following replacement: where is the volume of , , are the area and unit normal vector of one face of , respectively, and ( ) denotes a generic scalar (space vector) function. This leads to final explicit Euler time stepping scheme we implemented in the code: with the so-called Rusanov flux (Rusanov 1961): is the flux term evaluated along of the considered for the respective volume. Flux ± ( ì ) in Eq. (72) averages component-wisely over the flux within the two adjacent volumes. The average then is corrected (limited): max is the largest eigenvalue of the matrix ( ì ) acting on the gradient along if we write down the PDE along the face normal as d ì d It indicates the largest propagating speed of the quantities in the system. We close this subsection by briefly commenting that is subject to the Courant-Friedrichs-Lewy (CFL) condition with where < 1 is a problem-specified safety parameter. Our scheme employs a global time stepping scheme and thus uses the smallest global face length | |. It remains invariant over time as we fix the finest resolution in our simulations. The maximum eigenvalue max , however, changes over time and thus has to be recalculated after each step.

Implementation
Our code splits up the computational domain along the Peano spacefilling curve (SFC) into subdomains (Li et al. 2022; Weinzierl 2019) (Fig. 4): All patches are ordered along the SFC. We cut this sequence of patches into segments such that each rank gets exactly one segment hosting roughly the same number of patches. As the Peano SFC is continuous, the set of patches per rank form a connected subdomain of the computational domain which does not overlap with any subdomain handled on another rank. Per rank, we apply the SFC splitting once more such that each thread per rank obtains its own subdomain: The patches within the computational domain are first distributed among the ranks and each rank then distributes its patches once more among the threads. This gives us a two-level non-overlapping MPI+OpenMP parallelisation. Our realisation with patches supplemented with a halo of width one allows each thread to run through the mesh, and to update all of its patches independently of the other ones. After this mesh traversal, the halo layers are copied over for neighbouring patches of the same size, halo finite volumes overlapping with coarser resolution patches are updated due to a linear interpolation, while halo volumes overlapping with finer resolutions are updated through averaging over the finer volumes. We map the individual patch updates per thread onto a task formalism (Li et al. 2022) and process the patches along the MPI boundaries prior to other tasks such that the data transfer required for the halo updates can overlap with further computations (Charrier et al. 2020). We assume that the tasks can compensate for any geometric ill-balancing on the MPI level. We do not dynamically rebalance throughout the computation.
In our experiments, we use four nodes of Durham's COSMA 7 cluster with one MPI rank per compute node. Each node hosts a dualsocket Intel Gold 5120 CPU processor. Therefore, each rank splits up its domain into 28 further subdomains. Our experiments stick to = 3. While this setup yields a relatively low arithmetic load per patch compared to the overhead that we need to maintain the halo volumes, it ensures that we can use a rather aggressive coarsening towards the domain boundaries to reduce the overall computational burden.

Code Units
To solve the system of equations numerically, it is usually convenient to recast them by using dimensionless quantities. In the E H PE 2 implementation, we adopt the so-called supercomoving coordinates, which are used in other simulation codes such as (Teyssier 2002).
The original formulation of this coordinate system could be found in Martel & Shapiro (1998). Its idea is to apply the following rescaling of the variables: Here ,¯m ( ) are respectively the critical density today and mean density of matter at time ; is the comoving size of unit code length; d , and denote, respectively, the (physical) time interval, physical coordinate and peculiar velocity. We use the quantities with a tilde in our code, we therefore call them code unit in the following context. The supercomoving coordinate system factors out most of the effect from the Hubble expansion, and thus allows us to implement the original fluid equations Eq. (16-19) in a static space with just minor changes. For the special case = 5/3, the only change of the fluid equations is a re-calibration of the gravity term in Eq. (17), which now needs to be derived from the following code-unit Poisson equation: whereΦ is the Newtonian potential in code unit Solving Eq. (76) under spherical symmetry gives us the following solution of the Newtonian gravitational force˜≡ −dΦ/d˜(again, in code unit): where we have defined˜(<˜) to be the total "mass perturbation" within radius˜, i.e., the difference between the total mass therein and the mass in the same region were the density there equal to¯m. For other fluid equations, we only need to replace physical quantities with code quantities directly. For cases ≠ 5/3, extra terms are needed for supercomoving coordinates (although they are straightforward to derive), which we do not cover here.
The generalisation to calculate the modified gravitational force in the DGP model is straightforward: we multiply the fifth-force-to-Newtonian-gravity ratio given in Eq. (48) to Eq. (78) directly to obtain the fifth force in the DGP model. Most terms in Eq. (48) are constants or time-dependent functions, and the only term that needs to be rewritten in code unit is where we recall that is the initial radius of the fluid element located at at time , and ( , ) is the total mass enclosed within at the initial time . As no shell crossing happens during the evolution, the mass within the radius of this fluid element remains the same, which means: where (< , ) denotes the total mass enclosed in radius at time > . In our code implementation, the mass is calculated by counting volumes (see Section 3.4 below), and thus (< , ) = ≤ ( )ℓ 3 , where the subscript labels the volumes, ℓ is the cubic size of volumes , and ( ) is the density (all in physical units). Notice that we have: in the Einstein-de Sitter universe. Putting Eq. (81) back to Eq. (80), we get where in the second equality we have used 4 3 3 = ≤ ℓ 3 , while in the final equality we have replaced ℓ and with their code-unit expressions,l and˜, which does not change the ratio ℓ 3 / 3 , and used (<˜) ≡ ˜<˜(˜− 1)l 3 . Eq. (82) is the final code expression that we use in our simulation.

Simulation Settings
In this subsection, we discuss how we implement the spherical collapse scenario on E H PE 2. We describe the hyperbolic equations and grid setting that are used in the simulations, the initial conditions and boundary conditions, and how we calculate the total perturbed mass,˜(<˜), at arbitrary radius˜.

Equations and Grid setting
In the simulations, we implemented the original conservation form of the (gravity-driven) Euler equations in code unit: where˜,˜,˜,˜represent the density of mass, momentum, energy and pressure in code unit respectively,˜=˜˜is the force density with˜the gravitational acceleration, which is proportional to˜(< )/˜2. We consequently obtain ì = (˜,˜,˜) in Eq. (67). All simulations we presented in this paper use the same grid setup on a cubic box [−1.5, 1.5] 3 . The maximum refinement level within the tree formalism is 3, corresponding to a resolution of 243 3 patches on the finest level. Every patches contains 27 volumes again ( = 3). We coarsen this mesh once at a distance of 0.5 (in code units) away from the origin, and coarsen it once more at 0.7. Figure 5 illustrates the AMR refinement pattern we used for the simulation. The exact refinement pattern is chosen such that it covers the refinement radii. The safety parameter (CFL ratio) we use in Eq. (74) is = 0.3. The adaptive Cartesian grid used in our simulations, with patches and volumes that we describe in Section 3.1 therein. The patches with = 3 (i.e., every patch contains 3 3 volumes) are separated from each other in the visualisation with gaps for clarity. Three levels of the grid are shown here. Only one quarter of -plane taken from a slice of the simulation box perpendicular to the -axis is plotted. The diagonal lines are visualisation artefacts as we use the cubic finite volumes. The refinement transitions are conservative, i.e., they are slightly larger than the resolution transitions imposed by the refinement strategy. Right Panel: The density field (in code unit) on the same slice for a snapshot during a simulation. Some fluctuations of the density field could be seen out of the central peak, as we discuss them in Section 3.4.3.

Initial Conditions
The simulations shown in this paper start at scale factor = 0.001, and end around ≈ 0.3. The simulation domain is initially filled with collisional cold gas of = 5/3 in critical density (which is unity in code units). Our overdense seed, the spherical tophat, is placed at the origin and is set to have a radius˜= 0.05 and total perturbed mass = 0.15. The treatment of the initial conditions of the pressure, density and velocity is subtle. Although we should expect a pressureless infall for most regions in the simulation box at the beginning, we can not set a zero initial pressure numerically. Likewise, although it seems to be quite natural to set a zero initial velocity profile within our comoving coordinate system, we can not do this in our implementation, neither. Both of these would lead to a negative pressure in the first time step. This is because in this step the energy equation, Eq. (85), does not update the local energy given the zero momentum (i.e., both the flux and the source terms are zero in this equation). On the other hand, the momentum itself is updated normally according to Eq. (84) as its source term (the force density) is nonzero. Since we calculate the pressure using: the fact that˜is updated (mostly increased in magnitude) while˜is not during the first step can cause an accidental and unphysical drop of pressure at the end of this timestep, and frequently (for zero initial pressure, it is always) the pressure turns to be negative where gravity is strong, i.e., near the centre. This issue would be worse if we put a point mass as the overdense seed at the centre, like the one in (Teyssier 2002), because it leads to an extremely large magnitude of the gravity force in the adjacent volumes of the point mass.
To address this negative pressure issue, our solution is three-fold. Firstly, we stick to using a tophat overdensity rather than a point mass as our seed, though it harms the solution partially (see the section for results below). A tophat initial profile smooths the gravity field and reduces the magnitude of a potential negative pressure. Secondly, we set a very small but non-zero value for the pressure initially: it makes the system more robust to the pressure drop in the first time step, and can quickly converge to the correct pressureless solution outside the shock later in the simulation. Finally, we introduce a pre-set initial velocity profile. We assume our momentum field has evolved a small period of (physical) time before the simulation begins, according to the initial gravity field: such that the energy can get updated as well. These adaptions successfully solve the initial negative pressure issue without the explicit construction of consistent initial condition which does not yield unphysical solutions. The freedom of adjusting our initial conditions without harming the final self-similarity is expected given the convergence of the solution (Alard 2020), and we have explicitly checked that it is true for our simulation by tuning our initial pressure.

Boundary Conditions
Our setup to simulate spherical collapses requires free inflow boundary conditions. Because we expect ì to be almost stationary in comoving coordinates (or approaching the Hubble flow physically) as we move away from the centre of the computational domain, homogeneous Neumann boundary conditions can yield the free inflow as long as the computational domain is sufficiently large. However, such a large domain is computationally inefficient or even unfeasible, and it is also not clear whether 'large' is well-defined in an evolving system: the shock propagates outwards towards the border over time, thus making it a challenge to use homogeneous Neumann boundary conditions throughout the entire evolution. We therefore use the following hybrid scheme: where ì in and ì out denote, respectively, the solution vectors in the volumes on the inner and outer sides of the boundary (see Fig. 6). The boundary conditions in E H PE 2 are implemented by specifying how the quantities in ghost volumes out of the boundary ì out are calculated from ones in their direct neighbours within the domain ì in . In most times, we use the extrapolating boundary condition Eq. (89), where the superscript (1) means we use the first-order approximation of the gradient ∇ ì at approaching the domain boundary Ω, multiplied by the distance between the two volumes, ℓ in .
The different behaviours of these two types of boundary conditions are illustrated in Fig. 7. The linearly extrapolated boundary condition is more accurate than the homogeneous Neumann one specified by Eq. (88), but it underestimates the momentum inflow from beyond the boundary. As a result, the code-unit density at the boundary, in , will drop to under unity later in the evolution: this is unphysical because the density everywhere in this collapse scenario should be above the critical density. Whenever this happens, we switch to the homogeneous Neumann boundary condition, Eq. (88). The latter usually overestimates the inflow, and thus can provide some 'compensation'. After the density˜i n increases back to above unity, we continue using the extrapolating boundary condition again.
The Finite Volume scheme uses normal boundary conditions where the normal is axis-aligned. However, our solution is sphericalsymmetric. The boundary condition's normal alignment thus is erroneous. Even with Eqs. (88-89), we have to ensure that the domain remains sufficiently large compared to the area of interest, such that this misalignment becomes negligible. This naturally limits the maximum simulation time up to which our results are not polluted sig-nificantly by the tangential boundary errors, as the solution's steep gradient moves towards the domain boundary.
Similar arguments hold along resolution transitions. As we interpolate linearly along the resolution boundary, our solutions do not follow exact spherical symmetry: the mesh and its resolution transitions should be spherical, and we should interpolate linearly along a spherical transition. Yet, our grid is Cartesian. This "misalignment" results in fluctuations or finger patterns (Fig. 5, right panel). Our code has two ingredients to mitigate the resulting error: on the one hand, we use 2:1 balancing (Sundar et al. 2008), since a more aggressive resolution change would amplify any error. On the other hand, we ensure that the "first" (finest to second finest) resolution transition is sufficiently far away from the region of interest, i.e. the shock. In return, this implies that the maximum runtime yielding physically admissible results is bounded further, as long as we disable adaptive mesh refinement-a technique which is intrinsically limited, as the area of interest expands and thus eventually yields a regular grid with excessive memory footprint.

Mass Integration
Most of the terms in Eqs. (83-85) can be implemented in E H PE 2 directly as part of the Rusanov scheme on Cartesian meshes we describe above, because they are all localised variables, i.e. follow up the update pattern of any Finite Volume scheme. However, the gravitational force is not localised as we will need the total perturbed mass within radius . To get˜(<˜), we construct a mass array {˜} 0≤ ≤ max which stores the total perturbed mass values within radii {˜} 0≤ ≤ max . Herẽ max =˜= max is chosen to be the radius of the largest sphere in the simulation box: half of the domain length. The values of˜are calculated by accumulating the mass in all volumes that are withiñ per time step: wherel is the size of the accumulating volume. The plain summation is consistent with our choice of piece-wise constant Finite Volumes. During the subsequent time step, we apply the following interpolation rule per volume according to its radius˜for the required perturbed mass: The perturbed masses for volumes outside˜m ax are approximated by assuming that the density there is equal to that at˜m ax . During our simulations, the densities in those volumes depart little from unity and thus contribute little to the total perturbed mass. This approximation is therefore acceptable. More accurate schemes could be used in future simulations, such as using a scheme of density interpolation that can extend to the furthest corner of the simulation box. Within˜m ax , on the other hand, the accuracy of this interpolation rule depends on the size and arrangement of the sample array {˜}.
In our simulations, we use a sample array size of 200, and keep our sample radii {˜} invariant over time.  88), where the velocity field outside the boundary (i.e., to the right of the vertical line) is assumed to be a constant equal to the velocity value just inside the boundary (the red dashed line). In this case, the inflow from beyond the boundary is overestimated, and thus harms the quality of the boundary. Right panel: The first-order extrapolated boundary condition corresponding to Eq. (89), as indicated by the green dashed line. Its prediction of inflow is more accurate than the Neumann case but is underestimated. We combine these two boundary conditions in our simulations depending on the local density at the boundary.

SIMULATION RESULTS
In this section, we report the simulation results of spherical collapse scenarios in different gravity models using our new code. To make comparison to the theoretical predictions we got in Section 2, we will also show results that are rescaled following Eqs. (21-24), after we restored the quantities in physical unit using Eq. (75).

Einstein-de Sitter universe
We first show the simulation results in the Einstein-de Sitter universe.
Since gravity is standard, we can use Eq. (15) as the scaling radius. The rescaled profiles of physics quantities are plotted over the radius coordinates in the code unit (supercomoving coordinates), in Fig. 8. We illustrate five snapshots of the system (at scale factor ≈ 0.022, 0.031, 0.047, 0.076, 0.145) from the late part of the simulation when the corresponding Δ is relatively small. The system remains in stable evolution before the numerical issues we reported in the last section pollute the solution. A clear outward-propagating shock can be seen in the figure.
The same profiles of quantities are plotted again, but now against the rescaled radial coordinates , in Fig. 9. The theoretical self-similar lines from Section 2 are shown as black dashed lines for comparison. We can see a clear self-similarity here, as the rescaled simulated quantities have converged during the time period considered, when the scale factor increases by a factor of seven. The coloured vertical lines in the figures are the positions of the tophat edge at the time of the corresponding snapshots, within which the density and pressure solution deviate from the self-similar solution and flatten: this is expected as the gas within the tophat does not experience the full gravity from the mass perturbation anymore. The radius of this edge is shrinking in the rescaled plots over time because the turnaround radius that is used to define increases as time evolves.
The rescaled solutions agree with the theoretical predictions quite well, especially for the preshock solutions of the density and velocity. Yet, there are some deviations from the self-similar solution, notably a shift of the shock position. Because of this, the infall velocity of the gas just outside the shock is lower than the theoretical prediction. This is a common numerical artefact caused by volumes with finite widths, which cannot exactly resolve the infinitesimally thin shock. We have checked that the agreement with the self-similar prediction improves as we use finer volumes. A detailed convergence study is beyond scope here.
Another factor that may have contributed to the difference between theory and simulation is that the theoretical solution here is obtained under the assumption of Δ 1, and this is not well satisfied in the simulations. The different shells of gas have different initial radii and corresponding values of Δ, with the outer shells having larger and therefore smaller Δ, and vice versa. The outer shells also collapse to the shock at a later time. We note that the outer shells that collapse to the shock at later stages of the simulations usually have Δ ≈ 0.1, and the inner shells have even larger Δ. The difference between these values of Δ 1 might affect the accuracy of the simulation results. This claim is supported by the time convergence of the profiles in Figure 9 toward the self-similar solutions. However, it is not clear to what extent letting the simulation run for longer, so that shells with ever larger will fall to the shock, helps here, since some of the inaccuracy of the simulation results is due to numerical dissipation. Additionally, as we explained in Section 3.4.3, the maximum runtime yielding physically admissible results is bounded, and simulations after a longer time will begin to depart from the self-similar solution generally.
We next study the effect of hybrid boundary condition scheme introduced in Section 3.4.3. Figure 10 gives the tail part of the density profiles of three simulations which are identical except for the implementations of the boundary conditions. The three panels correspond to the three types of boundary conditions mentioned above, respectively homogeneous Neumann (outflow), pure linear extrapolation and the hybrid scheme. A clear abnormal uprising of density near the boundary can be seen in the homogeneous case (the first panel), as it overestimates the inflow from beyond the boundary. This effect would "propagate" inwards and eventually pollute the solution, making it unstable. On the other hand, the density drops to under unity (or the critical density in physical units) when we use the extrapolated boundary condition (the second panel), leading to a negative density later in the simulation. By using the hybrid scheme (the third panel), An outward-propagating shock is clearly visible in all three panels. The curves are sampled over the positive direction of the axis, but we have checked that for all the simulations we report in this section the solution only has a very weak dependence on the direction along which we extract it from the simulation domain, see Appendix A. Figure 9. (Colour Online) The rescaled density, velocity and pressure profiles from the the same simulation of the Einstein-de Sitter model, plotted against the rescaled radial coordinate, . The self-similar theoretical prediction (Bertschinger 1985) is shown as black dashed lines. The vertical dashed lines with colours indicate the locations of the tophat edge at the same five times as shown in Fig. 8, and the numerical solutions depart from the self-similar prediction within it. This location is moving inwards as the rescaling radius ta increases over time. Convergence over time to the theoretical solution can be observed in the plots.
we manage to keep a relatively stable and smooth density evolution near the boundary throughout the simulation.

DGP models
In this subsection we report the simulation results of the DGP model introduced in Sections 2.1 and 2.2.2, with = 1.0 and various values of the screening parameter . For a clear comparison with the standard gravity, we first rescale our modified gravity results using the same turnaround radius formulation Eq. (15), following what we did first in Section 2.3. As the rescaling radius is identical in the different gravity models, the differences after the rescaling also represent the difference in the real evolution, thus showing the effects of modified gravity and the screening mechanism.
The results at ≈ 0.076 for models with = 1.0 and = 0, 1, 5, 10, 50, 100 are summarised in Figure 11. Those results agree with what one should expect for an enhanced gravity force and presence of screening: for the non-screening case ( = 0), in which gravity is constantly enhanced in time and space, a stronger shock is observed and it also happens at a larger radius than in EdS. In the other cases, as the screening becomes stronger and stronger (i.e., increasing ), the results approach that of standard gravity in an EdS universe.  One may have noticed that we require a bigger to achieve a similar screening effect, compared to Figure 1. This is mainly due to the fact that the parameter Δ, which characterises the mean initial overdensity density within some given initial radius , takes different values at the different initial radii covered by a real simulation, while the theoretical profiles are obtained assuming a fixed Δ, e.g., Δ = 0.001. To get rid of the Δ dependence in our results, we use the idea of rescaling using the true turnaround radius as described in Section 2.3. The difference is that this time we do not need a "re-rescaling": after we restore the profile quantities in physical units, we find the real turnaround radius directly by its physical meaning, i.e., we locate the radius where the physical velocity crosses zero. This method can be applied to all models including the EdS, which we have checked explicitly to give the same result as in the subsection above. After we located this real turnaround radius for simulations with DGP model, we use this value for our rescaling. The result of the same simulations and same timestamp in this new rescaling scheme then are plotted as the solid lines in Fig. 12, and their theoretical predictions (as shown in Fig. 3) are overplotted as the dashed lines with the same colour scheme. In the figure titles, we have used primes to indicate the quantities calculated using the numerically determined turnaround radius, ta .
Just like the theoretical results we got in 2.3, the new rescaled solutions are close to that in EdS universe. They are broadly in line with the theoretical predictions as well. The shock in DGP model happens at a slightly smaller radius, and the velocity in the gas shell just outside the shock has a bigger magnitude. This result is possibly caused by the fact that the gravitational force is stronger in the DGP  Fig. 11 (solid lines) for the DGP model with different parameters (see legends), now rescaled using the real turnaround radius ta as described in Section 2.3. The quantities with a prime are calculated using this new rescaling radius. We also plot the theoretical self-similar predictions described in Section 2.3 for each case, as dashed lines with the corresponding colours. model, so that the collapse is also stronger and faster. The qualitative trend is also as expected, as the curves for the models with screening are between the ones of EdS and a constant enhancement of Newton's constant ( = 0). Given that the real physical evolutions of these models in the simulation are quite different (cf. Fig. 11), these results demonstrate the reliability of E H PE 2 engine to carry out both standard and modified gravity simulations, and support the idea that self similarity can be found (at least as a very good approximation) in more general gravity models beyond EdS as well.

DISCUSSION AND CONCLUSION
To summarise, we have studied the spherical collapse of collisional gas in both an Einstein-de Sitter universe and DGP gravity model in this paper. We have derived self-similar solutions, for the first time, for some special cases of the latter class of models. The existence of self-similar solutions in spherical collapse scenarios is nontrivial: for example, while the EdS model admits a self-similar solution, this is lost if the model includes a cosmological constant. This is even more true for modified gravity models, in which the law of gravity may be modified in complicated time-and spatial-dependent ways. Indeed, we have tried to search for self-similar solutions in several classes of modified gravity theories that feature certain screening mechanisms. Chameleon-type models (Khoury & Weltman 2004a,b) do not admit self-similar solutions, because the fifth force there is not only scale dependent but also environment dependent. We have not found selfsimilar solutions for K-mouflage-type models (Babichev et al. 2009;Brax et al. 2013) either: in this model, the fifth force is given by where K is a parameter describing the coupling strength of the scalar field with matter, which is usually taken as a constant or function of time. The radial gradient d /d can be schematically obtained by solving where (·) is a nonlinear function, and ta is again the EdS expression of the turnaround radius, Eq. (15). For the fifth force to also respect self-similarity, it should be possible to express it as where F is a function of only. This is satisfied if (·) is a linear function and K is a constant, but this simply corresponds to a model with a constant enhancement of Newton's constant, identical to the DGP variant with = 0 considered above. For general (·), one has to require K to depend on both ta (and through which also depend on the initial radius and overdensity Δ) and to satisfy the above condition. Even for the DGP model we considered above, demanding a self-similar solution places some constraints on certain details, in particular the requirement that becomes a time-dependent function which grows at the same rate as the horizon size of the EdS universe. The existence of self-similar solutions in specific models offers us a way to test our numerical code for models other than EdS. The self-similar solutions we obtained for the modified DGP model behave as one would expect for an enhanced gravity with the Vainshtein screening mechanism at work. For example, we see that the shock happens at a larger radii in the DGP variant with = 0, and the infall velocity is larger outside the shock, compared with EdS, as a result of a stronger gravitational collapse. For the other DGP variants where > 0, the results generally lie between EdS and = 0, indicating a suppressed fifth force, and the suppression effect is larger for larger . It is notable that, despite the substantial differences in the evolutions and solutions of the different gravity models considered, after the (more 'proper') rescaling using the true turnaround radius of individual models, the solutions in the different DGP variants are all very close to that in the EdS model with standard gravity (though their agreement is not perfect). We also notice that, after this proper rescaling, the self-similar solutions in the DGP models depend very weakly on Δ, as also happens in EdS. Apparently, we should test these observations for other types of gravity models too. if they hold there as well, this is an interesting indication that the properly rescaled solutions in different gravity models are close to each other, which in turn implies that self-similarity should hold approximately, even though not exactly, in generic models. We leave a more detailed exploration of this possibility to future work.
Behind our new physical insights is a new implementation of cosmological hydrodynamical simulations of the spherical collapse scenario for different gravity models, based on the publicly-available hyperbolic PDE engine E H PE 2. We have described various technical details in our implementation, including the initial and boundary conditions which must be properly set up in order to get stable and correct evolutions. We find that the numerical simulations of the same EdS and DGP models as introduced above yield good agreements with the theoretical predictions we derived. This thus not only supports our findings on the self-similarity in the considered models, but also serves as a validation of the reliability and correctness of our E H PE 2 implementation.
By comparing our theoretical predictions to the simulation results, we find that, although to a large degree the code is capable of handling the collapse scenarios in different gravity models, there are still some inaccuracies in the current simulation results, in particular at and around the shock. The observed shift and weakening of the shock are likely caused by numerical dissipation, which may be suppressed by increasing the spatial and temporal resolutions. There are several possible ways of doing this. First, we are currently using a simple Finite Volume formalism which employs a generic Riemann solver. This scheme can be further extended to higher-order formalisms, e.g., Discontinuous Galerkin methods in combination with Runge-Kutta schemes or ADER-DG (Zanotti et al. 2015). Those schemes are in principle compatible with our scenarios, straightforward to implement, and could work properly to enhance the resolutions, but it remains an open question if these methods are well-suited to capture the steep gradients near the shock. We could also directly increase the resolution of our simulations, but this comes at additional runtime cost. Future work will look at the outsourcing of the individual patches to GPUs. This will provide us with the opportunity to work with significantly finer resolutions and a much higher efficiency. For the temporal side, we need to check if local time stepping or subcycling help to reduce the vulnerability of the current explicit time stepping scheme to numerical dissipation. It may also help to use more accurate Riemann solvers, as the current Rusanov solver only 'reacts' to the biggest eigenvalue of the system, cf. Eq. (72), so that it does not preserve the characteristics of all five evolving quantities well if they propagate with different wave speeds.
With a working hydrodynamical simulation code at hand, where new models of gravity can be straightforwardly implemented, a natural next step is to run simulations for more realistic modified gravity models that do not have self-similar solutions, including the original DGP model, the K-mouflage model and the chameleon model. For the latter we may need to either add a multigrid solver for the scalar field, or adopt some approximate solutions such as the thinshell solution. In a future project, we will compare the collapse of collisional gas in these different models in detail. If the above speculation, namely the spherical solutions rescaled by the true turnaround radii of models are approximately the same in different cosmological models, turns out to be correct, then the differences in the physical solutions of these models can be largely ascribed to the differences in their turnaround radii, which might offer a simple way to model the modified gravity effects. In addition, we plan to add more physical processes, such as radiative cooling (e.g., Abadi et al. 2000), in the code, to understand how they interfere with the effects of a modified law of gravity. Altogether, these will hopefully offer us new insights into the behaviour of gas, and hence the galaxy formation process, in modified gravity models.