{'fraction_non_alphanumeric': 0.03770874767181236, 'fraction_numerical': 0.007339710199829529, 'mean_word_length': 4.455054244876873, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structure Formation and Backreaction in Growing Neutrino Quintessence

A dependence of the neutrino masses on the dark energy scalar field could provide a solution to the why now problem of dark energy. The dynamics of the resulting cosmological model, growing neutrino quintessence, include an attractive force between neutrinos substantially stronger than gravity. We present a comprehensive approach towards an understanding of the full cosmological evolution including the formation of large-scale neutrino structures. Important effects we account for are local variations in the dark energy and the backreaction on the background evolution, as well as relativistic neutrino velocities. For this aim, we develop a relativistic N-body treatment of the neutrinos combined with an explicit computation of the local quintessence field. At its current stage, the simulation method is successful until z ~ 1 and reveals a rich phenomenology. We obtain a detailed picture of the formation of large-scale neutrino structures and their influence on the evolution of matter, dark energy, and the late-time expansion of the universe.


I. INTRODUCTION
In today's cosmological standard scenario, a cosmological constant Λ is assumed to explain the observed accelerated expansion of the Universe. Although so far consistent with all major observational probes [1][2][3][4], it faces two fundamental and unresolved problems. Besides its disturbingly tiny value (the cosmological constant problem), it remains miraculous why Λ has become important just recently (the why now or coincidence problem).
It is thus tempting to think of alternative cosmological scenarios in which these two problems are alleviated. In this work, we study growing neutrino quintessence, which addresses both of the problems. Growing neutrino quintessence [5,6] relies on two assumptions. First, the dark energy component is described by a dynamical scalar field. Second, the neutrino masses depend on this field. Studying the background evolution, it is found that the coincidence problem can indeed be solved for a neutrino-cosmon coupling somewhat larger than gravity.
In order to explore the implications of this model, its evolution has to be understood also on the perturbation level. Linear perturbation theory, however, breaks down even at large scales [7]. This is due to large overdensities in the neutrino fluid becoming non-linear at z 2 on supercluster scales.
An understanding of the non-linear evolution, although of utmost importance for confronting the model with observational data, is still lacking. In first attempts, some aspects of the model have been studied with hydrodynamical and N -body methods [8,9]. Yet, three crucial aspects of the model have not been included so far. First, as shown by an analytical study of single non-linear neutrino overdensities, local variations of the neutrino mass can become very large [10]. Second, numerical results [9] show that neutrino particles are likely to reach highly relativistic velocities and thereby leave the Newtonian limit. * m.weber@thphys.uni-heidelberg.de Third, the local mass variations as well as the relativistic corrections change the average energy-momentum tensor of neutrinos; this backreaction effect alters the background evolution of the dark energy scalar field.
It is the purpose of this work to form the basis for an adequate N -body based simulation method capable of including all major effects of growing neutrino quintessence. For the first time, we implement an explicit computation of the local dark energy scalar field for every time step. This allows us to account for local mass variations and backreaction effects. We present a relativistic treatment of the neutrino dynamics as appropriate for the high velocities of neutrinos occurring during the structure formation process. The results of linear perturbation theory are used for the initial conditions. All these aspects are not merely small corrections, but decisive for every quantitative prediction.
The organization of this paper is as follows. We present growing neutrino quintessence and derive the equations of motion required for our simulation in Sec. II. In Sec. III, we discuss strategies how to account for the peculiar features of growing neutrino quintessence in a numerical simulation. The numerical results are presented in Secs. IV and V. At first, we have a look at the formation of neutrino structures (Sec. IV), and then explore the effects on the evolution of dark energy and matter (Sec. V). We summarize in Sec. VI.

A. Overview
Growing neutrino quintessence is a possible solution to both the cosmological constant and the why now problem of dark energy [5,6]. Since it is a quintessence model with a dark energy scalar field, the cosmon ϕ, dark energy evolves dynamically. During most of the cosmological evolution, its energy density decays similarly to the densities of the other species. The vanishing of the cos-mon potential V (ϕ) for large values of ϕ can be rooted in the approach to a fixed point with effective dilatation symmetry. The scaling solution with dark energy of the same order of magnitude as radiation or smaller can then explain the tiny overall size of the dark energy density by the large age of the Universe.
In contrast to standard quintessence scenarios, however, a coupling between the cosmon ϕ and the neutrinos provides a solution to the coincidence problem as well. Due to their small masses, cosmological neutrinos have become non-relativistic only recently. This event triggers the onset of dark-energy domination in growing neutrino quintessence.
A coupling between the cosmon ϕ and a matter species, here assumed to be the neutrinos, is expressed by the exchange of energy and momentum. The individual energy-momentum tensors do not satisfy a conservation equation, only their sum does. Denoting with T αβ the energy-momentum tensor of the neutrinos and with S αβ the energy-momentum tensor of the cosmon field, we have Since no known symmetry requires Q = 0, we generally have to expect a non-vanishing coupling. A specific form of the coupling proposed by early works [11,12] is with a dimensionless coupling parameter β and T ≡ T α α , using units where the reduced Planck mass M P = (8πG) −1/2 is set to unity. Writing T = −ρ ν + 3p ν , we see that T vanishes as long as the neutrinos are relativistic, w ν = p ν /ρ ν = 1/3. The coupling only becomes effective once the neutrinos have become non-relativistic. For large values of the coupling β, the coupling can stop the evolution of the cosmon. An almost constant value of the cosmon potential leads then to an onset of dark-energy domination at recent times, similar to the concordance model ΛCDM.
On a particle physics level, the coupling, Eq. (2), is realized as a dependence of the neutrino masses on the cosmon field. The coupling parameter β, in terms of the average neutrino mass m ν ≡ m ν (ϕ) [5], now is We further assume that the coupling parameter β is constant, i. e., it does not depend on ϕ. In this case, we have the crucial relation In scenarios with an expansion history similar to ΛCDM, the coupling β takes large negative values, typically of order β ∼ −10 2 . The possible couplings to other matter species are assumed to be negligible. The self-interaction of the cosmon is given by a potential V (ϕ), for which an example is the exponential potential [5], where α is a dimensionless model parameter with typical values α 10 to satisfy early dark energy constraints [13,14].
The coupling, strong enough to stop the cosmon evolution in the background, also crucially modifies the evolution of perturbations. In fact, if the neutrinos are non-relativistic, the cosmon perturbation δϕ is approximately a factor of 2β larger than the neutrino-induced gravitational potential. The perturbation δϕ describes an attractive force between the neutrinos of order |F | ≈ |β∇δϕ| ≈ 2β 2 |F gravity | [8]. Since realistic scenarios have β 2 ≫ 1, the extra force leads to a very rapid growth of perturbations in the neutrino fluid becoming non-linear at z nl ≈ 1-2, even on large scales [7].
Analyzing growing neutrino quintessence quantitatively thus requires adequate methods to study the nonlinear evolution. Perturbations δϕ in the quintessence field imply local mass variations δm ν by virtue of Eq. (4). These variations can significantly change the averaged energy-momentum tensor entering the background equations. This backreaction effect of structure formation on the background evolution will turn out to be crucial. Furthermore, the forces during the non-linear evolution are strong enough to accelerate the neutrinos to relativistic velocities. We thus have to work with a relativistic description. For these aims, we shall next derive and collect the necessary equations.

B. Fundamental equations
In a first step, we will present a fundamental definition of the model in terms of an action and derive the basic equations, Eqs. (1) and (2), introduced above. In principle, one has to describe all three neutrino species together. This can be done in growing neutrino quintessence [6]. In what follows, however, we adopt the simplified working hypothesis of the neutrinos having all the same mass. The three flavors of neutrinos enter then only in the initial number density of neutrinos. For our simulation, we can deal effectively with one species.
The dynamics of the cosmon field ϕ and the neutrino field ψ are described by a usual scalar field Lagrangian L ϕ and a Majorana Lagrangian L ν respectively, with the peculiarity that the neutrino mass term m νψ ψ is assumed to depend on ϕ, Here, we use the vierbein formalism to describe the neutrino field in curved spacetime [15]. The quantities γ α (x) are related to the usual Dirac matrices γ a (a = 0, 1, 2, 3) by virtue of the vierbein e α a (x): γ α (x) = γ a e α a (x), where the vierbein is related to the metric by g αβ = e α a e β b η ab , η ab = diag(−1, 1, 1, 1). The action reads where L 0 includes the remaining cosmological species and gravity. A Majorana constraint relatesψ to ψ. We obtain the energy-momentum tensor T αβ for the neutrino field by the usual definition We assume a diagonal metric, which we have, e. g., in the Newtonian gauge-otherwise, one employs a generalization with vierbeins. The energy-momentum tensor reads where we have symmetrized in the indices α and β (cf., e. g., [15]). The equations of motion following from L ν are the Dirac equations in curved spacetime, In the absence of the coupling, i. e. in the case m ν = const., these equations imply the energy-momentum conservation equation ∇ β T αβ = 0 for the neutrino field. In growing neutrino quintessence, however, the mass depends on the spacetime coordinates x via ϕ(x), and the derivative also acts on m ν . Inserting the equations of motion into ∇ β T αβ , we find We will give an interpretation ofñ ≡ −iψψ below. Using the definition of β, Eq. (3), we may write ∂ α m ν (ϕ) = −βm ν (ϕ)∂ α ϕ. In our discussion, an important role is played by the trace T of the neutrino energy-momentum tensor, This eventually shows the form of the energy-momentum exchange, cf. Eq. (2), A standard computation for the energy-momentum tensor of the cosmon, leads to and one verifies the conservation equation ∇ β T αβ + S αβ = 0 for the sum of the neutrino and cosmon energy-momentum tensors.
In the non-relativistic limit,ñ ν introduced above corresponds to the neutrino number density n ν . In contrast toñ ν , however, n ν does not transform as a scalar. For a Lorentz transformation, e. g., n ν picks up the volume contraction factor 1/γ. This will become more concrete when we represent the neutrino field by effective relativistic particles (see Sec. III B).

C. Cosmon dynamics
The coupling between neutrinos and the cosmon field has important impacts on the evolution of both species. Variation of the action with respect to ϕ yields the modified Klein-Gordon equation We split into background quantities (spatial averages only depending on time) and perturbations (spatially varying with vanishing mean), g αβ =ḡ αβ + δg αβ , ϕ = ϕ + δϕ, and T αβ =T αβ + δT αβ . We assume that the metric perturbations δg αβ and the cosmon perturbation δϕ can be treated in linear approximation and that their time derivatives are small. We choose the conformal Newtonian gauge (see, e. g., [16]), in which the metric reads with the conformal time η, the comoving coordinates x, the scale-factor a(η), and the two scalar potentials Ψ and Φ. It is now straightforward to evaluate Eq. (17). Separating into a background part, independent of the spatial position, and a linear perturbation part, we obtain two evolution equations. For the background quantities, we find where primes denote derivatives with respect to η, H ≡ a ′ /a is the conformal Hubble parameter, and we have usedT = −ρ ν + 3p ν . Obviously, the coupling is ineffective as long as w ν =p ν /ρ ν = 1/3, i. e. if the neutrinos are relativistic. Once the neutrinos turn non-relativistic, the right-hand side of Eq. (19) stops the further evolution ofφ. By this mechanism, the model can solve the coincidence problem [5,6].
In the perturbations, the equation reads where the spatial derivatives refer to comoving coordinates. In the fluid description, δT = −δρ ν + 3δp ν . For our purpose, however, it is more convenient to calculate δT directly from the distribution of particles in our simulation.

D. Neutrino dynamics
In this section, we investigate the motion of a neutrino particle with a cosmon-depending mass m ν (ϕ). We describe a neutrino as a classical, yet relativistic particle with world line ξ α (τ ) and four-velocity u α = dξ α /dτ , where τ denotes the particle's proper time, defined via dτ 2 = −g αβ dξ α dξ β . The energy-momentum tensor of this particle is given by where g is the determinant of the metric and δ 4 (x) denotes the four-dimensional Dirac delta function. The factor 1/ √ −g ensures the correct normalization of the Dirac delta function in curved spacetime by compensating the invariant volume form √ −g d 4 x. The one-particle action is constructed from the energy-momentum tensor, The equations of motion can be obtained by varying S ν with respect to the particle's path ξ(τ ). In the uncoupled case, m ν = const., this would give the standard geodesic equation. The modifications due to the cosmon-neutrino coupling are the same that we will find below by using energy-momentum conservation.
In what follows, we will use Eq. (21) in the energymomentum conservation equation, Eq. (14), in order to derive the equations of motion. Modifications to the standard geodesic equation enter in two ways, on the righthand side through the exchange term Q α = −β T ∂ α ϕ and on the left-hand side through the cosmon-depending mass m ν (ϕ).
Since u α u α = −1, the right-hand side simply becomes The left-hand side is It is straightforward to show The derivative acting on m ν (ϕ) is ∂m ν (ϕ)/∂τ = −β m(ϕ) u λ ∂ λ ϕ. Apart from this contribution, i. e. in the uncoupled case, Eq. (24) reproduces the standard geodesic equation. Thus, we find Comparing this with Eq. (23), we arrive at the equation of motion: which describes the deviation from the geodesic motion due to the cosmon-neutrino coupling. This equation of motion can also be obtained by a conformal transformation of the geodesic equation [9]. Let us give a brief interpretation of its terms.
• Γ α ρσ u ρ u σ describes the gravitational effects. In the component α = 0, this is the Hubble damping. The spatial components give rise to the gravitational force. In the Newtonian limit, this would be given by ∇Ψ. In the relativistic case, the situation is more complicated involving both potentials, Ψ and Φ.
• β ∂ α ϕ is the cosmon-mediated fifth force. In the Newtonian limit, this corresponds to an attractive force between neutrinos about 2β 2 times stronger than gravity [8].
• β u λ ∂ λ ϕ u α represents a velocity dependent force. It can be understood as a consequence of momentum conservation. Particles are accelerated when they move into a direction where they lose mass. It also modifies the universal damping term, which is no longer proportional to the Hubble parameter H, but to (H − βϕ ′ ).
It is instructive to consider the case of a static particle, u i = 0. Then, the time component of the right-hand side of Eq. (27) vanishes, and the spatial components reduce to β ∂ i ϕ. Thus, a static particle is only affected by the spatial variation ∇ϕ and not by ϕ ′ , and u 0 = (1 − Ψ)/a does not depend on ϕ.

III. STRATEGY AND METHOD
Although many efforts to describe the structure formation process in growing neutrino quintessence have been made (cf., e. g., [7][8][9]17]), none of the applied methods was capable of drawing a comprehensive picture, nor of providing reliable quantitative results once nonlinearities are strong.
These challenges motivate the development of a new method, specifically designed for growing neutrino quintessence. We describe the framework in Sec. III A. As a fundamental ingredient, we keep track of the perturbed cosmon scalar field allowing us to incorporate local mass differences of the neutrinos. We find a significant neutrino mass suppression once a large fraction of the neutrinos is bound in non-linear structures. This induces a backreaction on the evolution of the cosmological background (Sec. III C). The initial conditions are drawn from the linear results (Sec. III D). We discuss numerical issues in Sec. III E.

A. Framework
An important goal is to make predictions for the density contrasts δ ν (x), δ m (x) and the peculiar velocity fields v pec ν (x), v pec m (x) of neutrinos and matter respectively. These are linked to the gravitational potential and to various observables. These fields, however, do not carry all the necessary information needed to describe their evolution. They merely arise from the first moments of the full phase-space distribution functions f ν (η, In a non-linear evolution, higher moments of these distributions become important. It is thus most efficient to directly sample the distribution functions by a finite number of effective particles, N ν and N m respectively. These particles carry a comoving position x, a velocity v = dx/dη, and a rest mass M ν (neutrinos), or M m (matter). Our method thus bases upon an N -body scheme.
The motions of these particles depend on the cosmological background, the two gravitational potentials Ψ(x), Φ(x), and the cosmon ϕ(x) =φ(η) + δϕ(x). In the Newtonian limit and neglecting the local mass variation of neutrinos, both gravity and the fifth force can be described as two-body forces. If we want to go beyond these approximations, we need to know the explicit values of these fields.
We model the fields Ψ, Φ, and δϕ as discrete values on a three-dimensional grid in a volume V with N c cells. The volume is a cube with side length L, every cell has the volume (∆x) 3 = V /N c . We employ periodic boundary conditions. For the concrete values chosen, cf. Table I  Let us outline the basic steps of the algorithm we use. The details will be given in subsequent sections.
(a) Generate a realization of initial fields δ ν , δ m , v pec ν , v pec m from the linear spectra. (b) Sample the corresponding distribution functions with N ν effective neutrino and N m effective matter particles.
(a) Accelerate and move the effective neutrino and matter particles (cf. Secs. III B and II D).

B. Particles and fields
As stated above, we assume that the fields Ψ, Φ, and δϕ can be treated linearly and that their time derivatives, Ψ ′ , Φ ′ , and δϕ ′ , are small. This leads to simple algebraic equations in Fourier space that allow us to calculate the fields from a given distribution of particles. This method requires to carry out several Fourier transforms at each time step, which can be done efficiently by the use of Fast Fourier Transform routines. Let us collect the relevant equations.
The gravitational potential Φ, on subhorizon scales, is obtained from the Poisson equation [16], where k = |k| and δρ = δρ ν + δρ m + δρ ϕ . While we can obtain δρ ϕ from the cosmon field and its perturbation, we need to calculate δρ ν (and δρ m analogously) from the effective particles.
Since δρ ν = −δT 0 0 , we get the contribution of one particle at position ξ from Eq. (21), where we have introduced the Lorentz factor, and the determinant of the spatial metric, In the case of matter, we can approximate v 2 ≪ 1. In contrast, we keep the full relativistic equations for neutrinos since values of v 2 close to one can be reached once large neutrino structures form.
Finally, the discrete field value ρ ν (x) at a cell x is obtained by summing up the contributions of all the particles located inside the cell. Subtracting the mean ρ ν yields δρ ν (x). The Poisson equation, Eq. (28), can now be solved as follows. The sum δρ(x) = δρ ν (x) + δρ m (x)+ δρ ϕ (x) is Fourier transformed and inserted into Eq. (28). The resulting field Φ(k) and its gradient ikΦ(k) are then transformed back to real space, yielding Φ(x) and ∇Φ(x).
The second gravitational potential Ψ equals Φ as long as anisotropic stress is negligible, i. e. when T i j is small for i = j. Inspecting the energy-momentum tensor of a single particle, cf. Eq. (21), the off-diagonal components T i j are of order v 2 and thus only important for relativistic species. In our case, the only contribution comes from the neutrinos. Perturbation theory yields the following expression for the difference of the two potentials [16]: with the traceless components of the neutrino energymomentum tensor, The field Σ i j (x) is obtained from the neutrino particles by using the form of the one-particle energy-momentum tensor, Eq. (21). A straightforward calculation yields The third field, δϕ, is calculated from Eq. (20) in Fourier space. On the right-hand side, we need δT = δT α α . Again from Eq. (21), the one-particle contribution is Given gravitational potentials Ψ, Φ, and the cosmon perturbation δϕ, the acceleration of the effective particles is obtained from the corresponding equations of motion. In the case of neutrinos, this is Eq. (27). In the case of matter, the right-hand side vanishes, and the equation specializes to the usual geodesic equation.
For matter particles, we make the approximation v 2 ≪ 1, and simply find For the neutrinos, it is adequate to solve the equation of motion, Eq. (27), directly for the spatial components of the four-velocity u i = γv i (1 − Ψ)/a. Since this automatically respects u α u α = −1, the speed of light limit is robustly enforced. A differential equation for v itself could be unstable in this regard.
The gradients ∇Ψ, ∇Φ, ∇δϕ, occurring in the equations of motion are calculated in Fourier space and then transformed back to position space like the fields Ψ, Φ, δϕ themselves. This limits the accuracy on scales comparable to the cell size ∆x. In particular, forces between two particles located in the same cell are completely neglected. We will investigate the impact of this simplification in Sec. III E.
Fourier transformations are also used to estimate the power spectrum of perturbation variables, e. g. of the neutrino density contrast δ ν (x). The spectrum P ν follows from the definition, which we discretize on the grid.

C. Backreaction effects
The usual way to study cosmological dynamics is to perform a complete split between the background and the perturbation evolution. First, averaged energymomentum tensors and the unperturbed Friedmann metric are used to find the evolution of the background quantities. Second, the evolution of perturbations on the precalculated background is investigated.
This procedure neglects the modifications to the background evolution due to the presence of perturbations, the so-called backreaction. In the standard ΛCDM case, this is a reasonable approximation since gravitational backreaction is small [18]. In growing neutrino quintessence, in stark contrast, the backreaction plays a decisive role. For coupled quintessence models with couplings of super-gravitational strength, a significant impact of the non-linear structure formation on the background evolution is expected [9,10,18,19].
We now illustrate the importance of backreaction in our specific model and explain how we account for it in the simulation. The background dynamics of growing neutrino quintessence differ from ΛCDM because they include the joint evolution of the cosmonφ and of the neutrinos described by their energy-momentum tensor T α β . The equation describing this evolution is the modified Klein-Gordon equation, Eq. (19). Its right-hand side is given by the quantityT =T α α = −ρ ν + 3p ν . In the standard procedure, neglecting backreaction, one would estimateT from the evolution of the averaged energy-momentum tensor. Once the neutrinos have become non-relativistic, the trace is simply given by the average energy density,T = −ρ ν , whose evolution follows, cf. [7],ρ This equation tells us thatT will depend on the evolution of the averaged cosmonφ but does not take into account local variations δϕ.
These local perturbations, however, are very important when we regard the evolution of neutrino structures. A large part of the neutrinos is concentrated in non-linear structures where δϕ takes negative values. This leads to a systematic suppression of the masses M ν ∝ exp(−βδϕ). We shall see that Eq. (39) thus leads to a wrong estimate ofT .
The correct method to estimateT taking backreaction into account is to use the exact expression, which we obtain from the one-particle contributions, Eq. (36). Performing the spatial averaging, we get as a sum over the particles p, This expression underlines the influence of the perturbations. The quantityT is suppressed by two essential effects. First, we clearly see the dependence on the particles' actual masses determined by the local cosmon field. Second, once the fifth force has accelerated the particles to relativistic velocities, the Lorentz factors γ > 1 lead to an additional suppression. One could ask how large the inconsistency becomes if one neglected the backreaction effects. For this purpose, we run our simulation without modifying the background evolution due to the influence of perturbations. Technically, we use Eq. (39) for the background evolution. In order to quantify the inconsistency, we measure the average, Eq. (40), in every time step. The comparison between the two estimates of the same quantity −T = −T α α are shown in Fig. 1. In addition, we show the actual averageρ ν of the energy density. For early times, as long as perturbations are small, the two estimates agree well. But once non-linear structures form, the suppression of −T becomes obvious. This error carries over to the evolution ofφ by Eq. (19) and thereby to the whole subsequent evolution of the background.
It is thus important to incorporate the backreaction effects and to evolve the background and the perturbations in parallel. In every time step, we thus use the actual average, Eq. (40), for evolving the cosmon fieldφ by Eq. (19).
We further see the importance of relativistic corrections. Without these, we would have −T =ρ ν . Figure 1 shows that this approximation becomes invalid at late times. In order to determine the expansion rate H, we measure the actual averageρ ν in our simulation box at every time step, instead of using background equations. The rate H is then obtained from the Friedmann equation 3H 2 = iρ i a 2 .

D. Initial conditions
Our simulations use the results from linear perturbation theory to draw initial conditions. The choice of an initial redshift z i is motivated by two restrictions. First, as long as neutrinos are highly relativistic (for z 5), one cannot neglect the higher moments in their phasespace distribution function. Second, linear perturbation theory becomes invalid at z ≈ 2. We choose z i = 4 where the equation of state w ν has fallen to about 10 −2 while the neutrino perturbations are still small.
Although z i = 4 is a good choice for the neutrinos, matter perturbations have to be treated non-linearly much earlier. We thus start the N -body treatment of matter at z = 49. During the stage from z = 49 to z = 4, matter evolves according to gravity on a growing neutrino quintessence background.
Typically, N -body codes for CDM simulations use a displacement field to find growing mode initial conditions for the particles in the Zel'dovich approximation (see, e. g., [20,21]). In the more complicated case of growing neutrino quintessence, however, we prefer to use the results of the linear calculation to draw the initial conditions directly.

Initial random fields
Before we distribute neutrino particles, we have to draw random perturbation fields, namely, the density contrast δ ν (η i , x) and the peculiar velocity field v pec ν (η i , x). We focus on scalar perturbations, i. e., the peculiar velocity v pec ν is related to a scalar velocity perturbation v s via v pec ν (η, x) = ∇v s (η, x). Let us consider generic scalar perturbations δA(η, k), δB(η, k) in Fourier space. They are solutions of the perturbation evolution equations as given in [7]. These are linear and only depend on the absolute value k = |k|. Hence, they allow for basis solutions A k (η), B k (η). We consider the adiabatic mode and write The coefficient α(k) characterizes the concrete realization, whereas the basis solutions A k (η), B k (η) are merely attributes of the differential equations. The two-point statistics is then Assuming that, at an early cosmological time η prim , the perturbations were described by the primordial spectrum P prim (k), it follows We assume Gaussian primordial perturbations with a Harrison-Zel'dovich spectrum where A s denotes the scalar amplitude, n s the spectral index, and k pivot the pivot scale. A random realization of the fields δA(η i , x) and δB(η i , x) used as initial fields in the simulation is now obtained in three steps. First, α(k) is drawn as a Gaussian field with the statistics given in Eq. (43). Second, the linear code [7] is used to obtain the basis solution A k (η i ) and B k (η i ). Since the linear code works in synchronous gauge, a gauge transformation to the Newtonian gauge is necessary (cf., e. g., [16]). Finally, Eq. (41) determines the perturbation fields δA(η i , k) and δB(η i , k) in Fourier space, which can be transformed to real space.
The procedure above also applies to matter starting at z = 49. A consistent realization of the two species requires to use the same coefficients α(k) for both.

Discrete fields and particle distribution
In a numerical implementation, real and Fourier space are discretized. We will now explain how to draw initial conditions on a discrete grid {k i }. Therefore, we first have to derive a discretized version of Eq. (43). In particular, we work with a discrete Fourier transform (DFT) rather than the continuous version.
First of all, on the right-hand side, the Dirac delta function is approximated by the Kronecker delta, where we have used ∆k = 2π/L. On the left-hand side of Eq. (43), we replace the quantities α(k i ) by the discrete Fourier quantitiesα i related to real space by a DFT.
The relation between α(k i ) andα i becomes obvious when we discretize the Fourier integral. We write k i = (i 1 , i 2 , i 3 ) ∆k and x j = (j 1 , j 2 , j 3 ) ∆x with ∆x = L/n and the number n of cells per dimension, i. e. n 3 = N c . It follows where we have used the definition of the threedimensional DFT. We arrive at the discretized version of Eq. (43), expressing the variance of the Gaussian random number α i . Only approximately half of the N c numbersα i are independent. We have to respect the hermiticity condition ensuring that the transformed field α(x) takes on only real values. In this way, we obtain our initial random fields on the grid. Once the fields δ ν (x), δ m (x), v pec ν (x), and v pec m (x) are calculated, we have to represent them by a distribution of N -body particles. We freely choose a total number of N -body neutrino and matter particles, N ν and N m , respectively.
For the neutrinos, we have to account for thermal velocities v th , which, at z i = 4, are small but not negligible. They are described by a Fermi-Dirac distribution f 0 (v th ), cf., e. g., [16]. The total particle velocity is then v = v th + v pec . We approximate the phase-space distribution function of neutrinos by In this way, the averaged particle velocities at each cell reproduce the peculiar velocity field. The distribution function implies the number density At each grid point x, we distribute the rounded number ⌊n(x) (∆x) 3 ⌋ of particles and correct the error statistically by the addition of further particles. To reduce shot noise, particles are put at random positions within the pixel volume. The velocity v of a particle in a pixel x is drawn from the distribution f 0 (v − v pec ν (x)). In order to enforce the correct local average, we draw particles pairwise with opposite thermal velocities.
For matter, the procedure is analogous with the only difference that thermal velocities are negligible: We have checked this procedure by comparing the spectra estimated from the resulting particle distribution with those obtained by the linear code.

E. Numerical issues
The strategies presented in the previous sections take care, in principle, of all relevant effects in growing neutrino quintessence. As an important caveat, the method encounters numerical instabilities once very concentrated structures have formed. We will discuss this problem in the following section. Further, we will address the influence of the limited resolution on our results.

Instability for z 1
In Sec. III B, we collected several equations that allow us to estimate the potentials Ψ, Φ, and δϕ from the distribution of particles. These calculations, however, already require knowledge of the potentials, e. g. via the metric coefficients, cf. Eq. (32). We cannot resolve these mutual dependencies exactly.
As long as the time steps are chosen small enough, it is a good approximation to use, in these cases, the potentials of the previous step. This is the procedure we generally apply.
This can fail, however, if small errors accumulate in such a way that the evolution becomes significantly inaccurate or even unphysical. In fact, once the formation of neutrino structures is advanced, this problem occurs in the estimation of δϕ. This limits our capability to simulate growing neutrino quintessence for late cosmological times, z 1.
In more detail, Eq. (20) for δϕ contains a source term δT = δT α α , which, in turn, is sensitive to the local neutrino mass variations ∝ exp(−βδϕ). Formally, Eq. (20) can be written in the abstract form with a highly non-linear function f . Attempts to discretize the Laplacian and to apply Newton's method to directly solve this non-linear equation are numerically intractable due to the large number of cells.
Another approach is to consider the iteration If it converges to a fixed point δϕ * , this is clearly a solution of Eq. (51). As a starting point δϕ (0) , one may use the result of the previous time step.
The convergence of the sequence δϕ (n) , however, is not guaranteed. We apply this iterative scheme until z = 1, since, for later times, the sequence starts to diverge.
In order to get a rough idea how the behavior of the sequence depends on the physical conditions, we consider a very simplified configuration. We assume the idealized situation in which all the neutrinos in the volume V are concentrated in a dense structure of comoving size l small enough to be approximated by a Dirac delta shape. Let us fix the origin of the coordinates x at the position of the structure. The structure is assigned a total mass M (n) ∝ exp(−βδϕ (n) (0)) in the n-th iteration. For convenience, we take care only of the leading effects, i. e., we drop the metric perturbations and neglect V ,ϕϕ in Eq. (20).
The resulting expression for T reads where γ eff arises from the Lorentz factors of the neutrinos inside the structure. The iteration prescription, Eq. (52), in Fourier space yields, for k = 0, Transformed back to real space and evaluated at the structure position x = 0, where we have performed a cut-off of the Fourier integral at the scale of the structure size k max = π/l. This value determines the next mass estimate M (n+1) ∝ exp(−βδϕ (n+1) (0)). If we had started with the correct solution δϕ * , the mass M * would remain fixed in the iteration. Let us assume we had started with a small error, M (n) = M * (1 + ε (n) ). Then we find, at linear order in ε (n) , the size of the error in the next step, The iteration can only converge if the errors decrease, i. e., if the absolute value of the prefactor is smaller than unity.
Of course, the above discussion is only illustrative. Nevertheless, it clarifies that the strong coupling β and large concentrations M * /l of neutrinos in structures are working against the iteration.

Resolution
The equations of motion for the effective neutrino particles, Eq. (27), are not reproduced by mere two-body forces. Instead, they require the knowledge of the fields Ψ(x), Φ(x), and δϕ(x) on a grid. This differs significantly from pure CDM simulations, which obtain a high resolution by using two-body interactions for small-scale forces as in gadget-2 [22]. Since we do not employ an adaptive mesh, our resolution is limited by the constant size of a grid cell ∆x. While this clearly affects the accuracy of our method at small scales (comparable to ∆x), our approach still seems appropriate for studying the model on large scales at which neutrino structures first form.
As a check, we have run our simulation with a lower resolution, N c = 128 3 , corresponding to twice a bigger cell size ∆x ≈ 4 h −1 Mpc. Apart from N c , all other parameters as well as the initial random realizations have been chosen identical. In Fig. 2, we show the dimensionless neutrino spectrum, at redshift z = 1 obtained from the two simulation runs. As expected, the reduction of N c leads to a loss of power on small scales. On the other hand, we find satisfying agreement on large scales k 1/∆x ≈ 0.25 h/Mpc. These results indicate that we can obtain a robust picture of large-scale neutrino clustering from our simulations. When studying the properties of individual neutrino structures in Sec. IV B, the resolution of small scales will be more important. We will then again use the lowresolution simulation with N c = 128 3 to quantify the effect.

IV. NEUTRINO STRUCTURES
Using the strategies presented in Sec. III, we are able to study growing neutrino quintessence at the non-linear level. We first investigate qualitatively how the strong fifth force leads to the formation of large neutrino structures, Sec. IV A. Within these structures, the local mass variation becomes an important effect as we will see in Sec. IV B. The effects of neutrinos accelerated to relativistic velocities are presented in Sec. IV C. Since the cosmological neutrinos are not directly observable, the indirect influence of neutrino structures on cosmological observables is essential and will be the topic of Sec. V.
The model parameters used in our simulation are listed in Table II

A. Growth of structure
We run the simulation with the cosmological parameters listed in Table II and the simulation parameters of  Table I. Recalling that the cosmon-mediated fifth force felt by the neutrinos is a factor of about 2β 2 ≈ 5 × 10 3 stronger than gravity, it is expected that large non-linear neutrino structures will form. The numerical results we get in our simulation agree with this expectation. We can follow the precise time evolution of the structure formation process. We give an overview of the evolution from a = 0.25 to a = 0.5 in Fig. 3. The images show snapshots of the number density field n ν (x) of neutrinos in the simulation box (periodic boundary conditions) at intervals ∆a = 0.05. We will now discuss the main stages of neutrino structure formation.
Until a ≈ 0.3, no large non-linear structures have formed. The perturbations can still be described linearly, and our results reproduce the linear calculation [7]. We give an impression of these perturbations and their growth by showing a spherical, two-dimensional section of the three-dimensional field. Although mainly linear, the growth already is very fast. From a = 0.25 to a = 0.30, the perturbations have grown by a factor of about 3 to 4. The particle velocities, however, are still far from relativistic. The Newtonian limit would be applicable, and we will also see that backreaction effects are small.
From a ≈ 0.3 to 0.4, large-scale non-linear inhomogeneities are forming. Most cells in the simulation box are already empty of effective neutrino particles, which increasingly concentrate in large and dense filament-like regions. During this process, the velocities of effective neutrino particles increase and often reach relativistic values. The Newtonian limit breaks down, and the fully relativistic equation of motion becomes essential. The large density inhomogeneities lead to non-negligible lo-cal variations δϕ of the cosmon and in turn to important local mass variations. These modified masses together with the Lorentz factors γ > 1 cause backreaction effects according to Sec. III C. In fact, we shall find in Sec. V A that a ≈ 0.4, or z ≈ 1.5 is a characteristic time where the evolution of the cosmological background has started to feel the impact of large inhomogeneities.
After a ≈ 0.4, approximately spherical neutrino structures form (as already predicted in [17]), mainly at the intersections of the large filament-like non-linearities. In the course of the evolution, these structures become even increasingly concentrated and spherical. The spherical shapes can be regarded as a hint for virialization. In the simulation cells, the number densities of effective particles exceed the average value by factors up to 10 5 . The structures are, however, very different from matter structures. They produce local cosmon inhomogeneities δϕ resulting in locally varying neutrino masses. In the structures' centers, the neutrino mass is heavily suppressed, cf. Sec. IV B.
In Sec. III E, we explained that our strategy to calculate the cosmon perturbations δϕ becomes problematic once the concentration of neutrino structures becomes very large. For this reason, we cannot follow the subsequent cosmological evolution for a > 0.5, i. e. z < 1. We expect the overall picture to remain stable, i. e. a collection of neutrino structures with very large central overdensities. Due to the attractive fifth force, some of these structures are expected to merge if they are close enough to each other.

B. Individual structures
We have found several distinct neutrino structures in our simulation box at z = 1. In the following, we have a closer look at one of the most pronounced structures. In Fig. 4, we show a two-dimensional slice through our simulation volume located at the center of the structure. Since the structure is almost spherically symmetric, we can characterize its shape by a radial profile. We show the neutrino number density and the local neutrino mass as a function of the physical radial distance r ph from the center of the structure. A concentrated core is clearly visible with a neutrino number density contrast of 2×10 5 in the central cell.
A very interesting phenomenon of growing neutrino quintessence are the local neutrino mass variations. Since the cosmon-mediated fifth forth leads to a substantial neutrino clustering in regions of negative cosmon perturbations δϕ, the neutrino mass is typically strongly suppressed in overdense regions. This suppression is clearly visible inside our exemplary structure, cf. Fig. 4. We find that the neutrino mass at the innermost core of the structure is approximately one order of magnitude smaller than outside the structure. This result stresses the relevance of the local cosmon perturbations. It remains an open question how the structures evolve for z < 1. The idea that the neutrino mass profile inside the structures may become constant and independent from the evolution of the cosmological background has been raised in [10]. Our simulation method does not resolve the dynamics below the size of a cell, which is around 2 h −1 Mpc in comoving units. Hence, at small distances from the center, the quantitative results are affected by rather large errors. In order to estimate the uncertainties, we show the corresponding results obtained from a simulation with a twice as large cell size (dashed lines). The number density profile shows a similar shape but is significantly less concentrated. Still, the mass suppression reaches similar values.
We have seen that the neutrino structures are still becoming increasingly concentrated in the final stage of our evolution (cf. Fig. 3). This is most visible for the largest structures and less significant for smaller ones. As an example, we show how the number density profile of one of the smaller, but still very pronounced structures evolves from a = 0.45 to a = 0.5 in Fig. 5. We have verified that the structure is not undergoing merging processes in the considered time range. Apart from a moderate transfer of neutrinos from the outer regions of the structure to its inner core, no significant changes in the profile are visible.

C. Relativistic effects
The motion of the effective neutrino particles in our simulation respects the relativistic laws of motion, cf. Sec. II D, which remain valid for velocities approaching the speed of light. This is essential since the neutrinos feel a strong acceleration due to the fifth force. We observe in our simulations that the Newtonian laws of motion assuming small velocities become inadequate as soon as non-linear clustering begins, z ≈ 2 to 1.5, cf. Sec. IV A.
The importance of the relativistic treatment becomes manifest in the evolution of the average equation of state w ν =p ν /ρ ν . With the help of the methods explained in Sec. III B, we can calculate the average pressurep ν and energy densityρ ν of the neutrino fluid from the distribution of particles in our simulation. Performing the average over the contributions of single particles p to the densityρ ν = −T 0 0 , Eq. (30), yields The average pressurep ν then follows from the calculation ofT =T α α by Eq. (40),p ν = (ρ ν +T )/3. The evolution of the equation of state w ν is shown in Fig. 6. The figure shows a steep increase of w ν , grow- ing several orders of magnitude between z = 2.5 and z = 1, exceeding the value of w ν ≈ 0.1. The velocities of many neutrino particles in the box are close to the speed of light at this stage of the evolution. The figure also shows the evolution of w ν obtained by solving the homogeneous background equations, i. e. neglecting backreaction effects. If neutrinos were not accelerated by the cosmon-mediated fifth force, their equation of state would continue to decrease. The pronounced oscillatory features are attributed to the neutrino mass variation; the cosmon field oscillates around the minimum of its effective potential (cf. [6]).
We observe a striking discrepancy between the background calculation and the evolution of the actual average. The consequences for the evolution of the background cosmon fieldφ will be discussed in Sec. V A.
The presence of relativistic neutrinos is also interesting with regard to gravity. Our tools allow us to calculate the anisotropic stress induced by neutrinos, which is the source for the difference between the two gravitational potentials Ψ and Φ, Sec. III B. We find that the effect is most pronounced in the vicinity of neutrino structures and becomes negligible on very large scales. In Fig. 7, we show the field Φ − Ψ at z = 1 on the same twodimensional slice already used in Fig. 4. In the regions of neutrino overdensities, the field Φ − Ψ shows anisotropic patterns with amplitudes of the order 10 −5 . This is, compared to Φ and Ψ themselves, a 1% to 10% effect.

V. IMPACT ON DARK ENERGY AND MATTER
Once the dynamics of growing neutrino quintessence can be described reliably, a confrontation with observational data is the next step. At the current stage, our simulation method allows us to follow the cosmological evolution only until z ≈ 1, and we have not explored the parameter space of growing neutrino quintessence. It is thus not yet possible to provide constraints on the model parameters.
Nonetheless, it is insightful to discuss some remarkable effects linked to observations. This regards the expansion history of the Universe (depending on the perturbation evolution due to strong backreaction), Sec. V A, and the evolution of matter perturbations, Sec. V B.

A. Backreaction and quintessence
We argued in Sec. III C that two essential effects in the structure formation process lead to an important backreaction of the perturbations on the background evolution. These are the systematic suppression of neutrino masses within structures, m ν (ϕ(x)) < m ν (φ), and the relativistic Lorentz factors, γ > 1. Numerically, we have seen that these effects are indeed strong. The neutrino mass inside structures is suppressed by an order of magnitude, Fig. 4, and the neutrino equation of state w ν grows considerably at late times, Fig. 6. Consequently, the true averaged energy-momentum tensorT α β of the neutrino fluid significantly differs from the idealized calculation based on a homogeneous fluid of non-relativistic neutrinos.
The strong backreaction is not only interesting in its own right. In growing neutrino quintessence, the evolution of the dark energy is intimately connected, and hence sensitive, to the evolution of the neutrinos, cf. Sec. II. In particular, the cosmological event of the neutrinos becoming non-relativistic serves as a trigger stopping the further evolution of the cosmon and leading to an epoch of accelerated expansion. We have seen, however, that the fifth force accelerates neutrinos again to relativistic velocities in the course of the non-linear evolution. This reduces the neutrinos' capability of stopping the cosmon evolution. The onset of accelerated expansion is thus expected to shift to later times. More precisely, the source term in the background evolution of the cosmon in Eq. (19) is proportional to the traceT . This trace was shown to be a sum over particle contributions ∝ m ν /γ, Eq. (40). So, both the effect of relativistic particles, γ > 1, and the mass suppression go into the same direction.
For a quantitative investigation, we calculate the deceleration parameter q = −a ′′ a/a ′ 2 + 1. For the expansion of the Universe to accelerate, q has to take negative values. We show the evolution of q measured in our simulation, i. e., including backreaction according to Eq. (40), compared with the expectation based on the a-priori averaged background equations, cf. Eq. (39), in Fig. 8. We observe that the deceleration parameter, al- ready very close to zero without backreaction, is indeed still at a relatively large value q ≈ 0.2 at z = 1. Although our simulation ends at z = 1, we expect the phase of accelerated expansion to start later due to backreaction effects. This is because we can expect that the two main contributions, the mass suppression and the relativistic Lorentz factors, remain important in the subsequent evolution.
By means of Friedmann's equations, we can express the deceleration parameter q as a sum over contributions of the different species i, The cosmon contribution to the deceleration thus is which is a comparison of the kinetic energy and the potential energy (the potential is normalized to V (0) = 0.9 × 10 −120 M 4 P with the reduced Planck mass M P ). We can use V (ϕ) = V (φ) in linear approximation in δϕ. The difference between the actual average V (ϕ) and V (φ) is way below the percent level.
Due to Eq. (60), we can understand the evolution of the deceleration parameter, Fig. 8, by examining the background evolution of the cosmon. We show the equation of state w ϕ as well as the average valuesφ andφ ′ in Fig. 9.
In the lowermost plot, we see that once backreaction effects become important (at redshifts around z ≈ 1.5), the time derivativeφ ′ takes significantly larger values than it would be expected neglecting backreaction. This shows that, in fact, the neutrino fluid is less effective in stopping the cosmon evolution. This is also reflected in the plot ofφ, which continues to grow although, neglecting backreaction, a mere oscillation around a very slowly increasing value has set in at z 1.5. As a consequence, the equation of state w ϕ is, due to the backreaction, further away from the cosmological constant value w = −1.
We conclude that the study of non-linear structure formation in growing neutrino quintessence turns out to be crucial for the understanding of its correct background evolution as well. A realistic model with the correct fraction of dark energy today might need a readjustment of the model parameters.

B. Matter perturbations
The dynamics of matter perturbations in our simulation are only described by the gravitational force ∝ ∇Ψ and the Hubble damping. The effective matter particles are accelerated according to the Newtonian law of motion, Eq. (37), and we do not differentiate between baryons and dark matter. Matter is not coupled to the neutrinos or to the cosmon except via gravity.
The gravitational effects due to the neutrinos and the cosmon include their gravitational potential and their impact on the expansion history. In particular, the large neutrino structures investigated in Sec. IV induce gradients in the gravitational potential leading to an additional acceleration of matter particles. We thus expect that the amplitudes of matter perturbations are enhanced compared to the ΛCDM case [9,23]. In order to study this effect, we follow the evolution of the dimensionless matter spectrum, The concrete values of ∆ m (k) of course depend on the model parameters we have chosen for growing neutrino quintessence, Table II. Hence, in order to isolate the impact of non-linear neutrino structures on the growth of matter perturbations, we normalize ∆ m (k) by the case where matter only feels its own gravitational potential. This means, we have started another simulation run following the evolution of only matter particles in a universe with the expansion history of unperturbed growing neutrino quintessence. Since the expansion history is, for the chosen set of model parameters, similar to ΛCDM, we somewhat imprecisely label the corresponding amplitudes by ∆ ΛCDM m (k). The quotients ∆ m (k)/∆ ΛCDM m (k) at different redshifts are shown in Fig. 10. At early times, z > 2, the distribution of the neutrinos is still very smooth and matter fluctuations grow only due to their own gravitational potential. At later times, we observe an enhancement of matter fluctuations at large scales, where the neutrino fluid is clustering, cf. Fig. 3. At z = 1, the amplification exceeds 10% on the largest scales of our simulation, whereas it is still below the percent level at small scales.
By means of the continuity equation, the change δ ′ m of the matter density perturbation is linked to the velocity field v pec m . If the density perturbations show a steep increase, as we observe it on large scales, large peculiar velocities must be present.
We can measure the matter bulk flowv pec m in a subvolume V of our simulation defined as the average peculiar velocity,v Next, we introduce the variance of the bulk flow U 2 l on different comoving scales l as a measure of the expected bulk velocity in volumes of size l 3 . Technically, we divide the simulation box into n subvolumes V i of equal size l 3 = L 3 /n. We then obtain the average peculiar velocityv pec m,i within the boxes i and compute the variance according to We show the evolution of U l for two large scales in Fig. 11. Indeed, the enhancement of bulk flows is much more pronounced than the increase of density perturbations. This feature of growing neutrino quintessence has been studied in [9,23]. The increase of bulk flows reaches factors of about 1.5 to 2 at z = 1, and it is growing.
The large bulk flows are the consequence of large-scale gradients in the gravitational potential Ψ. Whereas matter perturbations are still linear on scales k 0.1 h/Mpc, the neutrino density field is highly inhomogeneous on these and larger scales. It comes as no surprise that the non-linear neutrino structures dominate the large-scale gravitational potential. In order to quantify this, we estimate the power spectrum P Ψ (z, k) of the gravitational potential. The dimensionless spectrum is a measure for the fluctuation variance of the cosmological gravitational potential on spatial volumes ≈ π 3 /k 3 . We compare the large-scale evolution of Ψ(z, k) (including neutrino structures) with Ψ ΛCDM (z, k) (neglecting neutrino structure formation) in Fig. 12. Since Φ = Ψ on large scales to a good approximation, cf. Sec. IV C, the corresponding plot for Φ(z, k) is visually identical. Once the non-linear effects of neutrino structure formation become important at z 2, the large-scale gravitational potential grows drastically compared to ΛCDM. Since the matter perturbations have only grown moderately under the influence of neutrino structures until z = 1, see Fig. 10, this effect is due to the neutrinoinduced gravitational potential. Consequently, Fig. 12  shows that the large-scale potential of the neutrino structures exceeds that of matter by about an order of magnitude. The large-scale gravitational potential Ψ(k), Eq. (64), at z = 1 reaches values between 10 −5 and 10 −4 in our simulation.
Whereas the large-scale gravitational potentials in ΛCDM are constant during matter domination and very slowly decay thereafter, growing neutrino quintessence shows a completely different behavior. During the formation of neutrino structures, the large-scale gravitational potentials become deeper. This is expected to have a drastic impact on the integrated Sachs-Wolfe effect, which is sensitive to the time derivative (Ψ + Φ) ′ . In this way, the evolution of the gravitational potentials might be decisive for scrutinizing the model. A continuing growth of the gravitational potentials for z < 1 would clearly be in conflict with observational results [24]. The evolution for z < 1, however, could be very different. The virialization of the neutrino structures might stop the growth of the neutrino-induced gravitational potential.

VI. SUMMARY
The cosmological dynamics of growing neutrino quintessence show a much higher complexity than the standard ΛCDM cosmology or models of uncoupled dark energy. This work has made a crucial step towards a consistent simulation of growing neutrino quintessence including all its relevant effects.
The standard methods of linear perturbation theory and (Newtonian) N -body simulations, while very successful in the ΛCDM case, rely on approximations in conflict with the effects of growing neutrino quintessence. Linear perturbation theory breaks down on all scales [7], and in Newtonian N -body simulations, the particles' ve-locities exceed the speed of light [9]. Furthermore, the local variation of the average neutrino mass due to perturbations in the dark energy scalar field could so far only be studied in idealized configurations [10]. None of the applied methods has shed light on the expected backreaction [19] of the highly non-linear neutrino perturbations on the evolution of the cosmological background.
We have thus decided to develop a new method coping with these challenges right from the start. This effort has led to an N -body based simulation, designed from scratch and adjusted to growing neutrino quintessence. In particular, we incorporate local variations in the dark energy scalar field and respect relativistic dynamics for the neutrino component. We run the simulation for an exemplary choice of model parameters. Our results show that all the effects that had to be neglected in previous approaches are indeed important.
The simulation confirms the formation of large-scale neutrino structures (already predicted in [17]), and we give a detailed description of their evolution (Sec. IV A). As a main result, we observe a significant neutrino mass suppression of an order of magnitude inside concentrated structures (Sec. IV B) agreeing with [10]. This is the consequence of significant inhomogeneities in the dark energy scalar field, the cosmon ϕ. Once these become too large, our method encounters numerical difficulties (Sec. III E), which is why we do not follow the cosmological evolution for redshifts z < 1.
We find that the velocities of neutrinos are accelerated to relativistic values during the process of structure formation. The average neutrino equation of state reaches values w ν ≈ 0.1 at z ≈ 1 (Sec. IV C). Further, we have shown that the relativistic motion induces a difference between the two gravitational potentials. In the vicinity of large neutrino structures, the difference Φ−Ψ typically amounts to ∼ 10 −5 .
Local mass variations and relativistic velocities of the neutrinos lead to strong backreaction effects (Sec. III C). We have demonstrated that the backreaction can notably modify the late-time expansion of the Universe (Sec. V A). It is likely that the onset of the accelerated expansion occurs later as compared to a homogeneous approximation.
The evolution of matter perturbations is affected by the neutrino-induced gravitational potential, which dominates on large scales (Sec. V B). As a consequence, large-scale bulk flows of matter are amplified by a factor close to two at z = 1. The density field, however, is much less affected. On large scales, we have found an effect of about 10%.
The current stage of the simulation method presented in this work is already promising, and the efforts to extend its range of applicability will continue. When the cosmological evolution can be followed until z = 0 for a collection of different model parameters, a confrontation of growing neutrino quintessence with observational data becomes possible.
So far, we are only able to identify several effects that might be interesting regarding observations. The largescale neutrino-induced gravitational potentials could be observable directly via gravitational lensing or indirectly via the resulting large-scale bulk flows of matter and the enhanced density power spectrum. The time evolution of the gravitational potential, on the other hand, leaves imprints in the integrated Sachs-Wolfe effect. In contrast to ΛCDM, we have observed increasing large-scale potentials in growing neutrino quintessence during the formation of neutrino structures.
The main motivation to search for alternatives to the cosmological constant scenario are the difficulties to understand the tiny value of Λ and the why now problem. Any competing model should avoid them in the first place. Growing neutrino quintessence offers a mechanism to solve the why now problem of dark energy. With a rich phenomenology beyond the standard scenario, the prospects are promising that the model eventually can be put to stringent tests once its non-linear evolution is completely understood.