EDP Sciences
Open Access
Issue
A&A
Volume 585, January 2016
Article Number A37
Number of page(s) 7
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/201526439
Published online 14 December 2015

© ESO, 2015

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Since 1998, it has been known that the universe expands at an accelerating rate that is consistent with the existence of a cosmological constant Λ (Riess et al. 1998; Perlmutter et al. 1999). The standard model for cosmology, Λ cold dark matter (ΛCDM), gives an excellent fit to most modern precision observations of large scale structures 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 extreme fine tuning. This is known as the cosmological constant problem and is a severe problem in modern physics (see Weinberg 1989, for an early introduction to this issue).

A viable solution to the cosmological constant problem is to assume that the particle physics vacuum energy is totally concealed on gravitational scales, while other mechanisms are responsible for the measured expansion. One way to search for such mechanisms consist of 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).

To be considered viable, an important requirement that any modification of gravity must fulfil 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. (2015). Screening mechanisms are needed because conventional GR is tested and confirmed to very high precision on solar system scales, meaning 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 investigate a specific form for modified gravity by using N-body simulations, namely the disformal model of gravity. Disformal models were first introduced by Bekenstein (1993), and have now been widely studied by applying them 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; Brax et al. 2013a; Bettoni & Liberati 2013; Neveu et al. 2014; Brax & Burrage 2014; Deruelle & Rua 2014; Tsujikawa 2015; Sakstein 2014, 2015; Koivisto & Urban 2015; van de Bruck & Morrice 2015). However, this is the first time disformal models have been studied on non-linear 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 already been investigated in the non-linear regime (Davis et al. 2012; Llinares & Mota 2014; Llinares et al. 2014), as well as several other conformally invariant scalar tensor theories (Boehmer & Mota 2008; Li et al. 2008, 2011, 2013; Baldi et al. 2010; Zhao et al. 2011; Brax et al. 2013b; 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 on the statistical properties (i.e. the power spectrum and halo mass function) of the simulated matter distribution on non-linear scales. The simulations are performed with a modified version of the N-body code Isis (Llinares et al. 2014), which is itself a modified gravity version of the code Ramses (Teyssier 2002).

The paper is organised as follows: Sect. 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 Sect. 4 contains the results from the analysis. We discuss the results and draw conclusions in Sect. 5.

2. The equations and the code

2.1. The model

The model is defined by the following scalar-tensor action (1)where g and are the Einstein and Jordan frame metrics, R is the Ricci scalar, and 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 quartic symmetron potential with the three free parameters μ, λ, and V0, (2)In disformal gravity the Jordan frame metric is related to the Einstein frame metric according to (3)The specific forms of A and B that we will study in this paper are as follows: where B0 and β are free parameters for the disformal coupling. The normalization constant φ0 is chosen to be the vacuum expectation value of the symmetron field . The mass scale M is a free parameter, which decides 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 B0 = 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). We use natural units where c = ħ = 1.

2.2. 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: (6)where we define (7)See Koivisto et al. (2012) for details on the derivation. Here, a dot represents a partial derivative with respect to cosmic time. The Einstein frame metric is assumed to be a flat Friedmann–Lemaître–Robertson–Walker metric with a scalar perturbation Ψ, specifically (8)The symbol a is the expansion factor, 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 matter 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, (9)As such, the new dimensionless field χ = φ/φ0 should stay in the range χ ∈ [−1,1], at least for a symmetron-dominated case when B0 is small. Also, for numerical convenience, we introduce the parameter aSSB, which defines the expansion factor at the time of spontaneous symmetry breaking, assuming a uniform matter distribution and no disformal term. We further define the density at which the symmetry is broken: (10)a dimensionless symmetron coupling constant (11)the range of the symmetron field in vacuum (12)and a dimensionless disformal coupling constant (13)By taking into account these definitions, we can rewrite Eq. (6)as (14)where we have also fixed the three free functions as stated above and defined (15)What remains to be done, before we can insert this equation into Isis, is to 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 . Finally, we introduce the variable (16)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 and . Setting γ2ρ = B = 0 and A ≈ 1 here will result in symmetron equations equivalent to Eqs. (22) and (23) of Llinares & Mota (2014).

These equations are solved using the leapfrog algorithm. The time scales associated with the oscillations of the field are much smaller that those associated with the movement of matter. Because of this, we use shorter leapfrog time steps for the scalar field than for matter, as described in Llinares & Mota (2014).

2.3. 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 (19)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 (20)To find the necessary barred Christoffel symbols, we use the equation from Zumalacárregui et al. (2013), (21)The resulting equation of motion for the dark matter particles is given by (22)Here it is possible to recognise the acceleration terms associated with perturbed FLRW geodesics in standard gravity, namely (23)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 that are proportional to . As a first attempt to simulate this model, we neglect the damping terms and focus our analysis on the impact of the extra terms associated with 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 on supercomoving form. The extra acceleration that results from the fifth force, using supercomoving time and with the dimensionless field χ, that we inserted into the code is given by (24)

2.4. Tests

To test that our modifications to the code Isis were properly implemented, we compare results that were obtained with our code in the pure conformal limit (i.e. when we set B0 = 0) to 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.

thumbnail Fig. 1

Comparison of the average field profiles around a massive halo. The top panel shows the field profile for the symmetron code from Llinares & Mota (2014; dashed red) and for our disformal code, in the pure conformal limit (solid black). The bottom panel shows the relative difference between the different codes.

Open with DEXTER

thumbnail Fig. 2

Relative difference between the matter power spectra from our disformal code, in the pure conformal limit, and the symmetron code from Llinares & Mota (2014).

Open with DEXTER

2.4.1. Cosmological comparison to symmetron results

We run two simulations using both the code presented in this paper and the one described in Llinares & Mota (2014) with the following symmetron parameters: θ = 1, λ0 = 1 Mpc, and aSSB = 0.5. For the disformal part of the model, we used B0 = 0, such that the disformal code should be able to reproduce the already published symmetron non-static 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 1283 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.

Figure 1 shows a comparison of the field profile for the most massive halo found in the simulation (with a mass of 2.5 × 1014M/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 are 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.

thumbnail Fig. 3

Sensitivity of the field profiles to differences in the initial conditions for the scalar field. The curves correspond to the average field profile around a massive halo. Different curves correspond to different amplitudes ϵ of the initial random field perturbations.

Open with DEXTER

2.4.2. Dependence on initial conditions for the scalar field

The initial conditions for the scalar field are generated in the same way as in Llinares & Mota (2014): 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 are 64 Mpc /h and 643 respectively. We do several runs changing the amplitude of the initial conditions with values ranging from 10-5 to 10-17. Figure 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 × 1014M/h. We find no significant deviations when varying ε. For the scientific simulations that we present in the following section, we chose to use ε = 10-13.

3. Cosmological simulations

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 matter distribution is 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 matter distribution and assume a flat ΛCDM background cosmology provided by the Planck collaboration: Ωm = 0.3175, ΩΛ = 0.6825, and H0 = 67.11 km s-1 Mpc-1 (Planck Collaboration XVI 2014). The number of particles is 2563, and the size of the box is 64 Mpc /h.

The initial values of the dimensionless scalar field χ at the initial time is calculated by assuming that 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 the previous section, we found that the evolution of the scalar field is not very sensitive to the assumed initial conditions.

Table 1

Model parameters for the different simulation runs.

Table 1 describes the model parameters employed in the simulations. The model Disformal A has what we consider to be standard parameters. Disformal B has an amplified disformal part owing to increased b0 and β. Disformal C has β = 0, which sets B(φ) to be constant. A constant disformal term might give some insight into the effect of the disformal term alone, since the equations show that the positive and negative parts of the disformal effective potential have different shapes when β ≠ 0. The Symmetron-like simulation is run with the code presented in this paper, but with the disformal part of the equations turned off by setting b0 = 0. In all of the modified gravity simulations, we use the symmetron model parameters (λ0,aSSB) = (1 Mpc /h,0.5,1).

After doing one simulation for each set of parameters, as described in Table 1, we found that Disformal B gave the most unexpected results – as will be presented in the next section. We therefore performed a total of four simulations of Disformal B, with different initial seeds for the scalar field, but with identical parameters and initial matter distribution. In two of these simulations, the scalar field ended up in the positive minimum of the effective potential (φ ≈ + φ0), while in two others, the scalar field ended up in the negative minimum (φ ≈ −φ0).

4. Results

4.1. Power spectrum and halo mass function

thumbnail Fig. 4

Relative difference of the matter power spectrum of the four models of modified gravity with respect to GR.

Open with DEXTER

We study the impact of disformal terms in the power spectrum of density perturbations. The estimation of the power spectrum is made using a Fourier based method. Discreteness 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 that were found by Davis et al. (2012) and Llinares & Mota (2014): 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 on the modifications made to GR by the symmetron model. In the case of the first simulation of Disformal B, the stronger disformal terms counteract the effects of the symmetron, bringing the power spectrum closer to GR.

Figure 5 shows the relative difference between all four simulations of Disformal B, and the GR power spectrum. The two disformal simulations where the field fell to the positive potential minimum, have the strongest suppression of power on all scales. For the two other disformal simulations, where the field stabilized around the negative potential minimum, there is some reduction on small scales (k> 3 h Mpc-1), but actually an increase on large scales (k ≈ 1 h Mpc-1).

thumbnail Fig. 5

Relative difference of the matter power spectrum of four different simulations of Disformal B with respect to standard gravity (GR). The legend indicates the sign of the potential minimum to which the scalar field fell.

Open with DEXTER

thumbnail Fig. 6

Cumulative mass function for each of the five simulations (top panel). The relative difference of the four modified models compared to GR is shown in the bottom panel.

Open with DEXTER

Additionally, we study the halo mass function. 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 Fig. 6. The behaviour is the same as that found in the power spectrum: the symmetron model increases the number of small haloes, while the strong disformal term in Disformal B acts in the opposite direction, diminishing the symmetron effect.

These findings on non-linear scales agree with earlier results from linear evolution in disformal theory. In the paper by van de Bruck & Morrice (2015), they find on the linear perturbation level, that disformal terms counteract the conformal coupling to some degree.

4.2. 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 × 1014M/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 top panel of Fig. 7 shows the mean profile of the scalar field. In the Symmetron-like, Disformal A, and Disformal C simulations, the scalar field chose to oscillate around the negative minimum (the symmetron potential has two different minima with different signs). In these cases we show −⟨ χ for a better comparison of the shape of the field profile. 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 towards 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 bottom panel of Fig. 7. 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 at rest 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 less than in the symmetron case in the low density region (r> 1.5rvir). The reason for this may be related to the speed of waves of the scalar field being less 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 ). Thus, inside massive haloes, the waves of the scalar field can cluster and give rise to extra oscillations. Further analysis must be made to confirm this hypothesis.

thumbnail Fig. 7

Radial field distribution for the most massive halo found in the simulations. The top panel shows the average field value at a distance r from the centre of the halo. The bottom 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.

Open with DEXTER

5. Conclusions

For the first time, we run 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 a conformal factor A = 1 + (φ/M)2. The disformal part is given by an exponential factor B = B0exp(βφ/φ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 (Llinares et al. 2014) 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 (Llinares & Mota 2014). We present results obtained from a set of five cosmological simulations.

In the set of simulations, we find that the stronger disformal terms can 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 symmetron case. However, at least in the region of the parameter space studied, the conformal part of the model 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 the first simulation we found that Disformal B was the only model where the field fell to the positive minimum of the symmetron potential. This could be very important, since the disformal factor B(χ) ∝ b0exp(βχ) will change by several orders of magnitude when χ → −χ. The reason for the field falling to one minimum or the other in a particular simulation is complex and can be attributed to the chaotic nature of the equation of motion for the field. We did a total of four simulations on Disformal B with varying initial field distributions. The results show that the power spectrum was reduced significantly only for positive field values. This indicates that the exponential shape of B(φ) is more important than the value of b0 to achieve the masking of the symmetron clustering. Furthermore, we will probably see new physics and observational signatures of this model in future studies of domain walls (see Correia et al. 2014, for a description of domain walls in asymmetric potentials).

To understand the differences in the distribution of dark matter 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, which are larger for the disformaly coupled models. This implies that the differences in the power spectra and mass functions when comparing the disformal models and the symmetron are not the result of the field value, but of 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 by only increasing the derivatives.

Possible observational consequences of the field oscillations in haloes are (time varying) changes in how photons are lensed by the oscillating disformally coupled field. Because light rays will follow geodesics that are dictated by the modified Jordan frame metric, anomalies in cluster lensing masses could indeed be a “disformal smoking gun”. This possibility will be studied in a future paper.

Acknowledgments

The simulations were performed on the NOTUR cluster HEXAGON, the computing facilities at the University of Bergen. C.L.L. and D.F.M. acknowledge support from the Research Council of Norway through grant 216756.

References

All Tables

Table 1

Model parameters for the different simulation runs.

All Figures

thumbnail Fig. 1

Comparison of the average field profiles around a massive halo. The top panel shows the field profile for the symmetron code from Llinares & Mota (2014; dashed red) and for our disformal code, in the pure conformal limit (solid black). The bottom panel shows the relative difference between the different codes.

Open with DEXTER
In the text
thumbnail Fig. 2

Relative difference between the matter power spectra from our disformal code, in the pure conformal limit, and the symmetron code from Llinares & Mota (2014).

Open with DEXTER
In the text
thumbnail Fig. 3

Sensitivity of the field profiles to differences in the initial conditions for the scalar field. The curves correspond to the average field profile around a massive halo. Different curves correspond to different amplitudes ϵ of the initial random field perturbations.

Open with DEXTER
In the text
thumbnail Fig. 4

Relative difference of the matter power spectrum of the four models of modified gravity with respect to GR.

Open with DEXTER
In the text
thumbnail Fig. 5

Relative difference of the matter power spectrum of four different simulations of Disformal B with respect to standard gravity (GR). The legend indicates the sign of the potential minimum to which the scalar field fell.

Open with DEXTER
In the text
thumbnail Fig. 6

Cumulative mass function for each of the five simulations (top panel). The relative difference of the four modified models compared to GR is shown in the bottom panel.

Open with DEXTER
In the text
thumbnail Fig. 7

Radial field distribution for the most massive halo found in the simulations. The top panel shows the average field value at a distance r from the centre of the halo. The bottom 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.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.