N-Body Simulations for Extended Quintessence Models

# N-Body Simulations for Extended Quintessence Models

Baojiu Li, David F. Mota and John D. Barrow DAMTP, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge CB3 0HA, UK Institute of Theoretical Astrophysics, University of Oslo, 0315 Oslo, Norway
###### Abstract

We introduce the -body simulation technique to follow structure formation in linear and nonlinear regimes for the extended quintessence models (scalar-tensor theories in which the scalar field has a self-interaction potential and behaves as dark energy), and apply it to a class of models specified by an inverse power-law potential and a non-minimal coupling. Our full solution of the scalar field perturbation confirms that, when the potential is not too nonlinear, the effects of the scalar field could be accurately approximated as a modification of background expansion rate plus a rescaling of the effective gravitational constant relevant for structure growth. For the models we consider, these have opposite effects, leading to a weak net effect in the linear perturbation regime. However, on the nonlinear scales the modified expansion rate dominates and could produce interesting signatures in the matter power spectrum and mass function, which might be used to improve the constraints on the models from cosmological data. We show that the density profiles of the dark matter halos are well described by the Navarro-Frenk-White formula, although the scalar field could change the concentration. We also derive an analytic formula for the scalar field perturbation inside halos assuming NFW density profile and sphericity, which agrees well with numerical results if the parameter is appropriately tuned. The results suggest that for the models considered, the spatial variation of the scalar field (and thus the locally measured gravitational constant) is very weak, and so local experiments could see the background variation of gravitational constant.

## 1. Introduction

The nature of the dark energy (Copeland, Sami & Tsujikawa, 2006) is one of the most difficult challenges facing physicists and cosmologists now. Although a cosmological constant (plus cold dark matter, to provide the concordance CDM paradigm) could be a solution – and is indeed consistent with virtually all current observations, it suffers from theoretical difficulties such as why its value must be so small yet nonzero, and why it becomes dominant only at the low redshift. In all the alternative proposals to tackle this problem, a quintessence scalar field (Zlatev, Wang & Steinhardt, 1999; Wang et al., 2000) is perhaps the most popular one (although a new proposal by Barrow & Shaw (2010) provides a completely new type of explanation that does not require new scalar fields). In such models the scalar field is slowly rolling down its potential, its energy density is dominated by the potential energy and almost remaining constant provided that the potential is flat enough. The flatness of the potential, however, means that the mass of the scalar field is in general very light and as a result the scalar field almost does not cluster so that its effects in cosmology are mainly on the (modified) background expansion rate.

One reason for the wide interest in quintessence models is that scalar fields appear in abundance in high-energy physics theories, in which they are often coupled to the curvature invariants or even other matter species, leading to the so-called extended quintessence (Perrotta, Baccigalupi & Matarrese, 2000; Baccigalupi & Perrotta, 2000; Baccigalupi, Perrotta & Matarrese, 2000) and coupled quintessence (Amendola, 2000, 2004; Jesus et al., 2008) models respectively. The former is just a special class of a scalar-tensor theory (Fujii & Maeda, 2003; Riazuelo & Uzan, 2002), with the scalar field being the dark energy. These two classes of generalised quintessence models have been studied in detail in the linear regime in the literature (Bean, 2001a, b; Mangano et al., 2003; Clifton, Mota & Barrow, 2005; Nunes & Mota, 2006; Pettorino, Baccigalupi & Mangano, 2005; Koivisto, 2005; Brookfield et al., 2006; Koivisto & Mota, 2007; Mota & Shaw, 2007, 2006; Lee, Liu & Ng, 2006; Mota et al., 2007; Bean et al., 2008; Bean, Flanagan & Trodden, 2008; Boehmer et al., 2008, 2010).

In recent years, studies of the cosmological behaviour of the coupled quintessence model in the nonlinear regime have also been made, either via semi-analytical methods (Manera & Mota, 2006; Mota & van de Bruck, 2004; Mota, 2008; Shaw & Mota, 2008; Mota et al., 2008a; Mota, Shaw & Silk, 2008; Saracco et al., 2009; Wintergerst & Pettorino, 2010), or using -body simulation techniques (Maccio et al., 2004; Nusser, Gubser & Peebles, 2005; Kesden & Kamionkowski, 2006a, b; Springel & Farrar, 2007; Farrar & Rosen, 2007; Baldi & Pettorino, 2010; Hellwing & Juszkiewicz, 2009; Keselman, Nusser & Peebles, 2009, 2010; Hellwing, Knollmann & Knebe, 2010; Baldi, 2010; Baldi et al., 2010). In these studies the effect of the scalar field is generally approximated by a Yukawa-type ’fifth force’ or by a rescaling of the gravitational constant or the particle mass, without solving the scalar field equation explicitly. Very recently, Li & Zhao (2009, 2010); Zhao et al. (2010); Li & Barrow (2010) gave a new treatment and obtained an explicit solution to the scalar field perturbation on a spatial grid. The new results confirmed that the approximations adopted in the old literature were good for the models considered there (where the scalar potential was not very nonlinear), but for highly nonlinear potentials they broke down.

For the extended quintessence (more generally scalar-tensor) models, investigations using -body simulations are rarer. The work of Pettorino & Baccigalupi (2008), for example, outlined a recipe which uses certain approximation, such as a rescaling of gravitational constant, and does not solve the scalar field equation of motion explicitly. In Rodriguez-Meza et al. (2007); Rodriguez-Meza (2008a, b), the authors approximated the effect of scalar field coupling as a Yukawa force. However, none of these previous works tries to solve the scalar field on a mesh directly, and this is what we want to do in this work.

The aims of this work are threefold. Firstly, we want to develop the formulae and methods that are needed to solve the scalar field explicitly, which could serve as the basis for future work, and to find the regime of validity of our method. Secondly, we want to understand whether the approximations adopted in the previous studies are good or not; given the severe limits in the computing power; if those approximations do work well, then one does not need to resort to an exact scalar field solver, which is considerably more economical. Finally, we want to study structure formation in the nonlinear regime for some specific models, and investigate both the scalar field effects on the clustering of matter and the spatial variation of the gravitational constant (which is common to scalar-tensor theories).

The organisation of this paper is as follows: In Sect. 2 we list the basic equations which are needed in -body simulations and give their respective non-relativistic limits. To prevent the main text from expanding too much, some useful expressions are listed in Appendix A, and the discrete versions of the resulted equations are discussed and summarised in Appendix B. In Sect. 3 we briefly describe the numerical code we are using (relegating further details to Knebe, Green & Binney (2001); Li & Zhao (2010)), and the physical parameters of our simulations. We also present some results regarding the background cosmology and linear perturbation evolution in our models, which could be helpful in the understanding of the -body simulation results (our algorithm for the background cosmology is summarized in Appendix C). Sect. 4 contains the -body simulation results, including key structure formation observables such as nonlinear matter power spectrum, mass function and dark matter halo profile, as well as the spatial variation of the scalar field. It also includes several checks of the approximations made in the literature. We finally summarise and conclude in Sect. 5.

We use the unit unless explicitly restoring in the equations. The metric convention is . Indices run while run .

## 2. The Equations

This section presents the equations that will be used in the -body simulations, the model parameterisation and discretisation procedure for the equations.

### 2.1. The Basic Equations

We consider a general Lagrangian density for scalar-tensor theories

 L = 12κ∗[1+f(φ)]R−12∇aφ∇aφ+V(φ)−Lf, (1)

in which where is the (bare) gravitational constant, is the Ricci scalar, is the coupling function between the scalar field and curvature, the potential for and the Lagrangian density for fluid matter (baryons, photons, neutrinos and cold dark matter). Note that is a fundamental constant of the theory.

Varying the associated action with respect to metric yields the energy-momentum tensor of the theory (note the tilde, which is used to distinguish it from the defined below):

 ~Tab = Tfab+∇a∇bφ−12gab∇cφ∇cφ+gabV(φ) (2) −1κ∗[f(φ)Gab+(gab∇c∇c−∇a∇b)f(φ)]

where is the Einstein tensor, and is the energy-momentum tensor for matter (including baryons, dark matter, neutrinos and photons, which we collectively refer to as ’fluid matter’, although in -body simulations we use discrete particles rather than a fluid).

As usual, we can rearrange the Einstein equation as

 Gab = κ∗~Tab (3)

so that it now looks like

 Gab = κ∗1+fTfab−11+f(gab∇c∇c−∇a∇b)f (4) +κ∗1+f[∇aφ∇bφ−12gab(∇φ)2+gabV] ≡ κ∗Tab.

Note the difference between and ; throughout this paper, we will use a superscript for normal fluid matter, and quantities without a superscript always mean the total effective ones [the final line of Eq. (4)]. It is sometimes useful to define an effective Newton constant . Neither nor is the gravitational constant measured in a Cavendish-type experiment, which we denote instead by and is given by

 κ⨁ = κ∗1+f2+2f+4(dfd√κ∗φ)22+2f+3(dfd√κ∗φ)2 (5)

where is added to make dimensionless, which is the convention we shall always follow below. itself is obviously not a constant and we measure only it present-day value, .

Varying the action with respect to the scalar field, , gives the scalar field equation of motion

 ∇a∇aφ+∂V(φ)∂φ+R2κ∗∂f(φ)∂φ = 0 (6)

Since we will follow the motions of dark matter particles in the -body simulations, so we also need their geodesic equation. The dark-matter Lagrangian for a point particle with mass is

 LCDM(y)=−m0√−gδ(y−x0)√gab˙xa0˙xb0, (7)

where is the general coordinate and is the coordinate of the centre of the particle. From this equation we derive the corresponding energy-momentum tensor:

 TabCDM=m0√−gδ(y−x0)˙xa0˙xb0. (8)

Taking the conservation equation for dark matter particles (which, unlike in (Li & Zhao, 2009, 2010), does not couple to any other matter species, including the scalar field ), the geodesic equation follows as usual:

 ¨xa0+Γabc˙xb0˙xc0 = 0, (9)

where the second term on the left-hand side accounts for gravity.

Eqs. (4, 6 , 9) contain all the physics needed for the following analysis, though certain approximations and simplifications might have to be made in due course to make direct connection to -body simulations.

We will consider an inverse power-law potential for the scalar field,

 V(φ) = Λ4(√κ∗φ)α, (10)

where is a dimensionless constant and is a constant with dimensions of mass. This potential has also been adopted in various background or linear perturbation studies of scalar fields (either minimally or non-minimally coupled); the tracking behaviour its produces makes it a good dark energy candidate and for that purpose we shall choose . Meanwhile, the coupling between the scalar field and the curvature tensor is chosen to be a non-minimal one:

 f(φ) = γκ∗φ2, (11)

where is another dimensionless constant characterising the strength of the coupling. Note that here again is added into and to make a dimensionless quantity . Although the exact value of is unknown, so is and we can solve for instead of , not caring about the exact individual values of and .

### 2.2. The Non-Relativistic Limits

The -body simulation only probes the motion of particles at late times, and we are not interested in extreme conditions such as black hole formation and evolution, so we can take the non-relativistic limit of the above equations as a good approximation.

The existence of the scalar field and its coupling to the curvature leads to several possible changes with respect to the CDM paradigm:

1. The scalar field has its own energy-momentum tensor, which could change the source term of the Poisson equation because the scalar field, unlike the cosmological constant, can cluster (though the clustering is often quite weak in scalar field models). Also, unlike in coupled scalar field models, here the term will appear in the Poisson equation.

2. The background cosmic expansion rate is in general modified, and can either slow down or speed up the rate of structure formation.

3. The two gravitational potentials in the conformal Newtonian gauge metric , in which and are respectively the conformal time and comoving coordinate, are no longer equal to each other (as in general relativity), but are instead related by (see below).

It therefore becomes clear that the following two equations, in their non-relativistic forms, need to be solved in order to obtain the gravitational force on particles:

1. The scalar field equation of motion, which is used to compute explicitly the value of the scalar field at any given time and position;

2. The Poisson equation, which is used to determine the gravitational potential and force at any given time and position from the local energy density and pressure, which includes the contribution from the scalar field (obtained from the equation of motion).

Note that unlike in the coupled scalar field models, there is no fifth force because there is no direct coupling to the particles. The scalar coupling to the curvature, however, does modify the gravitational potential so that gravity no longer follows Einstein’s prescription and so this is a modified gravity theory.

We now describe these two equations in turn. For the scalar field equation of motion, we denote by the background value of and write . Then using the expressions given in Appendix A we write

 a2∇a∇aφ = φ′′+2Hφ′+→∇2xφ−2ϕφ′′ (12) −(ϕ′+3ψ′+4Hϕ)φ′

in which with the conformal time, is the derivative with respect to the comoving coordinate , and . Then, with the background part subtracted, Eq. (6) can be rewritten as

 δφ′′+2Hδφ′+→∇2xδφ+[V,φ(φ)−V,φ(¯φ)]a2 −2ϕ¯φ′′−(ϕ′+3ψ′+4Hϕ)¯φ′ +12κ∗[Rfφ(φ)−¯Rfφ(¯φ)]a2 = 0,

in which a bar denotes the background value, and the subscript denotes derivatives with respect to . Note that has the same sign as .

In our -body simulations we shall work in the quasi-static limit, i.e., we assume that the spatial gradients are much greater than the time derivatives, . Therefore, the time derivatives in the above equation are dropped and we obtain the simplified version

 c2→∂2x(aδφ) =

in which due to our sign convention , and we have restored the factor in front of (the here and in the remaining of this paper is times the in the original Lagrangian unless otherwise stated). Note that here has the dimension of mass density rather than energy density.

To complete Eq. (2.2), we still need expressions for and , which are again obtained using the quantities in Appendix A:

 R = −6a2a′′a(1−2ϕ) (14) +1a2[6ψ′′+6H(ϕ′+3ψ′)−4→∂2xψ+2→∂2xϕ], ¯R = −6a2a′′a (15)

and so

 Rfφ−¯R¯fφ ≐ ¯fφδR+¯Rδfφ (16) ≐ −1a2¯fφ(4→∂2xψ−2→∂2xϕ)−6a2a′′aδfφ

where we have again dropped time derivatives of and since they are small compared with the corresponding spatial gradients, and , .

Since only but not appears in the Poisson equation (shown below) , we also want to eliminate the in the scalar field equation of motion. This is easy in general relativity, because there we have the simple relation , which unfortunately no longer holds in scalar-tensor theories. However, we could use the components of the Einstein equation () to get a new relation between and . Noting that our -body simulations probe the very late time evolution (when radiation is negligible) when the only significant source for () is the scalar field, and

 ∇i∇jf = −1a2∂i∂jf  (i≠j) (17)

where , we could write the component of Einstein equation as

 ∂i∂j(ϕ−ψ) = −c21+f∂i∂jf

which gives approximately

 ∂i(ϕ−ψ) = −c21+f∂if

and so

 4→∂2xψ−2→∂2xϕ ≐ 2→∂2xϕ+41+f→∂2xf (18) ≐ 2→∂2xϕ+4¯fφ1+¯fc2→∂2xδφ.

It is important to note that in the second line of Eq. (18) we have implicitly linearised the equation; this is valid only if is not strongly nonlinear and . It turns out that the model considered in this work satisfies these criteria (). If either or is highly nonlinear, then we might have ; in that case we should not approximate to even in the coefficients of the perturbation variables such as here, or write . The reason for the latter stricture is as follows: if is highly nonlinear, then might change a lot even if fluctuates a little, implying that for the linearisation to apply on our spatial grid we need very small grid sizes which are impossible; moreover, it becomes complicated to decide which solution we should linearise around, as the values of in that area which we look at might be very different from the background value . The strategy for this situation is simple: instead of writing , we difference directly, because we know the value of in every grid cell. This will ensure no linearisation error. In what follows, however, we shall use Eq. (18), which causes negligible linearisation error but simplifies the equations a lot. We shall also write in the coefficients of perturbation quantities such as and .

Substituting Eqs. (16, 18) into Eq. (2.2) and rearranging, we complete the derivation of the scalar field equation of motion in the weak field limit, ending up with

 ⎡⎢⎣1+2¯f2φκ∗(1+¯f)⎤⎥⎦c2→∂2x(aδφ) (19) =

for our general Lagrangian Eq. (1) and

 [1+8γ2κ∗¯φ21+γκ∗¯φ2]c2→∂2x(a√κ∗δφ) (20) = −γ√κ∗¯φ→∂2xΦ−6γ(H′+H2)(a√κ∗δφ) −ακ∗Λ4a3⎡⎣1(√κ∗φ)1+α−1(√κ∗¯φ)1+α⎤⎦

for the model specified by Eqs. (10, 11), where .

Next consider the Poisson equation, which is obtained from the Einstein equation in the weak-field and slow-motion limits. Here we use the component of the Ricci curvature tensor, which is given as

 R0 0 = −3(a′′a−H2)(1−2ϕ) (21) +3ψ′′+3H(ψ′+ϕ′)+→∂2xϕ

using the expressions in Appendix A. According to the Einstein equations,

 R0 0 = κ∗2(ρTOT+3pTOT)a2 (22)

where and are the total energy density and pressure, respectively. Using these two equations and subtracting the background part (which is just the Raychaudhuri equation), it is straightforward to find that

 →∂2xΦ = κ∗a32[(ρTOT+3pTOT)−(¯ρTOT+3¯pTOT)]. (23)

in which we have dropped terms involving time derivatives of and , because they are much smaller than in the quasi-static limit. Using the energy-momentum tensor expressed in Eq. (4), the above equation can be rewritten as

 →∂2xΦ ≐ κ∗a32¯ρm(δ1+f−11+¯f)−¯fφ2(1+¯f)c2→∂2x(aδφ) −κ∗a3[V(φ)1+f−V(¯φ)1+¯f] +a(11+f−11+¯f)(κ∗¯φ′2+32f′′)

for the general Lagrangian Eq. (1) and

 →∂2xΦ = −γ√κ∗¯φ1+κ∗¯φ2c2→∂2x(a√κ∗δφ) −[κ∗Λ4a3(1+γκ∗φ2)(√κ∗φ)α−κ∗Λ4a3(1+γκ∗¯φ2)(√κ∗¯φ)α] +[(1+3γ)κ∗¯φ′2+3γκ∗¯φ¯φ′′]a ×[11+γκ∗φ2−11+γκ∗¯φ2]

for the model specified by Eqs. (10, 11). In these equations is the background density for matter, , and we have used the definition of given in Appendix C. We have also neglected the contribution from to the total density and pressure, because in the quasi-static limit we have and (which is confirmed by the -body simulation results111According to Eq. (20) we have , implying that , so neglecting time derivatives of is just like dropping time derivatives of and , which we have already done to obtain the modified Poisson equation.).

Finally, the equation of motion of the dark matter particles is the same as in general relativity

 ¨x+2˙aa˙x = −1a3→∇xΦ (26)

in which is determined by the modified Poisson equation Eq. (2.2). The canonical momentum conjugate to is so from the equation above we have

 dxdt = pa2, (27) dpdt = −1a→∇xΦ. (28)

Eqs. (20, 2.2, 27, 28) will be used in the code to evaluate the forces on the dark-matter particles and evolve their positions and momenta in time. But before applying them to the code we still need to switch to code units (see Sect. 2.3), further simplify them and create the discrete version (see Appendix B).

### 2.3. Code Units

In our numerical simulation we use a modified version of MLAPM ((Knebe, Green & Binney, 2001)), and we will have to change or add our Eqs. (20, 2.2, 27, 28) to it. The first step is to convert the quantities to the code units of MLAPM. Here, we briefly summarise the main results.

The (modified) MLAPM code uses the following internal units (where a subscript stands for ”code”):

 xc = x/B, pc = p/(H0B) tc = tH0 Φc = Φ/(H0B)2 ρc = ρ/¯ρ, u = ac2√κδφ/(H0B)2, (29)

where denotes the comoving size of the simulation box, is the present Hubble constant, and is the matter density. In the last line the quantity is the scalar field perturbation expressed in terms of code units and is new to the MLAPM code.

In terms of , as well as the (dimensionless) background value of the scalar field, , some relevant quantities are expressed in full as

 V(φ) = Λ4(√κ¯φ+B2H20ac2u)α, f(φ) = 1+γ(√κ¯φ+B2H20ac2u)2, Vφ(φ) = −α√κΛ4(√κ¯φ+B2H20ac2u)1+α, fφ(φ) = 2γ√κ(√κ¯φ+B2H20ac2u), (30)

and the background counterparts of these quantities can be obtained simply by setting (recall that represents the perturbed part of the scalar field) in the above equations.

We also define

 λ ≡ κΛ43H20, (31)

which will be used frequently below.

Making discrete versions of the above equations for -body simulations is then straightforward, and we refer the interested readers to Appendix B to the whole treatment, with which we can now proceed to do -body simulations.

## 3. Simulation Details

### 3.1. The N-Body Code

Some of our main modifications to the MLAPM code for the coupled scalar field model are:

1. We have added a solver for the scalar field, based on Eq. (B7). It uses a nonlinear Gauss-Seidel scheme for the relaxation iteration and the same criterion for convergence as the default Poisson solver. But it adopts a V-cycle instead of the self-adaptive scheme in arranging the Gauss-Seidel iterations.

2. The value of solved in this way is then used to calculate the total matter density, which completes the calculation of the source term for the Poisson equation. The latter is then solved using a fast Fourier transform on the domain grids and self-adaptive Gauss-Seidel iteration on refinements.

3. The gravitational potential obtained in this way is then used to compute the force, which is used to displace and kick the particles.

There are a lot of additions and modifications to ensure smooth interface and the newly added data structures. For the output, as there are multilevel grids all of which host particles, the composite grid is inhomogeneous and so we choose to output the positions and momenta of the particles, plus the gravity and values of and at the positions of these particles. We also output the potential and scalar field values on the domain grid.

### 3.2. Physical and Simulation Parameters

The physical parameters we use in the simulations are as follows: the present-day dark-energy fractional energy density and ,  , , . Our simulation box has a size of  Mpc, where . We simulate four models, with parameters , , and respectively. In all those simulations, the mass resolution is ; the particle number is ; the domain grid is a cubic and the finest refined grids have 16384 cells on each side, corresponding to a force resolution of about kpc.

We also run a CDM simulation with the same physical parameters and initial condition (see below).

### 3.3. Background and Linear Perturbation Evolution

Since the coupling between the scalar field and the curvature produces a time-varying effective gravitational constant, and the scalar field contributes to the total energy-momentum tensor, we expect that cosmology in the extended quintessence models is generally different from CDM at the background and linear perturbation levels. A good understanding of this will be helpful in our analysis of the results from -body simulations, and this is the subject of this subsection.

Our algorithm and formulae for the background cosmology are detailed in Appendix C, and are implemented in MAPLE. We output the relevant quantities in a predefined time grid, which could be used (via interpolation) in the linear perturbation and -body computations.

Fig. 1 shows the time evolutions of some background quantities of interests. For ease of comparison we have chosen and to be the same in all models including the CDM one (for definitions of and see Appendix C), and as a result in the upper left panel the curves for different models converge at common righthand ends. We see increasing results in an earlier and slower growth of (). This indicates a larger dark energy equation of state parameter, , which is confirmed by the upper right panel. Physically, this is because, the larger is, the steeper the potential becomes and thus the faster the scalar field rolls. Notice that is also larger for positive , with being the same. This is because in Eq. (6) the Ricci scalar and for positive the term has the same sign as , thus helping the scalar field to roll faster. Because of its large predicted value of , the model is already excluded by cosmological data, but here we shall keep it for purely theoretical interest (i.e., to see how changing or changes the nonlinear structure formation).

We are also interested in how the expansion rate in an extended quintessence model differs from that in CDM, and the results for our models are shown in the lower-left panel of Fig. 1, which plots the as a function of . The rather odd behaviour of the models at low redshift is because of the complicated evolution of the scalar field (and the fact that we have chosen to be the same for all models, again for ease of comparison), while the high-redshift behaviour could be seen directly from Eq. (C4). In Eq. (C4) the energy density of the scalar field can be dropped at high , and so we have

 (HH0)2 ≈ 1+f01+fΩma−1 (32)

where we have also neglected the radiation for simplicity (which is valid after the matter-radiation equality). This shows that in extended quintessence models the gravitational constant relevant for the background cosmology is rescaled by . Because where is the present-day value of , and is monotonically increasing in time, so for our choice of [cf. Eq. (11)] we have for and for : thus models with have .

It turns out that the gravitational constant relevant for the growth of matter density perturbations is also different from the one governing the background cosmology. If we denote the matter density perturbation by , then it can be shown, using the linear perturbation equations, that on small scales the evolution equation for reduces to

 δ′′m+Hδ′m = GN3H202Ωmδma2 (33)

in which and is the conformal time (see Appendix C), and we have defined

 GN ≡ 1+f01+f2+2f+4(dfd√κ∗φ)22+2f+3(dfd√κ∗φ)2. (34)

Note that this quantity could also be directly read off from the modified Poisson equation Eq. (B4).

In the lower right panel of Fig. 1 we display the evolution for in the models considered. Again, is larger at earlier times for positive and smaller for negative , because of our specific choice of in Eq. (11), and the fact that is always increasing in time.

It is well known that a higher rate of background expansion means that structures have less time to form, and a larger speeds up the structure formation. These two effects therefore cancel each other to some extent, which results in a weaker net effect of an extended quintessence field on the large scale structure formation. This is confirmed by our linear perturbation computation depicted in Fig. 2. In the right-hand panels of this figure we have plotted the matter power spectra for different models at two different redshifts (0 and 49). It is interesting to note that on small scales the matter power is closer to that of CDM, despite the significant differences in background expansion rate and (cf. Fig. 1). Because of this, we shall choose CDM initial condition for our -body simulations for all our models, saving the effort of generating separate initial conditions for different models.

The left hand panels of Fig. 2 display the CMB power spectra for the models we consider. Again the difference from CDM is fairly small, and there is only a small shift of the CMB peaks even though the background expansion rate changes quite a bit. The latter is because peak positions are determined by the ratio of the sound horizon size at decoupling and the angular distance to the decoupling, and in our model both of these decrease/increase as the Universe expands faster/more slowly, their ratio does not change much.

To briefly summarise, the study of background cosmology and linear perturbation shows that a modified background expansion rate and a rescaled gravitational constant, the two most important factors affecting structure formation in extended quintessence models are opposite effects. It is then of interest to see how these two effects compete in the nonlinear regime.

## 4. N-body Simulation Results

This section lists the results of extended quintessence -body simulations. We shall start with a few preliminary results which both give some basic idea about the extended quintessence effects and serve as a cross check of our codes. Then we discuss the key observables for the nonlinear structure formation such as matter power spectrum, mass function and halo properties. We also comment on the halo profile of the scalar field and the spatial variation of gravitational constant.

### 4.1. Preliminary Results

As mentioned above, in both the linear and -body codes we compute background quantities via an interpolation of some pre-computed table. Because background cosmology is important in determining the structure formation, it is important to check its accuracy. For this we have recorded in Table 1 The age of the universe today for different models as computed by these two codes. The two codes are compatible with each other indeed.

model linear code -body code
CDM 13.680 13.678
13.639 13.638
13.408 13.408
13.513 13.513
12.097 12.096

Because one of the advantages of our -body code is that it solves the scalar field perturbation explicitly, it is important to check that the solution is with expectations. From Eqs. (B1, B2) it could be seen clearly that, if the contribution to the local density and pressure from the scalar field is negligible compared with that from matter, then the modified Poisson equation and scalar field equation of motion end up with the same source term (up to a -dependent coefficient). In this situation we expect

 u = −2γ√κ∗¯φ1+8γ2κ∗¯φ21+γκ∗¯φ2Φc, (35)

which means that is simply proportional to with a time-dependent coefficient. In Fig. 3 we have checked this relation explicitly: we select a thin slice of the simulation box, fetch the values for and at the positions of the particles (about 10000 in total) therein, and display them as scatter plots. The solid curve is the approximation Eq. (35) while the green dots are simulation results; we can see they agree very well with each other, showing that the above approximation is a good one. Note that the scalar field perturbation is generally less than , compared with the background value . This confirms that it is consistent to neglect the perturbation in scalar field density/pressure, drop terms such as and , and replace by in coefficients of perturbation quantities such as and . It also serves as a check of the numerical code.

As a final consistency check, let us consider the total gravitational force on particles. In extended quintessence models, this is given by Eq. (B2), and when the perturbation in the scalar field density/pressure is negligible (which is the case as shown above) we get

 ∇2Φc ≈ 32GNΩmH20(ρc−1) (36)

in which is given in Eq. (34). On the other hand, if we consider (naïvely) that gravity is described by general relativity, then we should neglect the on the right-hand side. Manipulating Eqs. (B1, B2) we obtain:

 1+γκ∗¯φ201+γκ∗¯φ2∇2(Φc+γ√κ∗¯φ1+γκ∗¯φ2u) ≈ 32ΩmH20(ρc−1).

Thus acts as the potential for naïve gravity (i.e., general relativity), and by differcing it we could obtain the naïve gravitational force. In Fig. 4 we show the scatter plot of the naïve gravity versus full gravity for the same particles as in Fig. 3 (green dots) as well as their approximate ratio (solid line). Again, the agreement is remarkably good.

### 4.2. Nonlinear Matter Power Spectrum

As we have seen above, the linear matter power spectrum for the extended quintessence model really does not show much useful information on small scales, and so we need to investigate whether nonlinear effects could change this situation and therefore potentially place more meaningful constraints.

Fig. 5 provides a positive answer to this question. Here we have plotted the fractional difference of the extended quintessential nonlinear matter power spectrum from that for CDM (remember that we use the same initial condition for all simulations). We can see that for the models with the differences are small even in the nonlinear regime, indicating that the scalar field really does not affect the matter distribution significantly if the potential is flat. However, for the cases in which the coupling strength remains the same, the difference could be as large as , guaranteeing an observable signature.

Furthermore, for negative (the purple curve) the extended quintessential power spectrum beats the CDM one on small scales, whereas for the positive case (the pink curve) it is just the opposite. As shown before, when , both the background expansion rate and the effective gravitational constant governing the structure formation decrease, boosting and weakening the collapse of matter respectively. In our cases the first effect has clearly taken over on small scales.

### 4.3. Mass Function

A second important observable is the mass function. This gives the number density of dark matter halos as a function of halo mass. For this we need to identify the dark matter halos from the output particle distribution of the -body simulations, and this determination is performed using a modified version of MHF (Knebe & Gibson, 2004), MLAPM’s default halo finder.

MHF optimally utilizes the refinement structure of the simulation grids to pin down the regions in which potential halos reside and organize the refinement hierarchy into a tree structure. MLAPM refines grids according to the particle density on them and so the boundaries of the refinements are simply isodensity contours. MHF collects the particles within these isodensity contours (as well as some particles outside). It then performs the following operations: (i) assuming spherical symmetry of the halo, calculate the escape velocity at the position of each particle, (ii) if the velocity of the particle exceeds then it does not belong to the virialized halo and is removed. Steps (i) and (ii) are then iterated until all unbound particles are removed from the halo or the number of particles in the halo falls below a pre-defined threshold, which is 20 in our simulations. Note that the removal of unbound particles is not used in some halo finders using the spherical overdensity (SO) algorithm, which includes the particles in the halo as long as they are within the radius of a virial density contrast. Another advantage of MHF is that it does not require a predefined linking length in finding halos, such as the friend-of-friend procedure.

Our modification to MHF is simple: because the effective gravitational constant in the extended quintessence models is rescaled by a factor [cf. Eq. (34)], the escape velocity of particles from a halo is also multiplied by this factor, and in MHF we have only changed the criterion for removing particles from virialised halos accordingly. In reality, because we are only interested in the halos in this work, is quite close to 1 and the effect of our modification is not large.

The mass functions for our simulated models are shown in Fig. 6. It shows that all extended quintessence models considered here, irrespective of their parameters, produce less massive halos than CDM, whereas (only) the model produces a larger number of less massive halos. These features are in broad agreement with those shown in the matter power spectra (Fig. 5) where all models show less matter clustering on the large scales, whereas (only) the model shows more power on small scales. The physical reason is again the competition between the modified background expansion rate and rescaled effective gravitational constant .

### 4.4. Halo Properties

In the CDM paradigm, it is well known that the internal density profiles of dark matter halos are very well described by the Navarro-Frenk-White (Navarro, Frenk & White, 1996) formalism

 ρ(r)ρc = βrRs(1+rRs)2 (37)

where is the critical density for matter, is a dimensionless fitting parameter and a second fitting parameter with length dimension. and are generally different for different halos and should be fitted for individual halos, but the formula Eq. (37) is quite universal.

We are thus interested in whether the halo profiles in an extended quintessential Universe are also featured by this universal form. For this we select the 80 most massive halos from each simulation and fit their density profiles to Eq. (37). The results show that the NFW profile describes the extended quintessential halos at least as well as it does for the CDM halos. Fig. 7 shows the fittings for two halos randomly picked out of the 80: one at Mpc with mass and the other at Mpc with a mass .

There are some interesting features in Fig. 7. Firstly, for the models with (the top panels) the halo density profile for extended quintessence models (green asterisks) is very similar to the CDM results (black crosses) and thus their fittings almost coincide. Secondly, for the model of , the chosen halos show more concentration of the density profiles in the scalar model than in CDM. Thirdly, the model of has just the opposite trend and suffers a suppression of density in large parts of chosen halos.

To verify that the above features are actually typical for the corresponding models, we have plotted in Fig. 8 the fitting results for all the 80 massive halos in all simulated models. Here in addition to the NFW concentration parameter , where is the radius at which the density is equal to 200 times the critical density and the NFW parameter, we have also shown the fitting errors for each halo.

We would like to point out several important implications of Fig. 8. Firstly, for all models the fitting error for the extended quintessential halos (lower green asterisks) is comparable to that for the CDM halos (lower black crosses), indicating that the density profiles for the former are equally well described by the NFW formula Eq. (37). Secondly, for the models with (the top panels) we can see that the fitted for the extended quintessential halos is comparable to that for CDM, which is in agreement with our finding in Fig. 7 that the density profiles for the chosen halos are almost the same as in the CDM prediction. Thirdly, for the model of , the halos tend to be more concentrated (i.e., with larger ) than in CDM. Fourthly, for the model , the halos tend to be less concentrated (i.e., with smaller ) than in CDM. The above three features show that our qualitative findings in Fig. 7 are quite typical. Finally, the halo masses in the model of are on average smaller than those in CDM, because the upper green asterisks in the lower right panel consistently shift leftwards with respect to the upper black crosses: this is consistent with the mass function result that this model produces less massive halos than CDM.

In summary, the halo density profiles for the extended quintessence models are well described by the NFW formula, but the existence of the scalar field and in particular its coupling to curvature do change the concentration parameters of the halos, so long as the potential is not too flat. It seems that the modified background expansion rate beats the effect of the rescaled effective gravitational constant here.

### 4.5. Halo Profile for Scalar Field Perturbation

We have already seen that the coupling between the scalar field and the curvature scalar causes time and spatial variations of the locally measured gravitational constant . It is then of our interest to ask how varies across a given halo and whether this could produce observable effects. This subsection answers this question, by giving an analytical formula and comparing it with numerical results.

Recall that Fig. 3 shows that to a high precision the scalar field perturbation is proportional to the gravitational potential [cf. Eq. (35)] everywhere. This means that if we could derive an analytical formula for in halos, then we know straightforwardly. Such a derivation has been done in (Li, Mota & Barrow, 2010) for a different model, but here we shall briefly repeat it for the extended quintessence model for completeness.

Assuming Eq. (37) as the density profile and sphericity of halos, we can derive , the circular velocity of a particle moving around the halo at a distance from halo centre, to be

 V2c(r) = GM(r)r (38) = 4πGβρcR3s[1rln(1+rRs)−1Rs+r]

where is the mass enclosed in radius , is the properly rescaled gravitational constant. Again, this equation is parameterized by and . From a simulation point of view, it is straightforward to measure and then use Eq. (38), instead of Eq. (37), to fit the values of and ; from an observational viewpoint, it is easy to measure , which could again be used to fit and .

The potential inside a spherical halo is then given as

 Φ(r) = ∫r0GM(r′)r′2dr′+C (39)

in which is the gravitational force and is a constant to be fixed using the fact that where is the value of the potential far from the halo.

Using the formula for given in Eq. (38) it is not difficult to find that

 ∫r0GM(r′)r′2dr′ =

and so

 C = Φ∞−4πGβρcR2s. (40)

Then it follows that

 Φ(r) = Φ∞−4πGβρcR3srln(1+rRs). (41)

If the halo is isolated, then and we get

 Φ(r) = −4πGβρcR3srln(1+rRs). (42)

However, in -body simulations, we have a large number of dark matter halos and no halo is totally isolated from the others. In such situations,