Visualization of cosmological density fluctuations with phase space analysis: case study: Brans–Dicke theory

Cosmological perturbation theory is a powerful tool to understanding the large-scale structure of the Universe. However, the set of field equations describes the general linear perturbations for the cosmological models is highly nonlinear and coupled where no analytical solution can be found. It is only after some simplification and numerical computation that we obtain limited solutions. On way around is to employ the phase space analysis and investigate the asymptotic stability of the model. The advantage of using this approach is that the system of field equations become simpler to solve numerically. The algorithm also determines the stability of the system. Here, we apply this approach in Brans–Dicke cosmology and study the attractor solutions after fitting the model with the SNeIa observational data. As a result, the model confirms small fluctuations of energy density. The model also predicts current universe acceleration confirmed by observation and benefits from phase space analysis that the new dynamical variables are independent of the model initial conditions.


Introduction
Over the past few decades, extensive research has been conducted in order to understand the formation of large scale structures in the universe. More recently, a number of surveys such as the Sloan Digital Sky Survey (SDSS) and the Two-Degree Field (2dF) have have been used to measure the redshifts of millions of galaxies and plot the three-dimensional distribution of matter in the universe. It has been found that smoothing the distribution on the largest scales (larger than several hundred Mpc scales) displays a homogeneous distribution of matter, in accordance with the assumptions that a e-mail: salehi.a@lu.ac.ir b e-mail: h.farajollahi@unsw.edu.au the universe in the large scale is homogeneous and isotropic. At smaller scales, however, there exists density fluctuations; over-dense and under-dense regions like clusters and voids which increase at smaller smoothing scales. In terms of density fluctuation amplitude, the rms is over order unity on average ∼ 10 Mpc scales. Noting that in the past two decades, the cold dark matter (CDM) and dark energy ideas have also contributed into the standard theoretical models of galaxy formations.
A fundamental assumption in structure formation is that it started to grow from weak density fluctuations in the early universe, then, amplified by gravity and eventually turned into the current structure formation. Measurements of the cosmic microwave background (CMB) by the WMAP satellite in combination with the 2dFGRS confirm the assumptions and allow an accurate determination of the geometry and matter content of the Universe about 380000 years after the Big Bang [1][2][3][4].
The perturbation models can in general be divide into scalar, vector and tensor perturbations, based on how they transform under spatial rotations and local time shifts. Based on these, one then study the evolution of small perturbations of homogeneous metric and matter fields. Scalar perturbation models have been widely studied in fourth order gravity mainly in the metric formalism [5][6][7][8][9][10][11][12][13][14]. These models originally developed for General relativity by Bardeen [15] and then expanded in the 1 + 3 covariant approach [16][17][18][19][20][21][22][23][24][25][26][27][28][29]. Many attempts have been implemented to numerically compute the matter power spectrum and its numerical analysis for classes of modified gravity theories, rather solving the full set of field equations [30][31][32][33]. Modifications to the linear order Einstein equations are thus introduced in terms of general functions of scale and time. In general, all studies have proved that the perturbation's growth depends on the scale. While in General Relativity the evolution of dust matter is scale-independent [18], the numerical investigation of perturbation theories shows that the matter perturbation growth is scale -dependent [34,35].
Perturbations in scalar-tensor theories have been considered from different approaches, for example: the reconstruction problem was addressed in [36], the density contrast evolution was studied in [20][21][22][23][24][25][26][27][28][29]37,38] and the second order perturbations were presented in [39]. In all these approaches the gravitational instability is a fundamental assumption. That is, the early universe has been almost perfectly smooth, with the exception of tiny density deviations with respect to the global cosmic background density and the accompanying tiny velocity perturbations from the general Hubble expansion. The minor density deviations vary from location to location. At one place the density will be slightly higher than the average global density, while a few Mega parsecs further the density may have a slightly smaller value than the average. The observed fluctuations in the temperature of the CMBR are a reflection of these density perturbations, such that the primordial density perturbations have been in the order of 10 −5 . Noting that the origin of this density perturbation field has as yet not been fully understood. The most plausible theory is that the density perturbations are as a result of quantum fluctuations in the very early Universe which during the inflationary phase expanded to macroscopic proportions.
In this work, we consider scalar metric perturbation. We also assume that the matter is made up of a single real scalar field. In the case of General Theory of Relativity that the field equations are extremely coupled and nonlinear, the perturbation become more complicated in which analytical solutions of the governing field equations are impossible to solve or give little understanding about the behaviour of the system. Our novel approach to this problem is to utilise phase space and stability analysis as a fascinating technique to study the dynamics of the system. We apply the technique in the Brans-Dicke (BD) theory as an alternative to General Relativity.
Noting that the gravitational constant in Brans Dicke theory is non-constant. It dependents on the scalar field where its contribution to the Lagrangian density is through its own kinetic term. The evolution of the scalar field as a source term affects the matter distribution in the universe. This in turn makes the gravitational constant dependent of the matter field as well. This provides a manifestation of the Mach's principle as interpreted by Dicke. (see [40] for a recent review). The authors in [41][42][43] have extensively studies different formulation of Brans Dicke theories and their applications in modern cosmology The authors in [44] have already studied the evolution of the density contrast for Jordan-Fierz-Brans-Dicke theories in a Friedmann-Lemaitre-Robertson-Walker Universe. Following their work , we are going to investigate the density perturbation of the model based on the phase plane analysis. In particular, we focus on the late time stable solu-tions and give a qualitative description of the model. This has been done by using proper transformations in which the complicated equations become simplified. In general, there are two approach to examine the stability analysis of the model by simultaneously solving the dynamical system and best fitting the stability parameters with the observational data.
1-we find the critical points and their eigenvalues in terms of parameters of the model, find the stability conditions for each critical point and then constrain the parameters with observation.
2-We first constrain the model parameters and initial conditions with observation. Then find the critical points and eigenvalues for the best fitted parameters and investigate the stability of the model [41,45].
This approach is more appropriate for complex systems that finding the critical points and eigenvalues in terms of model parameters is difficult. This is the approach we followed in this paper Our work is organized as follows: In Sect. 2, Brans-Dicke theory has been extracted from Scalar-tensor theories of gravity, we summarize the background field equations of BD theories in both Jordan and Einstein frame. In Sect. 3 following [44], we implement the metric perturbation and drive the perturbed equations. In Sect. 4, we best fit our model with the observational data. In Sect. 6, phase space analysis of the model is discussed. The perturbed second order, coupled and nonlinear equations are transformed into simple ones by introducing new dynamical variables. The stability of the model and attractor solutions are discussed. Section is devoted to discuss the solution in more details. We further derive equation of state parameter in terms of the best fitted stability and model parameters and test the model against observational data. Section 7 is the summary.

Extracting Brans-Dicke theory from Scalar-tensor theories of gravity
Scalar-tensor theories are usually formulated in two different frameworks; the Jordan Frame (JF) or the Einstein Frame (EF). The former defines lengths and times as would be measured by ordinary laboratory apparatus so that all observables (time, redshift, among others) have their standard interpretation. Also, the metric is minimally coupled to matter and the scalar field is coupled with the Ricci curvature. Here, we start with the scalar-tensor action in Jordan Frame (JF) [36] as, where G * denotes the bare gravitational coupling constant, R the scalar curvature of g μν , and g its determinant. The variation of action (1) gives, In the Brans-Dicke theory F( ) = and Z ( ) = ω/ and 2Z F + 3(d F/d ) 2 = 2ω + 3. Thus the field equations become, While the above equations are in Jordan frame (JF), it is usually more transparent to analyze the equations and the mathematical consistency of the solutions in the Einstein frame (EF). This is achieved by conformal transformation of the metric and a redefinition of the scalar field. Let us call g * μν and ϕ the new metric and scalar field, and redefine Action (1) then takes the form where g * is the determinant of g * μν , g μν * its inverse, and R * its scalar curvature. Note that the first term looks like the action of general relativity, but that matter is now explicitly coupled to the scalar field ϕ through the conformal factor A 2 (ϕ). The field equations deriving from action (12) take the simple form where, quantities referring to the Einstein frame will always have an asterisk and α(ϕ) ≡ d ln A dϕ is the coupling strength of the scalar field to matter sources. In Brans-Dicke theory, using Eqs. (9)-(11) and by considering constant value for ω, it will be as, is the trace of the matter energy-momentum tensor T Using Eqs. (8), (10) and (16), it can be concluded that ρ * = A 4 ρ and P * = A 4 P, where in Brans-Dicke theory , A = e αϕ . Also, from Eqs. (8) and (10), ds 2 = A 2 ds 2 * , where ds 2 = −dt 2 + a(t) 2 dl 2 and ds 2 * = −dt 2 * + a(t * ) 2 dl 2 are the JF and EF metrics respectively. This relation indicates that, dt = A(ϕ)dt * and a(t) = A(ϕ)a * (t * ). Thus, for Brans-Dicke the field equations in EF become Here, dot denotes derivative respect to t * . We Also can derive the equations by replacing the physical time t * with the conformal time η * . Since, dη = dt a and dη * = dt * a * , thus,H = dlna dη = a H and H * = dlna * dη * = a * H * . We also have where prime denotes derivative respect to conformal time η.
The conformal time η is the same in both frames η * ≡ η. Thus, the field equations for Brans-Dicke theory in EF will be where, c 2 s = P * ρ * andρ * = 4π G * ρ * a 2 * .

Perturbations in BD theories
We begin with the field equations for the action (12). The field equations obtained from action, neglecting the potential, are, From now on, let us drop the notation (*) for EF quantities.
Up to the first order perturbation, the equations of the gravitational field in BD theories (25), the scalar field Eq. (14) and the conservation Eq. (15) respectively become Let us consider the scalar perturbations of the flat FLRW metric in the longitudinal gauge where ≡ (η, x) and ≡ (η, x) are the Bardeen's potentials [15,46] and η is the conformal time. We also assume adiabatic perturbations and barotropic fluids where the perturbed stress-energy tensor can be expressed as where ρ is the perturbed energy density of the cosmological fluid, δ is the density contrast and v is the potential for velocity perturbations. In the following, we also assume that both perturbed and unperturbed matter obey the same equation of state.
Using the perturbed metric (29), the perturbed adiabatic and barotropic stress-energy tensor (30) and also assuming that the background equations hold, considering the (i j) components which yields ∂( − ) = 0 and indicates that both scalar potentials are equal as in GR. We finally yield the following set of equations in Fourier space: 4ϕ where means derivative with respect to the conformal time.

Constraining the cosmological parameters
The Eqs. (31)-(36) describe the general linear perturbations for the model. The differential equations are second order, nonlinear and coupled in a large number of variables and are impossible to solve them analytically. Our purpose in this section is first to simplify these perturbed equations without losing the generality by converting the second order equations to the first order ones. This makes the equations simpler to solve numerically and and prepare them for phase space analysis in the next section. The phase space is in particular useful in visualizing the behavior of the dynamical systems and studying the stability of the system.
To simplify these equations, let us introduce a new set of variables as: One then can rewrite the Eqs. (31)-(36) for the autonomous equations of motion as where N = ln(a) such that d d N = 1 H d dη and the new parameters ζ 1 ..ζ 6 are given by: From Eqs. (31)-(36), we can obtain the above parameters in terms of the new variables. By dividing both sides of Eqs. (32) and (34) by H 2 and using new variables (37), the parameters ζ 2 and ζ 3 become Also, by dividing both sides of Eqs. (33), (36) and (35) by H and using new variables (37), the parameters ζ 4 , ζ 5 and ζ 6 change to Next, we derive the parameter ζ 1 in terms of new variables and the other parameters. To do this, we first take derivative of Eq. (33) with respect to η: Then, dividing both sides of (53) by H 2 results Note thatρ = 4π Ga 2 ρ, thusρ =ρ ρ ρ + 2Hρ. Using this relation and Eq. (24), the Eq. (54) can be rewritten as Combining the above equation with Eq. (48), parameter ζ 1 becomes By substituting (48)-(56) into the Eqs. (38)-(45), we obtain a complete set of equations that describe the behavior of the system in terms of new variables. Also by dividing both sides of (31) by H 2 and some calculation, one can obtain a constrain as that reduces the number of independent equations to seven. Noting that in the new system of equations the parameter k is replaced by χ 2 H. Next we are planning to best fit the parameters H 0 and α in the system of Eqs. (38)-(45) by using observational data. The SNe Ia observations measure the apparent magnitude m of a supernova that is related to the distance modulus μ and luminosity distance d L of the supernova by where M is the absolute magnitude of the supernova. In flat FRW cosmology, the luminosity distance is given by The (59) can be expressed by two new differential equations as We have a system of Eqs. (38)-(45) plus the above two differential equations for luminosity, (60)-(61) that has to be solves simultaneously. We assumed that the universe is in mater dominated state where c 2 s = 0. We constrain the model parameters H 0 and α and initial conditions of the dynamical variables and luminosity differential equations with the observational data, using Supernova Type Ia observational data.
We conduct the χ 2 goodness of fit that given by where the sum is over the distance modulus of 557 Type Ia supernovae (SNeIa) sample ( [47,48]). In Fig. 1-RHS, the green, red and blue regions show the confidence region for the constraint on the parameters α and H 0 at the 68.3 %, 95.4 % and 99.7 % confidence levels. A 3-dim of the likelihood is also shown in the Fig. 1-LHS. In the following section we use the best fitted value of H 0 and α to conduct the stability analysis in the phase space for our model.

Phase space analysis
In stability analysis, and in the context of autonomous dynamical systems, the fixed points ( or critical points) are the exact constant solutions that describe the asymptotic behavior of the system. This can be done by simultaneously solving Eqs. (38)- (45) such that χ 1 = 0, . . . , χ 8 = 0 . As mentioned, imposing the constraint reduces the number of independent equations to seven. Employing stability analysis procedure, we find thirteen critical points given in Table 1.
Here, we used the best fit value of α = 3.35 and H 0 . Substituting linear perturbations χ 1 → χ 1 + δχ 1 , . . . , χ 7 → Table 1 Critical points for the perturbed system Stability requires the real part of the eigenvalues to be negative. As can be seen, in our model all of the elements of Ev 10 is negative so this critical point P 10 is stable. Also the 3 th and 4 th elements of the eigenvalue Ev 10 is complex with negative real part. This indicates that the nature of this point is stable focus and the trajectories are spiralling in towards the attractor.
In Fig. 2, the dynamical variables χ 1 , ..., χ 8 for the best fitted value of parameters α = 3.35 and h 0 = 0705 are plotted in 3-dim phase space. As can be seen, by small perturbation, the system start to oscillate and asymptotically approaches the stable state of attractor.
In Fig. 3, the 1-dim plots of the dynamical variables are given with respect to the redshift z for both best fitted and arbitrary initial conditions.The graphs show that the solutions are independent of the initial conditions. This is because, in the infinite time limit, the system approaches the attractor, independent of the initial conditions.
In Fig. 4, the solutions are sketched only for the best fitted model parameters and initial conditions. It shows that the system oscillates while approaching the attractor state ( the critical point P 10 ) in equilibrium. As can be seen in the graphs, over the course of time the oscillation decays, becomes less visible and accordingly the trajectory of the dynamical system approaches the focus closer and closer, as if spiralling into the critical point. Note that the attractor behavior reveals more as the system approaches the critical point. Asymptotically, the small perturbation eventually settles down and system approaches the attractor.

Discussion
In the last section by applying the phase space approach and using variables, {χ 1 , χ 2 , . . . , χ 8 }, we converted the set of coupled second order differential equations to first order equations and studied the stability of the model. The new variables are given in terms of original variables (H, , ϕ, δ, υ,ρ, k) and their first and second derivatives. As discussed, in stability analysis, we study the dynamics of these new variables instead of the original ones. In fact, one can simply show that the system whose oscillations determined by the initial conditions are: On the other hand, in phase space analysis, we plotted the new variables with respect to the redshift (Figs. 2, 3, 4). Noting that the trajectories of these variables are for the best fitted values of α and initial conditions. As explained, the trajectories are independent of the initial conditions as expected in phase space approach. By perturbing these variables, the system begins oscillating while the trajectories are departed from each other for different initial conditions. Eventually after a long time, where the system converges to the attractor at equilibrium, all the trajectories emerges to a unique final stable state.
Note that non of the variables (χ 1 , χ 2 , ..., χ 8 ) are observable. So we need to test our model against cosmological parameters that are directly or indirectly observable. For this purpose we sketch the equation of state parameter, w e f f that in matter dominated universe vanishes and its current value in an accelerating universe is about −0.82. Fig. 5 shows the EoS parameter for the best fitted values of α and initial conditions. As can be seen the parameter vanished in the matter dominated era and is equal to −0.82 today which is within the range of observation.The graph also shows that, by perturbing the system, the parameter socialites and in near future approaches −1, i.e. the universe will be in C DM phase.

Conclusion
This paper is devoted to study the conventional metric perturbation in the Brans-Dicke thoery using the phase space algorithm. in practice, we integrated two different perturbation model; the standard metric perturbation and the perturbation based on phase space analysis of the newly constructed variables from the original geometric, matter and other fields. The conventional scalar perturbation in the metric introduces the scalar field, , that together with the brans-Dicke scalar field ϕ, the matter density ρ and Hubble parameter constitute the model that can potentially explain the large scale structure of the universe.  In phase space formalism, the newly constructed variables are used to manifest the attractor characteristic of the model. It is shown that the scalar perturbation of , δ, υ and ϕ, strongly depend on the initial conditions while the newly constructed variables are independent of them. By infinitesimal perturbation of these new variables, the solutions approach the same state corresponding to the attractor state. Noting that for those initial conditions largely different from the one that create the attractor point, the solutions will be represented as an oscillation with large amplitude while for the initial conditions close to the one that create attractor point, the oscillation is with small amplitude. In standard metric perturbation, the complexity and nonlinearity of equations in addition to exceeding the number of unknown parameter's, and increasing the sensitivity of the solutions to the initial conditions make it difficult to to interpret the large scale structure of the universe. However, its application in phase space creates a novel mechanism that simplifies both the complexity of the problem and system analysis. here, we also study the stability of the system and visualise the behaviour of the dynamical variable. More importantly, by fitting the parameters of the model, we supported our result with observational data. Though, non of the original variables in metric perturbation or newly constructed variables in phase space are directly or indirectly observable, the EoS parameter resulted from model fitting is in good match with observational data. Also It is important to note that while phase space analysis provides a qualitative description of the model, the fitting algorithm create some quantitative variables compatible with observation.
Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: We have used the latest Union2 compilation related by the Supernova Cosmology Project(SCP) collaboration recently [48]].