{'fraction_non_alphanumeric': 0.029406864378540485, 'fraction_numerical': 0.01016327890703099, 'mean_word_length': 4.65975168945466, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tessellating the Universe: the Zel'dovich and Adhesion tiling of space

The adhesion approximation is a simple analytical model suggested for explanation of the major geometrical features of the observed structure in the galaxy distribution on scales from 1 to (a few)x100/h Mpc. It is based on Burgers' equation and therefore allows analysis in considerable detail. A particular version of the model that assumes the infinitesimal viscosity naturally results in irregular tessellation of the universe. Generic elements of the tessellation: vertices, edges, faces and three-dimensional tiles can be associated with astronomical objects of different kinds: clusters, superclusters and voids of galaxies. Point-like vertices contain the most of the mass and one-dimensional edges (filaments) are the second massive elements. The least massive are the two-dimensional faces and tiles (voids). The evolution of the large-scale structure can be viewed as a continuous process that transports mass predominantly from the high- to low-dimensional elements of the tessellation. For instance, the mass from the cells flows into faces, edges and vertices, in turn the mass from faces flows into edges and vertices, etc. At the same time, the elements of the tessellation themselves are in continuous motion resulting in mergers of some vertices, growth of some tiles and shrinking and disappearance of the others as well as other metamorphoses.


Introduction
We live in an expanding universe theoretically predicted by Russian mathematician A. Friedman in 1922 [14] and independently discovered by American astronomer E. Hubble in 1929 [15]. It means that in the perfectly homogeneous universe the relative velocity, v, between any two particles or objects would be described by the Hubble law v = H(t)r, where r is the distance between the objects, the velocity is directed along the straight line connecting the objects. In the real inhomogeneous universe this law can be applied only at sufficiently large distances. A positive function H(t) is called the Hubble parameter, it characterizes the rate of expansion of the universe. Its present value H 0 = H(t 0 ) ≈ 71 km/s/Mpc (Mpc stands for megaparsec, 1 Mpc = 10 6 pc ≈ 3,260,000 light years). Until the end of the millennium most cosmologists believed that the expansion of the universe was faster in the past and has been monotonically slowing down with time (i.e.,Ḣ < 0 at t > 0) due to the pulling effect of gravity. However the study of extremely distant supernovae stars revealed that the expansion of the universe recently (by cosmological standards) started to accelerate [31]. The acceleration of the universe can be easily explained if the most of mass in the universe is in the form of dark energy -a hypothetical form of energy uniformly filling all the space and experiencing a repulsive gravitational force. Although the physical nature of dark energy remains a mystery it is currently the most popular way of explaining the acceleration of the expansion of the universe. According to current measurements of the cosmological parameters dark energy accounts for 74% of the total mass-energy, 22% are due to dark matter -another puzzling component, while the ordinary matter often dubbed as baryonic matter makes only about 4% of the total mass. Dark matter is believed to consist of weakly interacting massive particles that interact with ordinary matter and themselves via gravity only for most intents and purposes. Although the exact kind of particles has been neither identified nor detected in physical experiments the hypothesis of dark matter is well founded. The gravitational interaction of dark matter is similar to that of ordinary matter. Dark matter plays a crucial role in the formation of the structure in the universe due to its dominance in density and therefore in generating attractive gravitational force. The universe is also permeated by the Cosmic Microwave Background (CMB) radiation that brings unique information about tiny initial perturbations of the uniformity of the universe that served as the seeds for the presently observed structure in the universe. However CMB radiation makes less than 0.01% of the total mass of the universe at present and has played a very limited role in the structure formation after the universe was roughly one thousandth of its present size. It has been also established that the geometry of the universe is very close to flat. The present universe is highly inhomogeneous on scales up to several hundred of Mpc as shown in Fig. 1. However, this figure does not show a true distribution of galaxies in physical space because the line of sight distances have been derived from the receding velocities via the Hubble law (Eq. 1). In the inhomogeneous universe the actual relation between the line of sight distance, r, and recession velocity, v obs , is given by the equation that involves an additional term, v p , called the peculiar velocity. The peculiar velocities are due to growth of density inhomogeneities in the universe that obviously requires a transport of mass from one place to another and therefore demands additional (with respect to uniform expansion) velocities. The recession velocities of galaxies, v obs , can be measured by making use of the Doppler effect, while the peculiar velocities cannot be easily measured. Therefore, the line of sight distance derived from the observed velocity r est = v obs /H 0 = r + v p /H 0 gives only a crude estimate of the true distances of galaxies. The density distribution in physical space is better illustrated (albeit only in statistical sense) by the N-body simulations of the structure formation in the currently popular cosmological model referred to as the Lambda Cold Dark Matter (ΛCDM) model shown in Fig. 2. In cosmology the term large-scale structure is commonly referred to a distribution of galaxies on the scales roughly from 1 to a few hundred Mpc. Galaxies can not probe smaller scales because of discreteness, and on larger scales the galaxy distribution is monotonically approaching homogeneity with the growth of the scale. The redshift surveys reveal spectacular abundance of structures often described as The scale is shown in the velocity units (km/s), the distance in Mpc can be obtained from the Hubble law (Eq. 1), e.g., 10,000 km/s correspond to 141 Mpc. Each point represents a galaxy. Adopted from [11]. filamentary, network, or bubble structure [10,8]. The origin of the large-scale structure is one of the most important problems in modern cosmology. Many fundamental issues in physics, cosmology and astronomy ranging from speculations on the physical nature of dark matter, to the measurement of angular anisotropy of the microwave background radiation as well as establishing the epoch of galaxy formation join together here, see e.g., [30,34].
Modern theory explains the formation of the structure in the universe as a result of growth of primeval tiny density fluctuations in the process known as gravitational instability. Primeval perturbations are assumed to arise as vacuum fluctuations during the very early stage when the universe was expanding exponentially (inflationary universe), see e.g., [28]. Afterwards, the density perturbations had a long evolution before they become galaxies, clusters of galaxies, superclusters and voids of galaxies. The formation of galaxies is a very difficult problem on its own. Many complex physical processes like star formation and supernova explosions are very important for understanding the galaxy formation. We shall discuss the mass distribution only assuming that galaxies are fairly good (though not perfect) tracers of mass on large scales.
As long as the density perturbations are small by amplitude their evolution is described by the linear theory of gravitational instability (see e.g. [30,34]). The linear theory is quite simple and therefore is well understood. In particular, it predicts the rate of growth of the perturbations, D(t): δ ≡ δρ/ρ ∝ D(t), where δρ = ρ −ρ and < δ 2 > 1/2 ≪ 1. However, in commonly used Eulerian form it is not very useful at the nonlinear stage when the amplitude of the density perturbations becomes large and the observable structures (sheets, filaments, and clusters of galaxies) form. The analysis of the evolution of perturbations at the nonlinear stage < δ 2 > 1/2 ≥ 1 becomes quite difficult. The most straightforward way to address the complexity of nonlinear evolution is to carry out three-dimensional N-body simulations (Fig.2). Usually in the simulations of this type the medium is assumed to consist of collisionless particles, in agreement with the hypothesis of dark matter in the form of weakly interacting particles. The trajectory of each particle in the simulation is numerically integrated in the gravitational field generated by all the rest particles. Boundary conditions are commonly assumed to be periodic that imitate the infinite universe.
Here we outline another approach to the problem of the large-scale structure in the universe. We present two approximate analytic solutions to the set of PDEs describing the growth of density inhomogeneities in an expanding universe. One suggested by Zel'dovich in 1970 [38] is known as the Zel'dovich approximation and the other is based on Burgers' equation [5] and is often referred to as the adhesion model [12,13] (see also [34]). A particular version of the adhesion model naturally describes the structure in the universe as an irregular tessellation or quazi-Voronoi tessellation [27]. The generic elements of the tessellation in 3D (vertices, edges, faces, and tiles) can be associated with the observed structures in the three-dimensional distribution of galaxies (clusters, superclusters and voids of galaxies).
Both the N-body simulations and approximations require initial conditions. In the linear regime the density fluctuations are assumed to be a realization of a Gaussian random field specified by the power spectrum and the amplitude. The current measurements of the angular fluctuations in the temperature of the microwave background radiation by WMAP (Wilkinson Microwave Anisotropy Probe) are in excellent agreement with the assumption of Gaussianity. [22]. The amplitude of the temperature fluctuations suggests that the scale where the density fluctuations have recently reached nonlinearity is about 3 -6 Mpc depending on exact definition of the scale of nonlinearity. This is also in a good agreement with the observations of the large-scale distribution of galaxies. The shape of the power spectrum of the initial (i.e., linear) density fluctuations is fully determined by the parameters of the cosmological model. At large scales (small k) it behaves as P (k) ∝ k while at very small scales it is approximately P (k) ∝ k −3 with smooth transition between these limits. It is approximately P (k) ∝ k −1.5 at the present scale of nonlinearity, see the discussion bellow and Fig.5.
Here we discuss the models that deal only with the dominating component of the mass able to cluster i.e., the dark matter. The luminous objects in the universe are made of baryons, therefore it is important to understand the dynamics of the baryonic component of mass as well. However, this problem is much more complicated and is beyond the scope of this review. We refer the interested reader to paper [17] that offers an interesting generalization of the models described here and addresses the dynamics of baryons.
The rest of the paper is organized as follows. The equation describing the evolution of density inhomogeneities are present in Sec. 2. Sections. 3 and 4 briefly describe the Zel'dovich approximation and the adhesion model respectively. Section 5 provides a short summary.

Basic Equations
The evolution of density inhomogeneities can be described by a system of three partial differential equations comprising the continuity, Euler, and the Poisson equation (see e.g. [30,34]). In order to exclude the uniform expansion of the universe comoving coordinates, x, and peculiar velocity, u p are commonly used. They are defined by the following relations A monotonically growing function a(t) called the scale factor describes the uniform expansion of the universe, and H(t) = d ln a/dt is the Hubble parameter introduced by Eq.1. The scale factor is completely determined by the parameters of the cosmological model. For instance, a(t) ∝ t 2/3 in the Einstein -de Sitter universe where Λ = 0, and Ω = 1, while in more realistic ΛCDM model it is a complicated but known function of time [6]. In the course of evolution neither peculiar velocities no gravitational potential reach the relativistic values, therefore the use of classical mechanics and Newton's gravity law is perfectly appropriate. In terms of the comoving coordinates and peculiar velocities the equations describing the gravitational instability in the expanding universe are as follows: the continuity equation the Euler equation and the Poisson equation where ρ and ρ are respectively the density and mean mass density; φ is the gravitational potential due to the inhomogeneities of density; G is the gravitational constant. The pressure is neglected since we study the medium interacting only gravitationally. Additional terms on the right hand side of the continuity and Euler equations (i.e., −3Hρ and −Hu p respectively) are due to the expansion of the universe and the factor 1/a is due to differentiation with respect to the comoving coordinates x: ∇ ≡ ∂/∂x i ≡ a · ∂/∂r i . The initial conditions are small density and smooth velocity perturbations imposed on a homogeneous density distribution.
As long as the amplitude of the density perturbations remains small their evolution can be analyzed in the linear approximation obtained by the linearization of the above equations. The exact solution of the linearized system has one growing mode which is the major object of our analysis and two decaying modes that can be ignored. The velocity in the growing mode is a potential vector field which is proportional to the gradient of the linear gravitational potential: v lin ∝ −∇φ lin . In the linear regime the spatial structure of the perturbations (in the comoving coordinates) remains unchanged and its amplitude is proportional to the growing solution δ ∝ D(t) also determined by the cosmological model [6]. For example, D(t) ∝ a(t) ∝ t 2/3 in the Einstein -de Sitter model, in general D(t) is a monotonically growing function of time.
Equations 4-6 describing the evolution of inhomogeneities in the universe become more convenient for further analysis after the following transformation of variables are made [33] In addition, we shall use D instead of time, t. The resulting set of equations becomes Here Ω 0 =ρ m /ρ c at present time (it is the ratio of the total mean mass density to the critical value we also assume the summation over repeated indices. The second equation in set (8) becomes even simpler after introducing the total derivatives One can easily check that in the linear regime (i.e., while < δ 2 > 1/2 ≪ 1) the only growing solution to the set of equations (8) is given by where q is the Lagrangian coordinates of a fluid element. The solution (10) gives the density contrast, velocity and potential in terms of the Lagrangian coordinates. In order to find these fields in the Eulerian space one needs to solve the equation of trajectories for q = q(x, D) and substitute it in solution (10). The field Φ 0 (q) is determined by the perturbation of density in the linear regime which is assumed to be a realization of a Gaussian random field. The power spectrum of the field that solely defines all its statistical properties is derived from observations of the fluctuations of temperature in the Cosmic Microwave Background radiation (see e.g., [28]); it is shown in Fig.5.
Although velocity v i remains constant for each particle in the linear regime the physical velocity u p varies with time (Eq.7). In the growing solution (Eq.10) gravitational force proportional to −∂ϕ/∂x i , exactly balances the drag force due to the uniform expansion of the universe Hu p ∝ v and thus the right hand side of Eq.9 vanishes to linear order.

The Zel'dovich Approximation
Zel'dovich (1970) derived the Lagrangian form of the linear theory and extrapolated it to the beginning of the nonlinear regime [38] (see also [34]). This extrapolation assumes that the initial perturbations are smooth in the sense that the initial power spectrum must decrease faster than k −3 at k → 0. In realistic cosmological models e.g., in the ΛCDM model the initial power spectrum P (k) ∝ ln(k)k −3 up to very large k corresponding to cosmologically irrelevant scales and thus the Zel'dovich approximation must be properly modified in order to be useful on the scale of the clusters and superclusers of galaxies. This will be discussed later.
The Zel'dovich approximation can be conveniently formulated as a mapping from Lagrangian space L{q} to Eulerian space E{x} defined by Eq.11. A two-dimensional illustration of the mapping is shown at four stages of the evolution in Fig.3. The stages are marked by σ =< δ 2 lin > 1/2 where δ ≡ [ρ(x, t) −ρ(t)]/ρ(t) = [η(x, D(t)) −η]/η is the density contrast. By convention, the parameter σ is computed in the linear theory and therefore σ ∝ D(t).
Utilizing conservation of mass ηd 3 x =ηd 3 q Zel'dovich has derived the density as a function of D(t) and the Lagrangian coordinate q where λ 1 (q), λ 2 (q) and λ 3 (q) are the eigenvalues of the deformation tensor d ij = ∂ 2 Φ 0 /∂q i ∂q j . In cosmology the initial condition is usually characterized by the spectrum P δ (k) (Fig.5); it is obviously related to the power spectrum of the potential as P δ = k 4 P Φ 0 . Similarly to the linear stage one can find the density distribution in the Eulerian space by solving Eq.11 for q and substituting it in Eq.12.
For realistic initial conditions it requires numerical calculations. It follows from Eq.12 that the first objects are formed around the peaks of the largest eigenvalue (we assume that λ i are ordered at every q: λ 1 ≥ λ 2 and λ 2 ≥ λ 3 ) and have very oblate shapes because generic peaks of λ 1 always have two other eigen values different from λ 1 . These objects are known in cosmology as "Zel'dovich's pancakes". This prediction of the Zel'dovich approximation differs from the extrapolation of the Eulerian linear theory that predicts the formation of the first Figure 3: The distributions of particles at four stages of evolution described by the Zel'dovich approximation in two dimensions. The panels are marked by σ ≡< δ 2 lin > 1/2 , the variance of the density contrast evaluated to the linear order, so σ ∝ D(t).
objects from the peaks of δ = λ 1 + λ 2 + λ 3 . In practice the difference in most cases is not large due to obvious correlation between δ and λ 1 that causes the peaks of δ be spatially close to the peaks of λ 1 . Three-dimensional gravitational N-body simulations are in a good agreement with the prediction of the Zel'dovich approximation [35]. Pancakes in the collisionless dark matter originate as the three-stream flow regions bounded by caustics, the surfaces of formally infinite density. The shape and other characteristics of the pancakes as well as of other generic types of structures (Fig.3) are determined by catastrophe theory [2].
The Zel'dovich approximation proved to be very good until orbit crossing when caustics form and the multi-stream flows occur i.e., up to the stage corresponding to σ ≈ 1 (see e.g. [34] and references therein). At σ > 1 Eq.11 predicts broadening of the multi-stream flow regions considerably faster than it happens in N-body simulations [9,18,29] as can be seen in Fig. 4. Compare two panels on the right in Fig.3 with two panels of Fig. 4. Despite the difference in the realizations of the initial Figure 4: The distribution of particles at two stages σ = 1 on the left and σ = 2 on the right obtained in a two-dimensional N-body simulation. Adopted from [32].
conditions one can see that the correct gravitational gravitational force in the N-body simulations results in much thinner pancakes than predicted by the Zel'dovich approximation. This discrepancy motivated the development of the adhesion model that will be discussed in the following section.
Here we briefly describe a simple modification that allows to utilize the Zel'dovich approximation in the case when the small scale power cannot be ignored [7,25].
N-body simulations both in 2D and 3D have shown that the largest structures depend very little on small-scale perturbations [3,23,24]. The most important factor that affect the shapes and general appearance of the large-scale structure are the initial waves just reached nonlinear regime, i.e., on scales where rms fluctuations are just becoming nonlinear [23]. This observation allows to use an auxiliary model that yields the large-scale structures similar to the original model but has no smallscale structures [7,25]. The scale that separates large-scale structures from small-scale structure is the scale of nonlinearity.
For instance, the initial power spectrum in the currently popular ΛCDM cosmological model is shown in Fig.5 by the solid line while the spectrum filtered at the present scale of nonlinearity is shown by the dotted line. If the realization of amplitudes and phases are same in both models then the auxiliary model with filtered power spectrum would have a very similar large-scale structure to one in the model with unfiltered power spectrum. The dotted line in Fig.5 shows the spectrum filtered with a Gaussian window; a simpler version of filtering when the power spectrum set to zero on scales with k > k nl gives a similar result. The auxiliary model with filtered power spectrum provides an option of using the Zel'dovich approximation and therefore better understanding of complex nonlinear processes that control the structure formation. The price of this simplification consists in missing all structures on smaller scales that are abundant in the model with unfiltered initial spectrum. The adhesion model is an attempt to remedy this drawback of the Zel'dovich approximation.
Finally, it is worth noting that both the Zel'dovich approximation (Fig.3, top right panel) and N-body simulations (Fig.4) suggest that the universe can be viewed as an irregular tessellation where 'tiles' are represented by voids, the regions of low mass and/or galaxy densities, where the flow consists of only one stream.

The Adhesion Model
The general idea of the adhesion model is very simple. We wish to use the Zel'dovich approximation everywhere except the regions of multi-stream flows. By adding a diffusion term to the Euler equation one can avoid the orbit crossing and therefore suppress the formation of the multi-stream flow regions. The advantage of this modification is the suppression of too fast growth of the thickness of pancakes observed in the Zel'dovich approximation but it is achieved at the coast of noticeable change of the dynamics inside multistream flow regions. A small coefficient of viscosity guarantees agreement with the Zel'dovich approximation before the orbit crossing.
Since the motion is potential one can introduce the velocity potential v = −∇Φ. Then assuming that the gravitational potential is approximately equal to the velocity potential ϕ ≈ Φ and also adding a viscosity term ν∇ 2 v to the right hand side of the second equation in set 8 one obtains the equation of nonlinear diffusion [12,13] ∂v ∂D Generally speaking the viscosity term needs not to be in the form of Eq.13 but choosing this particular form one obtains Burgers' equation that has an exact analytic solution [5]. For the potential motion v = −∇Φ Eq.13 can be solved by performing the Hopf-Cole substitution Φ(x, D) = −2ν log U(x, D).
As a result Eq.13 adopts the form of the linear diffusion equation Solving Eq.14 for the velocity one obtains where the "action" Figure 6: Simulation of galaxy distribution in redshift space using the adhesion approximation with small but finite viscosity. Adopted from [37].
In cosmology the adhesion model has been used in two forms: one assumes a small but finite value of the viscosity parameter ν and the other assumes it is infinitesimal. For small but finite ν the trajectory of a particle can be determined by solving the integral equation and the resulting density can be determined from the continuity equation [36,37,29,26]. An example of the structure obtained by this method is shown in Fig.6 where the authors [37] plotted the simulated galaxy distribution in a slice through the simulation box similar to the plots obtained from galaxy redshift catalogs as in Fig.1. An interesting attempt to explain the magnitude of the coefficient of viscosity from dynamics in multistream flow regions was made in [4].
We will discuss a different approach that assumes the infinitesimal value of the viscosity parameter ν → 0 in more detail because it has direct connection to tessellation and tiling. In this case the Figure 7: The illustration of the geometrical technique that finds the Eulerian position of a particle using the initial potential Φ 0 (q). For a given value of D the paraboloid (Eq.17) descends on the surface of the initial potential Φ 0 until it is tangent to the surface at point q by adjusting constant P 0 . The projection of the apex of the paraboloid denotes the Eulerian position of the particle. Adopted from [20]. integrals in Eq.15 can be evaluated using the method of steepest descents [12,13,20,21]. The result has a remarkable geometrical interpretation illustrated by Fig.7 in two-dimensional case. The initial Lagrangian coordinates q p and the current Eulerian coordinates x p of a particle at a chosen stage denoted by D, can be related simply by descending a paraboloid with apex at x = x p and adjusting correspondently the constant P 0 until it is in contact with the surface of the initial potential Φ 0 (q) at some point. The coordinates of the contact point are the Lagrangian coordinates of the particle q p . This procedure must satisfy an additional constraint: the paraboloid must not cross Φ 0 at any point. Thus, the rule is: the tangential contacts of the paraboloid with the potential are welcome but crossings are strictly forbidden. At the linear stage D is small and the curvature of the paraboloid is greater than the typical curvature of Φ 0 , thus this condition can be easily fulfilled. The Zel'dovich approximation is universally valid at this stage and the velocity filed is given by the following equation Figure 8: An illustration of the geometrical prescription of descending a parabola onto the initial velocity potential in order to find the Eulerian positions of particles and knots in one-dimenaional case. The particle having Lagrangian coordinate q 0 in the uppermost panel has the Eulerian coordinate of the paraboloid apex.
In the middle panel corresponding to a later stage the whole region between two points where the parabola is tangent to the initial potential is stuck into the knot m 1 moving with the velocity shown by the arrow, its mass is m 1 =η(q 2 − q 1 ). The knot m 2 is determined by the second paraboloid in the middle panel. The lower panel shows the knot M formed as a result of merging m 1 , m 2 and the mass between them. Adopted from [32].
where q(x, D) is the coordinate of the absolute maximum of the action S(x, D; q) at given x and D.
As time passes the curvature of the paraboloid decreases, so it becomes wider. As a result the points q that are not accessible without violation of the crossing rule as can be seen in onedimensional illustration in Fig.8. Both parabolas shown in the middle panel cannot descend further without crossing the initial potential Φ 0 . What has happened to the corresponding particles? The particles between two contact points i.e., between q 1 and q 2 for one parabola and between q 3 and q 4 for the other have been stuck into point masses m 1 and m 2 respectively due to orbit crossing. Their coordinates are the coordinates of the apices of the corresponding parabolas and the masses are m 1 =η(q 2 − q 1 ) and m 2 =η(q 4 − q 3 ). At later times the discrete masses formed earlier can grow by merging as well as by accreting the continuous medium between them. This is illustrated by the bottom panel. It is not hard to see that the transition from the state shown in the middle panel to one shown in the bottom panel requires the state when one parabola has three contacts with the The panel on the right shows the tessellation computed in the adhesion model from the same initial conditions. The tessellation is superimposed on the particle distribution. Adopted from [21]. potential simultaneously. At the next instant of time it becomes a little wider and looses the contact with the central peak.
In two dimensions the paraboloid can be tangent to the initial potential Φ 0 simultaneously in one, two or three points. The first type corresponds to the particles that form continuous medium. The points of the second type have been stuck into lines shown by thick points in Fig.7. These lines, the edges of tiles, are formed by the set of the apices of the paraboloids that are in contact with the initial potential at two points simultaneously. Finally, when the paraboloid is in contact with the surface of Φ 0 at three point simultaneously its apex marks the knot (vertex) where the lines (edges) meet. The edges and vertices form the boundaries of irregular tiles of tessellation in two dimensions. In the course of the evolution some tiles expand while the other collapse and cease to exist. This goes hand in hand with merging of the knots. At the instant of a merger the paraboloid can be in contact with the initial potential at more than three points simultaneously [1]. Figure 9 demonstrates the agreement between the tessellation computed in the adhesion model and the N-body simulations in two-dimensional space. Figure 10 illustrates the effect of small-scale power. The initial potentials in both panels were generated with the same amplitudes and phases but the amplitudes for the short waves k > k cr were set to zeros in the case shown in the left hand side panel. Therefore, the initial potential in the left hand side panel is a smooth version of the potential in the right hand side panel. The both models evolved to the same stage. One can see that while the overall appearance of the tessellation is similar in both panels the there more small-scale details in the panel on the right.
With a little imagination and help from the geometric model outlined above one can easily visualize the main features of the tessellation in three-dimensional space. The initial potential Φ 0 forms a three-dimensional hypersurface in four-dimensional space. A three-dimensional paraboloid (Eq.17) descends on the thre-dimensional hypersurface of the initial potential similar to one-and Figure 10: Two panels show the tessellation of the plane by the adhesion model superimposed on the particle distribution obtained in 2D N-body simulation. The initial potential Φ 0 in both panels have the same large-scale part but on the left hand side it does not have the small-scale part that is present in the potential on the right hand side. Adopted from [21].
two-dimensional cases. The rules are the same: it can touch the initial potential but must not cross it. Four generic situations are possible at a generic instant of time: it can have only one, two, three or four contact points simultaneously. If there is only one contact point then the corresponding particle is in free motion; the coordinates of the of contact point are the Lagrangian coordinates of the particle and the apex denores its Eulerian coordinates. In the remaining three cases the apex define a particle stuck on a surface, on an edge or at a vertex of the tessellation, if the number of contact points is two, three or four respectively. Again similarly to one-and two-dimensional cases at the critical instants of time when a tree-dimensional tile collapses and ceases to exist or two or more vertices merge the number of contact points can be greater than four [1].
The set of surfaces, filaments and knots define the tessellation of three-dimensioal space: voluminous low density tiles (voids) are separated by the faces (pancakes), the edges (filaments) are at the intersections of the faces, and the vertices (clusters of galaxies, halos in N-body simulations) are at the intersections of the edges. This tessellation can be viewed as the skeleton of the real structure.
The both variants of the adhesion model (with small but finite and infinitesimal viscosity) have been tested against the gravitational N-body simulations in two and three dimensions. Both the N-body simulation and the adhesion model used the identical initial conditions and were compared at several stages of the evolution. The geometrical version of the adhesion model was tested against the two-dimensional N-body simulations with the initial power law spectra P δ (k) ∝ k n with spectral indices n = 2, 0 and −2 and various cutoffs [21]. The version of the adhesion model utilizing a finite viscosity parameter ν has been quantitatively compared to fully nonlinear, numerical threedimensional gravitational N-body simulations in [26]. The results generally show a good agreement however they improve with the growth of the relative power on large scales. This is the case in the ΛCDM model as one can see from the initial power spectrum shown in Fig.5. Figure 11: Normalized initial gravitational potential at the centers of voids is plotted against the void diameter. Adopted from [32].

Summary
The adhesion model based on three-dimensional Burgers' equation (Eq.13) of the nonlinear diffusion has been found to provide a good explanation of the major large-scale features of the galaxy distribution in the universe. Its particular version assuming the infinitesimal coefficient of viscosity ν → 0 provides a natural geometrical construction that can be characterized as "quasi-Voronoi tessellation" [27]. Comparison with gravitational two-and three-dimensional N-body simulations has shown a good agreement including advanced nonlinear stage. The adhesion model provides a natural qualitative explanation of the origin of the large-scale coherent structures such as superpancakes and superfilaments, as a result of coherent motion of clumps due to large-scale inhomogeneities in the initial gravitational potential [19]. The initial gravitational potential also affects the sizes of voids as shown in Fig.11: the higher the peak of the initial potential the greater size of the voids formed around it. The formation and evolution of large-scale structure is described by the adhesion model as a two stage process [32]. During the first short stage matter falls predominantly onto pancakes and then moves along them towards filaments and then along filaments to the knots. At the end of the first stage the formation of the tessellation is completed and the bulk of the mass in the universe is located in knots (vertices) and filaments (edges) however a relatively small fraction of mass remains in pancakes (faces) and voids (three-dimensional tiles). The second stage is characterized by the deformation of the tiles and edges due to the dynamics caused by gravity . At this stage knots merge into more massive knots and small voids shrink and cease to exist giving space to growing large voids. Eventually almost all the mass concentrates in knots. Depending on the initial spectrum the knots themselves may move coherently in such a manner that they concentrate to superpancakes and superfilaments. This situation is likely to occur in the ΛCDMmodel at the present time due to the shape of the initial power spectrum (Fig.5). The superpancakes and superfilaments can be identified by applying the adhesion model to smoothed initial potential [21].