Chapter 8

Park’s Simulation of the Universe

I couldn’t wait to get home to Princeton after hearing in Rio about the Great Wall of galaxies. I told my graduate student at that time, Changbom Park (currently Professor at the Korean Institute for Advanced Study), that he should stop what he was doing and devote full time to making the world’s largest computer simulation. Time was short because everyone had seen that picture. Our department chairman, Jerry Ostriker, had just bought our department a supercomputer—so doing a new, largest simulation was possible. In only three weeks, Changbom wrote an N-body code to do the calculation. This was remarkably fast. The new simulation would involve 4 million particles. In 1975 Sverre Aarseth, Ed Turner, and I had created a 4,000-body simulation of galaxy clustering in the universe, breaking the existing record. Changbom would publish his thesis in 1990; thus, in 15 years, the number of particles had increased by a factor of 1,000. We joked that in another 15 years (by 2005), we should be doing 4 billion particles! (It seemed fantastic at the time, and yet Volker Springel and his colleagues would actually do a 10-billion-particle simulation by 2005.) Our extrapolation was consistent with Moore’s law, the empirical observation that computer power seems to double every 1.5 years. After 4.5 years (three doublings) we should see an increase of a factor of 8. Five years gives an increase of about a factor of 10, and in 15 years, if it continues in the same way, that should be three factors of 10, or a factor of 1,000. Perhaps no field has benefited more from the computer revolution than astronomy. Computers allow us to analyze ever-larger observational samples, and they also allow us to do ever-larger simulations to model them.

Changbom would be using Jim Peebles’s theory of cold dark matter and the Big Bang inflationary model as the basis for his initial conditions. He would lay down small-density perturbations as waves. Recall our analogy of ocean waves, where a wave crest represents a region of slightly higher than average density and a wave trough represents a region of slightly lower than average density. A single wave would be sinusoidal, with crests and troughs lined up as parallel furrows. The wave would have an amplitude (in density), a wavelength (distance between successive wave crests), and a direction. Changbom was going to lay down the waves at random. Each wave would be given a random direction and a random phase; that is, the wave crests and troughs would point in a random direction, and they would be shifted a random distance along the direction of the wave, so that a particular point would be equally likely to be part of a wave crest or a wave trough. The cold dark matter theory coupled with inflation would tell Changbom the mean-square amplitude of the waves at a given wavelength; then he would pick the amplitude of the particular wave randomly from a bell-shaped (Gaussian) distribution with that mean-square amplitude. This gives what is called a Gaussian random-phase distribution—the initial conditions predicted by inflation. This distribution gives a spongelike topology of high- and low-density regions in the initial conditions, as we have discussed. Inflationary theory tells us how much power (square of the amplitude of the waves) there is on average in fluctuations as a function of their wavelength. We call this the power spectrum. Just as a spectrum gives us information on how much light is present as a function of wavelength, the power spectrum tells us how much power there is on average in density fluctuations of different wavelengths. Inflation predicts the amplitude of fluctuations on different scales.

Using Peebles’ cold dark matter theory, one could follow the growth of the fluctuations through the action of gravity up through the period of recombination and after. Fluctuations then grew by gravity. Using this theory as a basis, Changbom was able to start his simulations when the universe was about a factor of 24 smaller than it is today. All the density fluctuations at this point were small with respect to unity (i.e., a few percent), and they were already growing due to gravity. Two million matter particles were placed on a 128 × 128 × 128 cubical grid and displaced slightly to produce the density fluctuations. They were given the individual initial velocities that they would have attained due to the growth of those fluctuations. Changbom assumed a density of matter that was 40% of that required to ever halt the expansion of the universe, close to the value of 30.8% measured by the Planck Satellite Collaboration in 2014. Most of this matter is in the form of cold dark matter, with the rest contributed by ordinary atoms. The computer then calculates the movement of these matter particles due to their mutual gravitational interactions. Essentially, Isaac Newton takes over. Slightly overdense regions become denser still as their excess gravity draws other particles in, while underdense regions grow into voids. The computer follows the complicated nonlinear process as clusters of galaxies and dense filaments form. The computer program used a fast Fourier technique to cut the computation time in calculating the gravitational forces on the cold dark matter particles. This technique had been pioneered by R. W. Hockney and J. W. Eastwood, Kevin Prendergast and R. H. Tomer, George Efstathiou, and Adrian Melott, among others. Standard computations of gravitational forces (such as those Peebles and Aarseth used) required of order N² operations per time step, whereas the fast Fourier technique required only of order N ln N operations. For N = 2 million, for example, N² = 4 × 10¹², whereas N ln N = 2.9 ×10⁷, for a considerable savings.

In the initial conditions, Changbom also introduced a second, biased set of 2 million particles placed to represent locations where individual galaxies were likely to form. It is recognized that galaxies are more likely to form in regions of higher cold dark matter density. The ordinary matter tends naturally to fall into such regions and form galaxies there. Likewise, galaxies are even less likely to form in voids where the matter density is low. (Observationally, today we find fluctuations of order unity in the number of galaxy counts in spheres of radii 39 million light-years, whereas fluctuations in the dark matter density on the same scale are only 80% as large.) So, in regions of higher-than-average matter density in the initial conditions, more “galaxy” particles were placed, and in regions of lower-than-average matter density fewer galaxy particles were placed. James Bardeen and his colleagues had developed formulas for doing this. These galaxies were then allowed to move due to the gravitational influence of the matter particles. As the matter particles formed clusters, the galaxy particles fell into them. As the universe simulation expanded by a factor of 24 to reach the present epoch, the full range of galaxy clustering occurred.

The simulation was divided into 256 × 256 × 256 cells on which the gravitational field was calculated. The boundary conditions on the box were periodic, front matching the back, left side matching the right side, and top matching the bottom, so Changbom was effectively modeling an infinite universe made of identical stacked boxes. No place in the cubical box was special. This is fine for modeling structures smaller than the cubical box size. The whole simulation modeled a cube 1.49 billion light-years on a side at the present epoch, large enough to encompass structures as large as the Great Wall. But would such a structure appear?

Finding a Great Wall in the Simulation

We were quite excited to look at the first pictures. Changbom just sliced the cube in thick slices like a loaf of bread—and there it was, a structure just as magnificent as the Great Wall! We found it right away looking at the 3D sample. The next step was to produce a simulated fan-shaped survey that would mimic the original 6°-thick fan published by de Lapparent, Geller, and Huchra. We simply placed an “Earth” in our simulation box and considered how the universe would “look” from that vantage point. We surveyed galaxies in a fan-shaped slice of the universe of the same size and shape as they had observed. Changbom gave the simulated galaxies brightnesses taken from the distribution found in the real universe and included only the simulated galaxies that would have been bright enough for de Lapparent, Geller, and Huchra to have seen from this Earth. The simulations are compared with the real de Lapparent, Geller, and Huchra slice in Figure 8.1.

The agreement was extraordinary! Both pictures showed large, empty voids surrounded by thin filaments and a prominent filament, or wall, extending from left to right all the way across the survey. This agreement was not the result of any extensive trial-and-error search through the data cube. We simply found a long filament immediately in the 3D data cube, placed an observation point at approximately the right distance from it, and took a snapshot picture.

The simulation also includes a great cluster, just like the Coma cluster, which appears as a sharp dagger pointed at Earth, sitting at the right-hand end of the simulated Great Wall. In the de Lapparent, Geller, and Huchra slice, the Coma cluster of galaxies appears as a comparable downward-pointing dagger in the Great Wall right at the center of the fan. We also made thicker slices, and these showed Changbom’s simulated Great Wall even more dramatically—just as occurred in the case of the Great Wall of Geller and Huchra.

Figure 8.1. De Lapparent, Geller, and Huchra’s slice (a) is compared with Changbom Park’s simulated slice (b). Park’s computer simulation using Cold Dark Matter and biased galaxy formation produced a Great Wall and large voids remarkably like those seen in the de Lapparent, Geller, and Huchra observational slice. (Credit: Changbom Park, Monthly Notices of the Royal Astronomical Society, 242: P59, 1990)

Changbom Park’s thesis, published in the Monthly Notices of the Royal Astronomical Society in 1990, proved that a cold dark matter inflationary model could produce structures as dramatic as the Great Wall of Geller and Huchra. It reproduced with remarkable fidelity the pattern of filaments and voids that they had observed. Gravity alone could produce these fantastic structures from the small initial fluctuations predicted by inflation. Exotic explanations, like nuclear explosions or cosmic string wakes, were not needed. One thing we knew about our simulations was that their galaxy clustering had a spongelike topology—clusters connected by high-density filaments and voids connected by tunnels. At this point we had only Geller and Huchra’s 2D slice, not a fully 3D sample to work with, so we could not check their sample’s 3D topology yet, but they were continuing to add to their slice, making it thicker and thicker, so eventually we would have a nice 3D volume to check for topology. For the time being, we knew that our simulations, which had a spongelike topology of large-scale structure, did produce the look of the Geller and Huchra slices perfectly. One did not require a froth of bubbles with a Swiss cheese topology to produce the Geller and Huchra slices.

The Great Attractor

Changbom Park and I continued to make larger and larger simulations and check them against new observations. In 1998, a group of seven astronomers (Donald Lynden-Bell, Sandra Faber, David Burstein, Roger Davies, Alan Dressler, R. Turlevich, and Gary Wegner) had discovered something they named the Great Attractor. This group of astronomers would become known as the “Seven Samurai.” To measure distances of elliptical galaxies, they made use of an observed relation: the higher the velocity dispersion (or range of velocities) of the stars orbiting within an elliptical galaxy, the bigger in physical size the elliptical galaxy would be. The Seven Samurai could measure the velocity dispersion of the stars by measuring the broadening of the spectral lines in the galaxy. Different stars traveling at different speeds relative to us would create different redshifts for the spectral lines of those stars, smearing out, or broadening, the spectral lines seen in the spectrum of the galaxy as a whole. Once they knew the velocity dispersion of the stars in the galaxy, they could deduce its physical size. If they knew its physical size and then measured its angular size in the sky, they could deduce its actual distance. An object of a given physical size has a smaller angular size if viewed from a larger distance. If you know how big a “stop sign” really is, its angular size as you look out your windshield will tell you how far away it is. In this way the team could estimate the distances to many elliptical galaxies. If the distance to the galaxy differed from that deduced via Hubble’s law from its redshift, it meant that the galaxy had a peculiar, or individual, motion over and above the general Hubble expansion. Our own galaxy, for example, has a peculiar velocity of about 550 kilometers/second relative to the cosmic microwave background. These peculiar velocities are, in the standard Big Bang inflationary cosmological model, produced by the extra gravitational attraction of nearby regions of above-average density (versus the lessened gravitational attraction of other nearby regions of below-average density). The Seven Samurai found that galaxies showed large peculiar velocities, all pointing to a Great Attractor located at a distance from us of about 200 million light-years in the constellation of Centaurus. They calculated that this must be due to an excess mass in this region of 5 × 10¹⁶ solar masses. They noted a concentration of galaxies seen in this region that was 20 times larger than the nearby Virgo Cluster (about 50 million light-years from us). The appearance of a Great Attractor was surprising. If the voids in the Geller and Huchra slices were produced by explosions, one might expect instead to find “Great Repulsors”—regions in the centers of voids from which galaxies were fleeing. Instead, when velocity flows in the universe were investigated, the first thing to appear was a Great Attractor.

Changbom and I decided to look for peculiar velocity flows in our simulations. We measured the gravitational potential energy at all points in our large cubic simulation. The low point in the gravitational potential—the most gravitationally bound spot in the entire simulation—should be a place toward which galaxies migrate. We knew the peculiar velocities of all the galaxies in our simulation, and so we put a slice through the most gravitationally bound spot and put peculiar velocity vector arrows on each galaxy showing where it was headed. We found a dramatic thicket of long arrows pointing right at the most gravitationally bound spot in the simulation—a Great Attractor! But why were there no Great Repulsors? We looked for the spot in the simulation that was highest in gravitational potential energy, the least gravitationally bound spot in the entire simulation. This should be a place from which galaxies were fleeing. But when we looked at a slice through the simulation that included this spot, we saw no thicket of long arrows headed away from it: no Great Repulsor. Why? It was simply because the least bound spot was located in the middle of a giant void, with no galaxies nearby to show the motion! It was, therefore, the lack of tracer galaxies nearby that made it unlikely for astronomers to discover a Great Repulsor. But a Great Attractor was demonstrably easy to find. Thus, the computer simulations had explained another, seemingly contradictory, discovery.

The velocity flows around the Great Attractor have now been mapped in an elegant way by R. Brent Tully, Helene Courtois, Yehuda Hoffman, and Daniel Pomarède (2014). At each position in space the peculiar, or individual, velocities of individual galaxies can be plotted as arrows using a surveying technique developed by Ed Bertschinger and Avishai Dekel (1989).¹ Tully and his colleagues then connect up the arrows to draw flow lines showing the direction of motion. These velocity flows are due to the gravitational attraction of neighboring groups and clusters of galaxies.

We may understand this by thinking about rivers and their watersheds on Earth. If you let a drop of water fall on the ground in a natural terrain, it will roll downhill, pulled by gravity. If you keep following its path, you will flow into a river and eventually arrive at the sea. The set of points where water can fall and eventually find its way into the sea at the mouth of the Mississippi river, for example, is called the watershed of the Mississippi. It is the total area drained by the Mississippi River. Tully and his colleagues have mapped the peculiar velocity flow lines of galaxies in our neighborhood. They show a giant watershed of flow lines leading to the Great Attractor, as if they were rivers (see Color Plate 5). We are part of this watershed. They have named this region the Laniakea Supercluster. (Laniakea means “immeasurable heaven” in Hawaiian.) We are at the outer edge of this region, which is 510 million light-years across.

We have a new addition to our cosmic address: Earth, Solar System, Milky Way, Local Group, Virgo Supercluster, Laniakea Supercluster. The Virgo Supercluster is just our nearby branch of the Laniakea Supercluster. Does this mean that our galaxy will eventually fall into the Great Attractor? No, the individual velocity of our galaxy moving us in its direction is small relative to the average velocity of separation we are experiencing relative to it, due to the overall expansion of universe. Moreover, the expansion of the universe is accelerating, so the Great Attractor is actually fleeing from us ever more rapidly. But the Laniakea Supercluster marks the extent of the gravitational influence of the Great Attractor.

The Coma cluster is part of a different, neighboring watershed called the Coma Supercluster. Tully and colleagues identified two other independent watersheds: the Shapley Supercluster and the Perseus-Pisces Supercluster. We can divide the volume of space into these gravitational watersheds. This presents a new challenge for N-body computer simulations. We can use the same velocity flow techniques to map out watersheds in the N-body simulations and see if they look just like the ones we find in the universe. That will be an interesting area for future research. Given our ability to find Great Attractors in the simulations so easily, I expect N-body simulations to pass this test as well.

Rodlike Simulation

In 1990, Broadhurst and colleagues did two deep, pencil-beam surveys pointed in opposite directions in the sky to produce a rod-shaped survey with a length of about 10 billion light-years centered on Earth. They found a succession of walls in the survey that appeared to be approximately equally spaced, about 600 million light-years apart. That did not seem to fit in well with the idea of random initial conditions produced by inflation. Changbom and I immediately did a long cosmological rodlike simulation. It simulated a volume of 12 billion light-years × 490 million light-years × 490 million light-years (see Figure 8.2).

In this rod-shaped volume we also discovered a series of walls, and in simulating 12 narrow pencil-beam surveys, we found one survey that was even more nearly periodic in appearance than the survey by Broadhurst and colleagues. Thus, the periodic appearance of walls in the observational data was not statistically abnormal at the 95% confidence level in a cold dark matter inflationary cosmology. Random initial conditions with cold dark matter and inflation naturally produced a in time. The double arrows indicate a length of 973 million light-years. This section of the sequence of walls in pencil-beam surveys that were just as dramatic as those found observationally.

Figure 8.2. One-eighth of the Park-Gott simulation of a rod-shaped volume in the universe. Galaxies at the top are at a redshift of 0.235; galaxies at the bottom are at a redshift of 0.384. Scanning from top to bottom, one is looking further out in the universe and further back in time. The double arrows indicate a length of 973 million light-years. This section of the rod-shaped volume is 1.46 billion light-years tall, 487 million light-years wide, and 97 million light-years thick. One sees giant voids and great walls. (Credit: Changbom Park and J. Richard Gott, III, Monthly Notices of the Royal Astronomical Society, 249: 288, 1991)

Overall, the simulations have been remarkably successful. Thin slices of the universe show voids most prominently; the same holds true for simulations. Thicker slices show Great Walls—likewise for simulations. Velocity flows in the universe show Great Attractors, just as in the simulations. Deep pencil-beam surveys show a succession of Great Walls, as do the simulations. In cosmology, the larger the observational samples have become, and the larger the simulation volumes have become to match them, the more spectacular has been the agreement between the two. This is a sign we are on the right track.