Copyright 1991-2015 John Bruce Davies

This research was awarded an Honorable Mention in the 1991 Gravity Research Foundation essay competition.


By examining a fluid dynamic matter-dominated Universe just after the separation of matter from radiation, we show that the governing non-linear equations have solutions that produce clumped distributions of matter. In the limit of zero flow velocities, the system of equations reduces to a single non-linear Klein-Gordon wave equation for the gravitational potential. Known analytical solutions of these equations yield clumped distributions of matter. The full set of nonlinear equations is then numerically solved for an arbitrary region of the Universe with initial conditions being those at the time of separation. Clumping of matter separated by low-density voids is obtained through self-organization of the initial flow perturbations. Predicted void size is of the correct order as expected from scaling arguments on observed large scale structures. Computed flow velocities correlate with the galactic streaming motions observed in clusters.


At large scales, much of the luminous matter in the Universe is clumped closely together in galactic superclusters with large low-density voids between the clusters.

Such matter distributions are fossil relics from the early Universe. The origin of this clumping has become clearer in the last two decades. In 1992, the COBE satellite observed temperature variations in the Cosmological Background Radiation (CMB) at the time of decoupling from matter at a Universe age of about 380,000 years. These variations were mapped in detail by WMAP in 2003 which also gave an age of the universe at 13.7 billion years. In 2009, the Planck satellite is expected to increase the level of detail sufficiently to clarify whether there are any specific patterns of hot and cold in the early universe and what is the specific statistical distribution of such variations.

Increasingly over the last few years, observational data has shown that the galaxies are not uniformly distributed throughout space-time. At all scales out to the horizon, the clumps, sheets and filaments of luminous matter are seen to be separated by large voids. Bahcall (1988) has shown that the large galactic supercluster structures observed today were mainly formed in the early Universe because the streaming velocities are considered too small to have affected this fossilized matter distribution.

Subsequent galaxy formation is usually modeled by linear perturbations of the fluid dynamical equations governing the matter fields. These density perturbation models are however unable to explain the observed clumped distribution of luminous matter.

Because the full set of governing equations are highly nonlinear their solutions can be expected to demonstrate nonlinear behaviour. Because matter density at the time of decoupling is well below relativistic constraints, we use the Laplace law relating density and gravitational potential in the Newtonian limit. We extend this in order to incorporate time-dependence of the potential, as per Chapter 2. This and the relevant polytropic equation of state together with the conservation laws and appropriate boundary and initial conditions govern the matter-dominated era. Inhomogeneous mass distributions are shown, analytically and numerically, to be solutions of this nonlinear system.

This set of nonlinear equations does not have obvious exact analytic solutions. Assumption of a spatially constant density and Hubble velocity relation gives a power-law dependence in time of Universe expansion, Zeldovitch (1977). Our approximation in order to obtain a single governing equation for analysis is to ignore the fluid dynamical velocity fields as per Bahcall's argument. Assuming a polytropic equation of state and combining the reduced equations results in a single equation governing gravitational potential. This is the well-known Klein-Gordon equation that for most equations of state is non-linear. Its known solutions possess amplitude dependence of the dispersion relation, peaking and breaking of waves leading to possible shocks, unstable waves growing in time, waves that split in two, and solitary waves of constant shape. Physical phenomena governed by similar non-linear Klein-Gordon equations, Whitham (1974), include grain boundaries in crystal systems and cracks in media around explosions, c.f. Davies (1980), both of which have defect distributions with their clustering separated by large voids.

In order to solve the full nonlinear set, numerical simulation is used to model a subregion of the expanding Universe. Standard methods are used to integrate forward in time from the instant at which matter decoupled from radiation. Random velocities due to thermal mixing are input into the system of equations with subsequent self-organisation of the matter into clusters. It is found that positive density anomalies tend to clump together and become separated by large low density zones. Computed sizes of low density zones are of the correct order, after scaling with Universe expansion, as those of the voids observed in the present galactic cluster distribution. Predicted streaming velocities are of the same order as those observed in galactic motions.

Governing Equations

The continuity equations are based on the Newtonian limit of the fields at particular points in space and time. Conservation of mass is expressed as:

(1) δρ/δt + Δ.(ρu) = 0

where ρ is density and Δ is the divergence operator and u is the vector of velocity. Conservation of momentum is similarly expressed as:

(2) ρ[δu/δt + (u.Δ)u] = -ρΔΨ + ΔP

where Δ is the gradient operator, Ψ is gravitational potential and P is the pressure. Conservation of energy gives for the temperature of the gas:

(3) ρ[δ/δt+(u.Δ)](CvT) = Δ.(K ΔT) - P Δ.u

where T is temperature, Cv is specific heat at constant volume, K is thermal conductivity, and viscous dissipation is ignored.

In order to incorporate expansion and its effect on the gravitational potential, the explicit time dependence of the gravitational field must be incorporated. It is thus necessary to use the non-relativistic limit of Einstein's field equation which corresponds to the time-dependent Newton's Law, c.f. Weinberg (1972), p. 156, i.e:

(4) [-1/c2.δ2/δt2 + Δ2]Ψ = -4πGρ

Davies (1988) has shown, in a more general derivation, that this relation is a consequence of the linearization of the Bianchi Identities in the most general space-time.

Klein-Gordon Solutions

In order to obtain an approximate solution, the limit of zero velocity for the flows is taken. This results in equation (2) becoming that of hydrostatic equilibrium:

(5) ρΔΨ = ΔP

The applicable equation of state that relates pressure and density in this system is taken to be polytropic:

(6) P = kρn

where n is the usual adiabatic gas exponent and k the corresponding constant. On differentiating the equation of state and substituting into the hydrostatic equation, we obtain:

(7) nkρ(n-2)Δρ = ΔΨ

which, on integration, yields:

(8) nk/(n-1). ρ(n-1) = Ψ + C

where the constant of integration, C, is taken to be zero for the reference state. Substituting equation (8) into (4) gives the fundamental time-dependent equation governing the gravitational potential:

(9) [-1/c2.δ2/δt2 + Δ2]Ψ = -4πG[(n-1)/nk.Ψ]1/(n-1)

This is the standard form of a Klein-Gordon equation in 3 spatial dimensions. Note that the left-hand side is linear while the right-hand side changes for the various equations of state. For monatomic hydrogen, n = 5/3, and the right-hand side is to the power 1.5.

Whitham (1974) shows that only in the linear case is the dispersion relation independent of the amplitude of the wave. Thus it can be expected that amplitude dependence of the dispersion relation will apply to these nonlinear waves in the early Universe. The near-linear approximate solution of the non-linear Klein-Gordon equation is obtained by expanding the r.h.s. in a Taylor series. It consists of the primary linear wave that travels at the velocity of light and corresponds to expansion of the Universe. Superimposed on this primary expansion are secondary waves which correspond to material motions and flows, Whitham (1974). Analytic solutions are known for certain nonlinear functions, Matsuno (1987), but, in general, Taylor expansions of similar nonlinear functions have yielded complex waves that exhibit inhomogeneous distributions of matter with clumping.

Numerical Simulations

In order to obtain total solutions, finite flow velocities are incorporated and the complete original set of nonlinear equations is solved numerically. This allows the important feedback effects of the advective terms to be incorporated and also prediction of void sizes and streaming velocities.

A finite arbitrary subregion of spacetime is modeled by standard difference techniques applied to the above primitive equations, (1 - 4). Initial conditions are taken as those at the time when radiation and matter decouple, c.f. Weinberg (1972), Table 15.4. However, the equations are valid in earlier times where mass gravity effects are due to element nuclei that are being formed. Initial conditions are found to be those with the gravitational potential and its time derivative set to zero. Initial densities are allowed to be variable with expansion-scaled dependence on the present density range expected for the Universe.

As usual, Zeldovitch (1977), pressure effects on the continuity laws are small compared with gravitational effects in this early Universe. However, the stochastic effects of thermal mixing are incorporated into the solution by introducing a random component into the velocities. This is a standard approach in analyses of the self-organising behaviour of nonlinear systems, c.f. Nicolis and Prigogine (1989). It is found that even small random velocities are sufficient to produce self-organization of the material flows into clumps of higher density separated by large low-density voids.

Figure 1 shows a time sequence, with intervals of 10,000 years, of density solutions across a typical cross-section using a random velocity disturbance of uniform distribution. The size of each square is one tenth the radius of the Universe at separation time. As time progresses, it is readily seen that clumping increases with fewer but larger voids (blue) separating the higher density (red) clumps. Figure 2 shows the various matter fields at a particular time after initiation. It is seen that all fields tend to clump with density, temperature and pressure having similar patterns. However, even though the velocity fields also clump, there appears to be little spatial correlation between them and with the other thermodynamic fields.


The main result from this nonlinear treatment of the early material Universe is that clumping of matter occurs rapidly within a time span on the order of a hundred thousand years. This results in a complex clumped distribution of matter and its associated fields and flows. The inhomogeneities so formed will remain as a fossilized relic throughout the subsequent expansion as the streaming velocities are too low to alter appreciably these clumped patterns over the age of the Universe, Bahcall (1988). An unexpected result of the numerical technique is that stable solutions are only possible with the initial conditions of a zero gravitational potential and a zero time derivative of the potential. It is also found that the initial density used in the numerical solution has little effect on the clumping phenomena.

Both the simplified Klein-Gordon solutions and the more complete numerical treatment yield clumping separated by low-density voids. Within the context of the restrictions on the numerical solution, it is seen that voids have dimensions of the order of 1/10 to 1/100 of the radius of the Universe at the time of separation. Observed voids in the present-day Universe have dimensions of up to a few hundred million lightyears which is of the order of 1/40 of the present observable Universe radius. Thus the observed relic structures are of the order predicted by this nonlinear system.

The temperature at the time of separation is of order a few thousand degrees. Random thermal velocities at this time are assumed to be less than 10-4 that of light. During the clumping phase, the flows aggregate in regions of both inward and outward flux velocities which can attain values up to the order of 1/1000 of that of light. Streaming velocities of the galaxies within clusters is observed to be of this order and therefore, assuming constancy of flow velocities during Universe expansion, correlate with the flux velocities predicted for the relic structure.

Evidence for large structures and their properties in the Universe, including giant congregations of galaxies and clusters, Lynden-Bell et al. (1988), the "Great Wall", a huge sheet of galaxies, Geller and Huchra (1990), possible periodicity of the galaxy distribution on a 100-megaparsec scale, Broadhurst et al. (1990), and infrared surveys of galaxy distibution, Saunders et al. (1991). Devlin et al. (2009) have shown that most of the far-infrared background light comes from massive galaxies at 1<z<4, just after the decoupling of matter from radiation. These occur in clusters and superclusters and imply that the active growth stage of most galaxies that are seen today is well behind them. As Small comments in the same issue of Science: ‘almost all of the sources detected are extremely luminous and very distant starbursts, potentially representing the formation of massive galaxies in a single burst of star formation’.

As Lindley (1991) has remarked, the disparate nature of the various large structures means that one should not hope to predict them in any specific way: ’it is sufficient that one's theory harbours some reasonable probability that such a structure might arise in the volume of the observable Universe’. This nonlinear treatment of the early Universe satisfies these conditions for producing large-scale structures with the self-organisation into clumps and voids as an inherent property of the governing nonlinear gravitational and fluid-mechanical equations..


Bahcall, N.A., 1988, Ann. Rev. Astro. Astrophys., 26, 631.

Broadhurst, T.J. et al., 1990, Nature, 343, 726.

Davies, J.B., 1980, ‘The Geometry of Faulting’,

E.O.S., Trans. Am. Geop. Un., 61, 1051.

Davies, J.B., 1988, ‘New Curvature-Torsion Relations through Decomposition of the Bianchi Identities’, Foundation of Physics, 18, 5, 563.

Davies, J.B., 1991, “Nonlinear Clumping in the Early Universe”, Honorable Mention, Gravity Research Foundation Essay Competition.

Geller, M. and Huchra, J., 1990, Science, 246, 897.

Lilly, S.J., 1988, Ap. J., 333, 161.

Lindley, D., 1991, Nature, 349, 14.

Lynden-Bell, D. et al., 1988, Astrophys. J.,326, 19.

Matsuno, Y., 1987, J. Math. Phys., 28, 2317.

Nicolis, N. and Prigogine, I., 1989, ‘Exploring Complexity’, Freeman Press.

Saunders, W. et al., 1991, Nature, 349, 32.

Weinberg, S., 1972, ‘Gravitation and Cosmology’, John Wiley & Sons, Inc., N.Y.

Whitham, G.B., 1974, ‘Linear and Non-Linear Waves’, Wiley Press.

Zeldovitch, Ya. B., 1977, Ann. Rev. Fluid Mech., 9, 215.