Cosmological simulations with disformally coupled symmetron fields

We investigate statistical properties of the distribution of matter at redshift zero in disformal gravity by using N-body simulations. The disformal model studied here consists of a conformally coupled symmetron field with an additional exponential disformal term. We conduct cosmological simulations to discover the impact of the new disformal terms in the matter power spectrum, halo mass function, and radial profile of the scalar field. We calculated the disformal geodesic equation and the equation of motion for the scalar field. We then implemented these equations into the N-body code ISIS, which is a modified gravity version of the code RAMSES. The presence of a conformal symmetron field increases both the power spectrum and mass function compared to standard gravity on small scales. Our main finding is that the newly added disformal terms tend to counteract these effects and can make the evolution slightly closer to standard gravity. We finally show that the disformal terms give rise to oscillations of the scalar field in the centre of the dark matter haloes.


Introduction
It has been known since 1998 that the universe expands at an accelerating rate which is consistent with the existence of a cosmological constant Λ (Riess & et al. 1998;Perlmutter & et al. 1999). The standard model for cosmology, ΛCDM, gives an excellent fit to most of the modern precision observations on large scales and of the cosmic microwave background, but it does not explain what the source of Λ is. Attempts to calculate the vacuum energy from particle physics yields answers that are several orders of magnitude off from the measured cosmological value of Λ. A cancellation of that many terms is very improbable and would require an extreme fine tuning. This is the cosmological constant problem, and is a severe problem in modern physics. See Weinberg (1989) for an early introduction to the problem.
A viable solution to this problem might be that the particle physics vacuum energy is totally concealed on gravitational scales and other mechanisms are responsible for the measured expansion. One way to search for such mechanisms consist in introducing a slight modification to standard general relativity (GR) in such a way that the equations for gravity will give rise to accelerated expansion on large scales. There are innumerable models for modified gravity, see Clifton et al. (2012) for a review.
An important requirement that any modification of gravity must fulfill to be considered viable is that it should reduce to GR on solar system scales. This is achieved through so called screening mechanisms, described in detail by Joyce et al. (2014). Screening mechanisms are needed because conventional GR is tested and confirmed to very E-mail: robert.hagala@astro.uio.no high precision on solar system scales, therefore any modifications to gravitational physics must give similar results within very tight constraints on these scales. See the review by Reynaud & Jaekel (2009) for an overview over solar system tests and constraints.
In this paper we will investigate through N -body simulations a specific form for modified gravity, namely a disformal model of gravity. Disformal models were first introduced by Bekenstein (1993), and now have been widely studied with applications to inflation, dark energy and dark matter (Kaloper 2004;Koivisto 2008;Zumalacárregui et al. 2010;Koivisto et al. 2012;van de Bruck et al. 2013;Bettoni & Liberati 2013;Neveu et al. 2014;Brax & Burrage 2014;Deruelle & Rua 2014;Tsujikawa 2014;Sakstein 2014Sakstein , 2015Koivisto & Urban 2015). However, this is the first time disformal models have been studied on nonlinear scales.
We consider a model where the scalar field has both conformal and disformal couplings to matter. For the conformal part, we use a symmetron potential (Pietroni 2005;Olive & Pospelov 2008;Hinterbichler & Khoury 2010). While for the disformal part we use an exponential term.
Conformally invariant Symmetron models have been already investigated in the non-linear regime (Davis et al. 2012;) as well as several other conformally invariant scalar tensor theories (Boehmer & Mota 2008;Li et al. 2008;Baldi et al. 2010;Zhao et al. 2011;Li et al. 2011;Li et al. 2013;Hellwing et al. 2014). The aim of this paper is to investigate for the first time the effects of adding a disformal term to the symmetron field. We will focus our analysis mainly in the statistical properties (i.e. the power spectrum and halo mass function) of the simulated matter distribution on non-linear scales. The simulations are done with a modified version of the N -body code Isis , which is itself a modified gravity version of the code Ramses (Teyssier 2002).
The paper is organised as follows: Section 2 describes the equation of motion of the scalar field and the associated geodesic equation; it also summarises how these are implemented into Isis and tested. Section 3 describes the cosmological simulations that are used for the analysis and section 4 contain the results from the analysis. We conclude and discuss in section 5.

The model
The model is defined by the following scalar-tensor action where g andḡ are the Einstein and Jordan frame metrics, R is the Ricci scalar, andL m is the Lagrangian density of matter (computed using the Jordan frame metricḡ whenever applicable). The field potential V (φ) can have many different forms, but we choose the symmetron Mexican hat potential with three free parameters µ, λ and V 0 , In disformal gravity the Jordan frame metricḡ is related to the Einstein frame metric according tō The specific forms of A and B we will study in this paper are as follows: where B 0 and β are free parameters for the disformal coupling. The normalization constant φ 0 is chosen to be the vacuum expectation value of the symmetron field φ 0 ≡ µ √ λ . The mass scale M is a free parameter, deciding the interaction strength of the conformal coupling.
This specific choice of A, B and V gives a symmetron model with an additional disformal term described by B (φ). The model reduces to the symmetron model when fixing B 0 = 0. Both the matter coupled symmetron part A and the disformal part B can lead to two separate screening effects (Hinterbichler & Khoury 2010;Koivisto et al. 2012). Note that we use natural units where c = = 1.

The equation of motion for the scalar field
The equation of motion for the scalar field that results from varying the action defined in eq.1 is: where we have defined See Koivisto et al. (2012) for details on the derivation. Here, the dots are partial derivative with respect to cosmic time. The Einstein frame metric is assumed to be a flat FRW metric with a scalar perturbation Ψ: a is the expansion factor, H =ȧ a is the Hubble parameter, and ρ is the total matter density. Note that we have not neglected the non-static terms in the scalar field since in this particular model there is a coupling between the density and the time derivatives, whose effects are still not fully understood in the non-linear regime of cosmological evolution.
For convenience, we normalize the field to the vacuum expectation value of the symmetron field, Such that the new dimensionless field χ = φ/φ 0 should keep in the range χ ∈ [−1, 1], at least for a symmetrondominated case when B 0 is small. Also for numerical convenience, we introduce the parameter a SSB , defining the expansion factor at spontaneous symmetry breaking assuming the conformal field and a uniform matter distribution. We further define the density at which the symmetry is broken a dimensionless coupling constant and the range of the field in vacuum By taking into account these definitions, we can re-write the equation of motion for the scalar field (eq. 6) as where we have also fixed the three free functions as stated above and defined What remains before we can insert this equation into Isis, is a switch to supercomoving time -the time variable used by Isis/Ramses, defined by Martel & Shapiro (1998) -and to split this second order differential equation into a set of two coupled first order differential equations. Supercomoving time τ is related to the cosmic time t by dτ = 1 a 2 dt. Finally, we introduce the variable which leads to the following set of first order differential equations for q and χ: where the primes denote derivatives with respect to supercomoving time.
The supercomoving variables with tildes are defined as H ≡ a 2 H andΨ ≡ a 2 Ψ. Setting γ 2 ρ = B = 0 and A ≈ 1 here will result in equations equivalent with equations (22) and (23) of .
These equations are solved using the leap-frog algorithm. The time scales associated to the oscillations of the field are much smaller that those associated to the movement of matter. Because of this, we use shorter time leapfrog steps for the scalar field than for matter in the same way as Llinares & Mota (2014) did.

The geodesic equation
In this theory of gravity, dark matter particles move in geodesics determined by the Jordan frame metric. Generally the geodesic equation reads whereΓ µ αβ are the Christoffel symbols that correspond to the Jordan frame metric. Assuming nonrelativistic particles we neglect quadratic terms in the velocity. The geodesic equation then takes the form To find the necessary barred Christoffel symbols, we use the equation from Zumalacárregui et al. (2013), The resulting equation of motion for the dark matter particles is given bÿ Here it is possible to recognise the acceleration terms associated with perturbed FRW geodesics in standard gravity, namelÿ where the second term is the standard gravity force and the third term is the Hubble friction. These terms are already included in Ramses, so we will focus on the other terms.
The effect of modified gravity on the trajectory of the particles is through the addition of a fifth force as well as new damping terms proportional toẋ. As a first attempt to solve this equations, we neglect the damping terms and focus our analysis on the impact of the extra terms associated to the fifth force on a slowly moving matter distribution. The only friction term that we keep is the one associated with the expansion of the universe, which is automatically taken into account in the solution when writing the equation as a first order equation. The acceleration due to the fifth force using supercomoving time and with the dimensionless field χ that we inserted into the code is given by

Tests
In order to test that our modifications to the code Isis are properly implemented, we compare results that are obtained with our code in the pure conformal limit (i.e. when we set B 0 = 0) with those that were presented by Llinares & Mota (2014) using a pure disformal symmetron non-static code. We also test the dependence of the final results on the initial conditions that are used for the non-static solver.

Cosmological comparison to symmetron results
We run two simulations using both the code presented in this paper and the one described in   model we used B 0 = 0, such that the disformal code should be able to reproduce the already published symmetron nonstatic results. Both simulations are done with the same initial conditions, which are generated with the code Grafic (Bertschinger 2001) using a box size of 64 Mpc/h and 128 3 particles. The comparison of the simulations is focused on the power spectrum of over-density perturbations and the field profiles of the most massive halo. Fig.1 shows a comparison of the field profile for the most massive halo found in the simulation (with a mass of 2.5 × 10 14 M /h). Because the field oscillates, we compare mean values in time that are calculated by averaging over the time interval from a = 0.995 to a = 1. The differences in the field profiles from the two codes are always below 0.5% and is sufficiently small to be attributed to the initial random values of the field. A comparison between the power spectra at redshift z = 0 from both simulations is presented in fig.2. The relative difference between the power spectra is well below 0.5% in the whole domain of simulated frequencies.

Dependence on initial conditions for the scalar field
The initial conditions for the scalar field are generated in the same way as in : the initial values for χ are drawn from a uniform distribution around zero with a small dispersion (χ 0 ∈ [−ε, ε]). We test the sensitivity of the redshift zero scalar field profiles against changes in the amplitude of the initial conditions. To this end, we run simulations with the only aim of tracking the evolution of the scalar field. As we are only interested in the field profiles, we use standard gravity in the geodesic equation to ensure identical matter distributions. The box size and number of particles is 64 Mpc/h and 64 3 respectively. We do several runs changing the amplitude of the initial conditions with values ranging from 10 −5 to 10 −17 . Fig.3 shows the time averaged scalar field profile that corresponds to the most massive object found at redshift z = 0 for all these simulations, which has a mass of 1.2 × 10 14 M /h. We find no significant deviations when varying ε. For the scientific simulations that will be presented in the following section, we choose to use ε = 10 −13 .

Cosmological simulations
In order to quantify the effects of the disformal terms in the cosmological evolution, we run a set of cosmological simulations using Newtonian gravity, a pure symmetron model and the symmetron model plus the disformal terms. The initial conditions are generated with the package Grafic (Bertschinger 2001) with standard gravity. The approximation that we make when not including modified gravity in the initial conditions is fully justified by the fact that because of screening effects, modified gravity starts to act only after the redshift of symmetry breaking has passed, which we choose to be at a much later time. All simulations use exactly the same initial conditions and assume a flat ΛCDM background cosmology provided by the Planck simulation b 0 β GR (no modifications) --Symmetron-like 0 0 Disformal A 1 1 Disformal B 2 2 Disformal C 1 0  The initial values of the unitless scalar field χ at the initial time is calculated assuming the scalar field is fully screened at the initial time of the simulation, and thus taken from a uniform random distribution with χ 0 ∈ [−10 −13 , 10 −13 ]. The initial time derivative of the scalar field is assumed to be zero. As discussed in previous section, we found that the evolution of the scalar field is not very sensitive to the assumed initial conditions. Table 1 describes the model parameters employed in the simulations. The simulation Disformal A has what we consider standard parameters. Disformal B has an amplified disformal part due to increased b 0 and β. Disformal C has β = 0, which sets B (φ) to be constant. A constant disformal term might give some insight about the effect of the disformal term alone, since the equations show that the positive and negative parts of the symmetron effective potential have different shapes when β = 0. The Symmetron simulation is run with the code presented in this paper, but with the disformal part of the equations turned off by setting b 0 = 0. In order to be able to compare with previous results, we use the same symmetron model parameters that were assumed in : (λ 0 , a SSB , θ) = (1 Mpc/h, 0.5, 1).

Power spectrum and halo mass function
We study the impact of disformal terms in the power spectrum of overdensity perturbations. The estimation of the power spectrum is made using a Fourier based method. Dis- creteness effects are corrected using the method proposed by Jing (2005). Figure 4 shows the relative difference between the modified gravity and GR power spectra.
We find that the pure conformal model has the same effects found by Davis et al. (2012) and : the fifth force has no effect at large scales (which makes the model consistent with the observed normalization of the perturbations) and gives an excess of power at small scales. Regarding the disformal terms, we found that in the simulations Disformal A and C, the disformal part of the model has no significant impact in the modifications made to GR by the symmetron model. In the case of Disformal B, the stronger disformal terms counteract the effects of the symmetron, bringing the power spectrum closer to GR.
Additionally, we study the halo mass function. In order to extract this quantity from the simulation data, we identify the dark matter haloes with the halo finder Rockstar (Behroozi et al. 2013), which uses a 6D friends-of-friends algorithm. The cumulative mass functions for all the simulations are presented in figure 5. The behaviour is the same as found in the power spectrum: the symmetron model increases the number of small haloes, while the strongest disformal terms act in the opposite direction, diminishing the symmetron effect.

Field profiles
We present the radial profile of the scalar field for the most massive halo found in the simulation at redshift z = 0, which has a mass of M = 5.1 × 10 14 M /h. In the case of static simulations, it is enough to extract information from the last snapshot that is output by the simulation code. However, in the case of non-static simulations like those presented in this paper, the scalar field has oscillations which makes it difficult to make a comparison between different simulations (the scalar field will often be in a different phase at a given time). To overcome this problem, we calculate a mean value in time taking into account several oscillation periods. The mean is calculated over the interval from a = 0.995 to a = 1. We assume here that the variations produced by the displacement of matter during that interval are much smaller than those related to the oscillations of the field.
The upper panel of figure 6 shows the mean profile of the scalar field. In the simulation Disformal A, the scalar field chose to oscillate around the negative minimum (the symmetron potential has two different minimum with different sign), in that case we show − χ . The fact that the potential has two minima should lead to the formation of domain walls (Llinares & Mota 2013;Llinares & Pogosian 2014;Pearson 2014), which we do not find in our simulations owing to the small size of the box. Domain walls might have formed in our box at an earlier stage, however, and collapsed before redshift zero.
The scalar field profile tends to the vacuum value far from the halo and presents a gradient responsible for a fifth force when approaching the halo. The innermost region of the cluster is screened, and thus its corresponding field value tends to zero in that region. The profile shows almost no variation among the different models. However, we find significant differences between standard and disformal symmetron models in the dispersion of the field profiles, which is a measurement of how large the amplitudes of the scalar field oscillations are. This information is shown in the lower panel of figure 6. In the symmetron case, we find that the dispersion of the profile goes to zero in the centre of the cluster, which means that the scalar field is completely still there and only has oscillations in the outer regions and the voids surrounding the halo. In the disformal case, the situation is different. We find that the amplitude of the oscillations does not go to zero in the centre of the cluster and that, in fact, the amplitude can even increase in the central region with respect to the values found outside of the halo. We also see that the dispersion is slightly smaller than in the symmetron case in the low density region (r > 1.5r vir ). The reason for this may be related to the fact that the speed of sound of the scalar field is smaller in high density regions (see eq.6, which has the form of a wave equation where the factor of (1 + γ 2 ρ) in front ofφ can be regarded as 1/c 2 s ). Thus, inside massive haloes the waves of the scalar field can cluster, giving rise to extra oscillations. Further analysis must be made to confirm this hypothesis.

Conclusions
We run for the first time cosmological simulations with disformally coupled scalar fields. The model includes both a conformal and a disformal coupling to matter. For the conformal part, we choose a symmetron potential and conformal factor A = 1 + (φ/M ) 2 . The disformal part is given by an exponential factor B = B 0 exp (βφ/φ 0 ). The aim of the paper is to test the effects of the disformal factor on the already known results for the symmetron model.
We present the formalism and modifications that we made to the code Isis  to be able to simulate disformal gravity. The disformal code presented here is based on the non-static solver for scalar fields that is part of the Isis code . We present results obtained from a set of five cosmological simulations.
In the set of simulations, we find that the disformal terms counteract some of the clustering effects of the symmetron field: the power spectrum of density perturbations and the halo mass functions are both smaller than in the  Fig. 6. The upper panel shows the absolute value of the radial field distribution for the most massive halo found in the simulation for the symmetron and disformal simulations. The lower panel shows the dispersion σχ obtained over several oscillations of the scalar field, which is a good estimation of the amplitude of the oscillations. The field is time smoothed from a = 0.995 to a = 1. symmetron case. However, at least in the region of the parameter space studied, the conformal part of the Lagrangian is dominant. Therefore, the end result is nevertheless a small scale increase of both the power spectrum and the mass function with respect to GR.
In order to understand the differences in the distribution of the scalar field that appear because of the new terms, we study the field profiles that correspond to the most massive halo found in the simulations. We find almost no differences in the field profile, but we do find differences in the amplitude of the oscillations of the field, being larger for the disformaly coupled models. This implies that the differences in the power spectra and mass function when comparing the disformal models and the symmetron are not due to the field value, but the field oscillations. While the gradient of the scalar field is independent of the presence of the disformal terms, it is important to keep in mind that the expression for the fifth force includes the time derivatives of the scalar field. Hence, the fifth force can be modified with respect to the purely conformal fifth force only by increasing the derivatives.