text:
# Simulating the formation of cosmic structure

## 1 Introduction

The origin of structure in the Universe is a central problem in Physics. Its solution will not only inform our understanding of the processes by which matter became organized into galaxies and clusters, but it will also help uncover the identity of the dark matter, offer insights into events that happened in the early stages of the Big Bang and provide a useful check on the values of the fundamental cosmological parameters estimated by other means.

Because of its non-linear character, lack of symmetry and general complexity, the formation of cosmic structure is best approached theoretically using numerical simulations. The problem is well posed because the initial conditions – small perturbations in the density and velocity field of matter – are, in principle, known from Big Bang theory and observations of the early Universe, while the basic physical principles involved are understood. The behaviour of the dark matter is governed primarily by gravity, while the formation of the visible parts of galaxies involves gas dynamics and radiative processes of various kinds. Using cosmological simulations it is possible to follow the development of structure from primordial pertubations to the point where the model can be compared with observations.

Over the past few years, there has been huge progress in quantifying observationally the properties of galaxies not only in the nearby universe, but also in the very distant universe. Since the clustering pattern of galaxies is rich with information about physics and cosmology, much effort is invested in mapping the distribution of galaxies at different epochs. Two large ongoing surveys, the US-based Sloan Digital Sky Survey (York et al. 2000), and the Anglo-Australian “2-degree field galaxy redshift survey” (2dFGRS, Colless et al. 2001), are revolutionizing our view of the nearby universe with order of magnitude increases in the amount of available data. Similarly, new data collected in the past five years or so have, for the first time, opened up the high redshift universe111In cosmology, distances to galaxies are estimated from the redshift of their spectral lines; higher redshifts correspond to more distant galaxies and thus to earlier epochs. to detailed statistical study (Steidel et al. 1996).

The advent of large computers, particularly parallel supercomputers, together with the development of efficient algorithms, has enabled the accuracy and realism of simulations to keep pace with observational progress. With the wealth of data now available, simulations are essential to interpret astronomical data and to link them to physical and cosmological theory.

## 2 Building a model

To build a model of large-scale structure, four key ingredients need to be specified: (i) the content of Universe, (ii) the initial conditions, (iii) the growth mechanism, and (iv) the values of fundamental cosmological parameters. I now discuss each of these in turn.

### 2.1 The content of the Universe

Densities are usually expressed in terms of the cosmological density parameter, $`\mathrm{\Omega }=\rho /\rho _{crit}`$, where the critical density, $`\rho _{crit}`$, is the value that makes the geometry of the Universe flat. The main constituents of the Universe and their contribution to $`\mathrm{\Omega }`$ are listed in Table 1.

The main contribution to the extragalactic radiation field today is the cosmic microwave background (CMB), the redshifted radiation left over from the Big Bang. These photons have been propagating freely since the epoch of “recombination”, approximately 300,000 years after the Big Bang. The CMB provides a direct observational window to the conditions that prevailed in the early Universe. The Big Bang also produced neutrinos which today have an abundance comparable to that of photons. We do not yet know for certain what, if any, is the mass of the neutrino, but even for the largest masses that seem plausible at present, $`0.1`$eV, neutrinos make a negligible contribution to the total mass budget (although they could be as important as baryons). The abundance of baryons is now known with reasonable precision from comparing the abundance of deuterium predicted by Big Bang theory with observations of the absorption lines produced by intergalactic gas clouds at high redshift seen along the line-of-sight to quasars (Tytler et al. 2000). Baryons, the overwhelming majority of which are not in stars today, are also dynamically unimportant (except, perhaps, in the cores of galaxies).

Dark matter makes up most of the matter content of the Universe today. To the now firm dynamical evidence for its existence in galaxy halos, even more direct evidence has been added by the phenomenon of gravitational lensing which has now been detected around galaxy halos (e.g. Fischer et al. 2000, McKay et al. 2001, Wilson et al. 2001), in galaxy clusters (e.g. Clowe et al. 2000), and in the general mass field (e.g. Van Waerbeke et al. 2001 and references therein). The distribution of dark matter in rich clusters can be reconstructed in fair detail from the weak lensing of distant background galaxies in what amounts virtually to imaging the cluster dark matter. Various dynamical tests are converging on a value of $`\mathrm{\Omega }_{\mathrm{dm}}0.3`$, which is also consistent with independent determinations such as those based on the baryon fraction in clusters (White et al. 1993, Evrard 1997), and on the evolution in the abundance of galaxy clusters (Eke et al. 1998, Borgani et al. 2001). Since $`\mathrm{\Omega }_{\mathrm{dm}}`$ is much larger than $`\mathrm{\Omega }_b`$, it follows that the dark matter cannot be made of baryons. The most popular candidate for the dark matter is a hypothetical elementary particle like those predicted by supersymmetric theories of particle physics. These particles are referred to generically as cold dark matter or CDM. (Hot dark matter is also possible, for example, if the neutrino had a mass of $`5`$ eV. However, early cosmological simulations showed that the galaxy distribution in a universe dominated by hot dark matter would not resemble that observed in our Universe (White, Frenk and Davis 1983).)

A recent addition to the cosmic budget is the dark energy, direct evidence for which was first provided by studies of type Ia supernovae (Riess et al. 1998, Perlmutter et al. 1999)222The possibility that dark energy might be the dynamically dominant component had been anticipated by theorists from studies of the cosmic large-scale structure (see e.g. Efstathiou et al. 1990), and was considered in the first simulations of structure formation in cold dark matter universes (Davis et al. 1985).. These presumed ‘standard candles’ can now be observed at redshifts between 0.5 and 1 and beyond. The more distant ones are fainter than would be expected if the universal expansion were decelerating today, indicating that the expansion is, in fact, accelerating. Within the standard Friedmann cosmology, there is only one agent that can produce an accelerating expansion. This is nowadays known as dark energy, a generalization of the cosmological constant first introduced by Einstein, which could, in principle, vary with time. The supernova evidence is consistent with the value $`\mathrm{\Omega }_\mathrm{\Lambda }0.7`$. Further, independent evidence for dark energy is provided by a recent joint analysis of CMB data (see next section) and the 2dFGRS (Efstathiou et al. 2002).

Amazingly, when all the components are added together, the data are consistent with a flat universe:

$$\mathrm{\Omega }=\mathrm{\Omega }_b+\mathrm{\Omega }_{\mathrm{dm}}+\mathrm{\Omega }_\mathrm{\Lambda }1$$
(1)
### 2.2 The initial conditions

The idea that galaxies and other cosmic structures are the result of the slow amplification by the force of gravity of small primordial perturbations present in the mass density at early times goes back, at least, to the 1940s (Lifshitz 1946). However, it was only in the early 1980s that a physical mechanism capable of producing small perturbations was identified. This is the mechanism of inflation, an idea due to Guth (1981), which changed the face of modern cosmology. Inflation is produced by the dominant presence of a quantum scalar field which rolls slowly from a false to the true vacuum, maintaining its energy density approximately constant and causing the early Universe to expand exponentially for a brief period of time. Quantum fluctuations in the inflaton field are blown up to macroscopic scales and become established as genuine adiabatic ripples in the energy density. Simple models of inflation predict the general properties of the resulting fluctuation field: it has Gaussian distributed amplitudes and a near scale-invariant power spectrum (Starobinskii 1982).

After three decades of ever more sensitive searches, evidence for the presence of small fluctuations in the early universe was finally obtained in 1992. Since prior to recombination the matter and radiation fields were coupled, fluctuations in the mass density are reflected in the temperature of the radiation. Temperature fluctuations in the CMB were discovered by the COBE satellite (Smoot et al. 1992) and are now being measured with ever increasing accuracy, particularly by detectors deployed in long-flight balloons (de Bernardis et al. 2000, Hanany et al. 2000, Leitch et al. 2002). The spectrum of temperature fluctuations is just what inflation predicts: it is scale invariant on large scales and shows a series of “Doppler” or “acoustic” peaks which are the result of coherent acoustic oscillations experienced by the photon-baryon fluid before recombination. The characteristics of these peaks depend on the values of the cosmological parameters. For example, the location of the first peak is primarily determined by the large-scale geometry of the Universe and thus by the value of $`\mathrm{\Omega }`$. Current data imply a flat geometry, consistent with eqn. 1.

The spectrum of primordial fluctuations generated, for example, by inflation evolves with time in a manner that depends on the content of the Universe and the values of the cosmological parameters. The dark matter acts as a sort of filter, inhibiting the growth of certain wavelengths and promoting the growth of others. Following the classical work of Bardeen et al. (1986), transfer functions for different kinds of dark matter (and different types of primordial fluctuation fields, including non-Gaussian cases) have been computed. In Gaussian models, the product (in Fourier space) of the primordial spectrum and the transfer function, together with the growing mode of the associated velocity field, provides the initial conditions for the formation of cosmic structure.

### 2.3 Growth mechanism

Primordial fluctuations grow by gravitational instability: overdense fluctuations expand linearly, at a retarded rate relative to the Universe as a whole, until eventually they reach a maximum size and collapse non-linearly to form an equilibrium (or ‘virialized’) object whose radius is approximately half the physical size of the perturbation at maximum expansion. The theory of fluctuation growth is lucidly explained by Peebles (1980).

Although gravitational instability is now widely accepted as the primary growth mechanism responsible for the formation of structure, it is only very recently that firm empirical evidence for this process was found. Gravitational instability causes inflow of material around overdense regions. From the perspective of a distant observer, this flow gives rise to a characteristic infall pattern which is, in principle, measurable in a galaxy redshift survey by comparing the two-point galaxy correlation function along and perpendicular to the line-of-sight. In this space, the infall pattern resembles a butterfly (Kaiser 1987). This pattern has been clearly seen for the first time in the 2dFGRS (Peacock et al. 2001)333Strictly speaking the ‘butterfly’ pattern does not prove the existence of infall since the continuity equation would ensure a similar pattern even if velocities were induced by non-gravitational processes. However, it can be shown that such velocities, if present, would rapidly decay..

### 2.4 Cosmological parameters

After decades of debate, the values of the fundamental cosmological parameters are finally being measured with some degree of precision. The main reason for this is the accurate measurement of the acoustic peaks in the CMB temperature anisotropy spectrum whose location, height and shape depend on the values of the cosmological parameters. Some parameter degeneracies exist but some of these can be broken using other data, for example, the distant Type Ia supernovae or the 2dFGRS (eg. Efstathiou et al. 2002). The CMB data alone do not constrain the Hubble constant, but there is a growing consensus from the HST key project (Freedman et al. 2001), and other methods, that its value, in units of 100 km s-1 Mpc,-1 is $`h=0.70\pm 0.07`$. In addition to $`h`$ and the other parameters listed in Table 1, the other important number in studies of large-scale structure is the amplitude of primordial density fluctuations which is usually parametrized by the quantity $`\sigma _8`$ (the linearly extrapolated value of the top-hat filtered fluctuation amplitude on the fiducial scale of 8 $`h^1\mathrm{Mpc}`$). The best estimate of this quantity comes from the observed abundance of rich galaxy clusters which gives $`\sigma _8\mathrm{\Omega }^{0.6}=0.5`$, with an uncertainty of about 10% (Eke, Cole & Frenk 1996, Viana & Liddle 1996, Pierpaoli et al. 2001).

## 3 Cosmological simulations

Operationally, the problem of the cosmic large-scale structure can be divided into two parts: understanding the clustering evolution of the dark matter and understanding the gaseous and radiative processes that lead to the formation of galaxies. Specialized simulation techniques have been developed to tackle both aspects of the problem. The evolution of the dark matter is most often calculated using N-body techniques, implemented through a variety of efficient algorithms, such as P3M (Particle-particle/particle-mesh; Efstathiou et al. 1985), AP3M (the adaptive mesh version of P3M; Couchman et al. 1995) and hierarchical trees (Barnes & Hut 1986, Springel et al. 2001, Stadel 2000). Gaseous and radiative processes are followed by combining a hydrodynamics code with an N-body code. Numerical hydrodynamic techniques used in cosmology include Eulerian methods (Cen 1992), Lagrangian codes based on Smooth Particle Hydrodynamics (SPH) (Gingold & Monaghan 1977), and hybrid codes (e.g. Gnedin 1995, Pen 1998). These techniques have different strengths and weaknesses, but they all give similar results in the simplest cosmological problems where a detailed comparison has been performed (Frenk et al. 1999).

There has been a rapid growth in the size and power of cosmological simulations in the two and a half decades since this technique was introduced into the subject by Peebles (1970). One way to measure this growth is by the number of particles employed in the simulations. The size of the largest simulations has grown exponentially, in a manner reminiscent of the well-known “Moore’s law” that describes the increase in cpu speed with time, except that the advent of massively parallel supercomputers led to a sudden order-of-magnitude jump in size towards the end of the past decade. The largest simulations carried out to date are the 1-billion particle “Hubble volume,” N-body simulations performed by the Virgo consortium, an international collaboration of reseachers in the UK, Germany and Canada.

### 3.1 Large-scale structure

Figure 1 illustrates the spatial distribution of dark matter at the present day, in a series of simulations covering a large range of scales. Each panel is a thin slice of the cubical simulation volume and shows the slightly smoothed density field defined by the dark matter particles. In all cases, the simulations pertain to the “$`\mathrm{\Lambda }`$CDM” cosmology, a flat cold dark matter model in which $`\mathrm{\Omega }_{\mathrm{dm}}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`h=0.7`$. The top-left panel illustrates the Hubble volume simulation: on these large scales, the distribution is very smooth. To reveal more interesting structure, the top right panel displays the dark matter distribution in a slice from a volume approximately 2000 times smaller. At this resolution, the characteristic filamentary appearance of the dark matter distribution is clearly visible. In the bottom-right panel, we zoom again, this time by a factor of 5.7 in volume. We can now see individual galactic-size halos which preferentially occur along the filaments, at the intersection of which large halos form that will host galaxy clusters. Finally, the bottom-left panel zooms into an individual galactic-size halo. This shows a large number of small substructures that survive the collapse of the halo and make up about 10% of the total mass (Klypin et al. 1999, Moore et al. 1999)

For simulations like the ones illustrated in Figure 1, it is possible to characterize the statistical properties of the dark matter distribution with very high accuracy. For example, Figure 2 shows the 2-point correlation function, $`\xi (r)`$, of the dark matter (a measure of its clustering strength) in the simulation depicted in the top-right of Figure 1 (Jenkins et al. 1998). The statistical error bars in this estimate are actually smaller than the thickness of the line. Similarly, higher order clustering statistics, topological measures, the mass function and clustering of dark matter halos and the time evolution of these quantities can all be determined very precisely from these simulations (e.g. Jenkins et al. 2001, Evrard et al. 2002). In a sense, the problem of the distribution of dark matter in the $`\mathrm{\Lambda }`$CDM model can be regarded as largely solved444However, the innermost structure of halos like those in the bottom-left of Figure 1 is still a matter of controversy..

In contrast to the clustering of the dark matter, the process of galaxy formation is still poorly understood. How then can dark matter simulations like those of Figure 1 be compared with observational data which, for the most part, refer to galaxies? On large scales a very important simplification applies: for Gaussian theories like CDM, it can be shown that if galaxy formation is a local process, that is, if it depends only upon local physical conditions (density, temperature, etc), then, on scales much larger than that associated with individual galaxies, the galaxies must trace the mass, i.e. on sufficiently large scales, $`\xi _{gal}(r)\xi _{dm}(r)`$ (Coles 1993). It suffices therefore to identify a random subset of the dark matter particles in the simulation to obtain an accurate prediction for the properties of galaxy clustering on large scales. This idea (complemented on small scales by an empirical prescription in the manner described by Cole et al. 1998) has been used to construct the mock versions of a region of the APM galaxy survey and of a slice of the 2dFGRS displayed in Figures 3 and 4 which also show the real data for comparison in each case. By eye at least, it is very difficult to distinguish the mocks from the real data.

A quantitative comparison between simulations and the real world is carried out in Figure 5. The symbols show the estimate of the power spectrum in the 2dFGRS survey (Percival et al. 2001). This is the raw power spectrum convolved with the survey window function and can be compared directly with the line showing the theoretical prediction obtained from the mock catalogues which have exactly the same window function. The agreement between the data and the $`\mathrm{\Lambda }`$CDM model is remarkably good.

### 3.2 Galaxy formation

Understanding galaxy formation is a much more difficult problem than understanding the evolution of the dark matter distribution. In the CDM theory, galaxies form when gas, initially well mixed with the dark matter, cools and condenses into emerging dark matter halos. In addition to gravity, a non-exhaustive list of the processes that now need to be taken into account includes: the shock heating and cooling of gas into dark halos, the formation of stars from cold gas and the evolution of the resulting stellar population, the feedback processes generated by the ejection of mass and energy from evolving stars, the production and mixing of heavy elements, the extinction and reradiation of stellar light by dust particles, the formation of black holes at the centres of galaxies and the influence of the associated quasar emission. These processes span an enormous range of length and mass scales. For example, the parsec scale relevant to star formation is a factor of $`10^8`$ smaller than the scale of a galaxy supercluster.

The best that can be done with current computing techniques is to model the evolution of dark matter and gas in a cosmological volume with resolution comparable to a single galaxy. Subgalactic scales must then then be regarded as “subgrid” scales and followed by means of phenomenological models based either on our current physical understanding or on observations. In the approach known as “semi-analytic” modelling (White & Frenk 1991), even the gas dynamics is treated phenomenologically using a simple, spherically symmetric model to describe the accretion and cooling of gas into dark matter halos. It turns out that this simple model works suprisingly well as judged by the good agreement with results of full N-body/gas-dynamical simulations (Benson et al. 2001b, Helly et al. 2002, Yoshida et al. 2002).

The main difficulty encountered in cosmological gas dynamical simulations arises from the need to suppress a cooling instability present in hierarchical clustering models like CDM. The building blocks of galaxies are small clumps that condense at early times. The gas that cools within them has very high density, reflecting the mean density of the Universe at that epoch. Since the cooling rate is proportional to the square of the gas density, in the absence of heat sources, most of the gas would cool in the highest levels of the mass hierarchy leaving no gas to power star formation today or even to provide the hot, X-ray emitting plasma detected in galaxy clusters. Known heat sources are photoionisation by early generations of stars and quasars and the injection of energy from supernovae and active galactic nuclei. These processes, which undoubtedly happened in our Universe, belong to the realm of subgrid physics which cosmological simulations cannot resolve. Different treatments of this “feedback” result in different amounts of cool gas and can lead to very different predictions for the properties of the galaxy population. This is a fundamental problem that afflicts cosmological simulations even when they are complemented by the inclusion of semi-analytic techniques. In this case, the resolution of the calculation can be extended to arbitrarily small mass halos, perhaps allowing a more realistic treatment of feedback. Although they are less general than full gasdynamical simulations, simulations in which the evolution of gas is treated semi-analytically make experimentation with different prescriptions relatively simple and efficient (Kauffmann White & Guiderdoni 1993, Somerville & Primack 1999, Cole et al. 2000)

The outcome of an N-body dark matter simulation in a $`\mathrm{\Lambda }`$CDM universe in which the visible properties of the galaxies have been calculated using the semi-analytic model of Cole et al. (2000) is illustrated in Fig. 6 (Benson et al. 2001a). Galaxies form mostly along the filaments delineated by the dark matter. Red galaxies predominate in the most massive dark matter halos, just as observed in real galaxy clusters. This segregation is a natural outcome of hierarchical clustering from CDM initial conditions. It reflects the fact that the progenitors of rich clusters form substantially earlier than a typical dark matter halo of the same mass. Fig. 7 shows the galaxy luminosity function which describes the abundance of galaxies of different luminosities. The theoretical predictions, shown by the line, agree remarkably well with the observations but this should not be regarded as a spectacular success of the theory because the free parameters in the semi-analytic star formation and feedback model have been tuned to achieve as good a match as possible to this specific observational dataset. In particular, the feedback model has been tuned to produce a relatively flat function at the faint end.

Having fixed the model parameters by reference to a small subset of the data such as the galaxy luminosity function, we can ask whether the same model accounts for other basic observational data. The galaxy autocorrelation function, $`\xi _{gal}(r)`$, in the simulations is plotted in Fig. 2 above. On large scales, it follows $`\xi _{dm}(r)`$ quite closely, but on small scales it dips below the mass autocorrelation function. This small scale “antibias” has also been seen in N-body/gasdynamical simulations of the $`\mathrm{\Lambda }`$CDM cosmology (Pearce et al. 1999, 2001, Dave et al. 1999), and in dark matter simulations that resolve individual galactic halos (Klypin et al. 1999). The galaxy autocorrelation function in the simulations of Benson et al. (2000) agrees remarkably well with the observational data (see also Kauffmann et al. 1999a). This is a genuine success of the theory because no model parameters have been adjusted in this comparison. The differences between the small-scale clustering of galaxies and dark matter result from the interplay between the clustering of dark matter halos and the occupation statistics of galaxies in halos which, in turn, are determined by the physics of galaxy formation. This conclusion, discussed in detail by Benson et al. (2000), has led to the development of an analytic formulation known as the “halo model” (e.g. Seljak 2000, Peacock & Smith 2000, Berlind & Weinberg 2002).

Another genuine prediction of the model is the dependence of the strength of clustering on the luminosity of different subsamples. It can be seen in Fig. 6 that the brightest galaxies are concentrated in the most massive clusters, leading one to suspect that their autocorrelation function must be stronger than average. This is indeed the case, as illustrated in Fig. 8 which compares the variation of the clustering length (defined as the pair separation for which $`\xi (r)=1`$) of galaxy samples of different intrinsic luminosity in the simulations of Benson et al. (2001a) with the observational data obtained from the 2dFGRS by Norberg et al. (2001a). The agreement between theory and observations is remarkable considering that there are no adjustble parameters in this comparison. The reason for the strong clustering of bright galaxies is related to the colour-density relation seen in Fig. 6: the brightest galaxies form in the highest peaks of the density distribution which, in initially Gaussian fields, are more strongly clustered than average peaks which produce less extreme galaxies.

The patch of model universe illustrated in the top panel of Fig. 6 is shown at the earlier epoch corresponding to redshift $`z=3`$ (when the universe was only about 20% of its current age) in the bottom panel of this figure. The galaxies are now blue, reflecting the colour of their younger stellar population. There are fewer galaxies in this plot than in the $`z=0`$ slice. In fact, this is the epoch when the first substantial population of bright galaxies formed in the simulation. As Baugh et al. (1998) argued, the properties of these model galaxies resemble those of the “Lyman-break” galaxies discovered by Steidel et al. (1996), even though different models make somewhat different predictions for their exact properties (Somerville et al. 2001). Most models, however, predict that the brightest galaxies at $`z=3`$ should be strongly clustered (Kauffmann et al. 1999b) and, indeed, the models of Baugh et al. (1998) correctly anticipated that the Lyman-break galaxies would have a clustering length comparable to that of bright galaxies today (Adelberger et al. 1998). This too should be regarded as a significant success of this kind of modelling in the $`\mathrm{\Lambda }`$CDM cosmology. As Fig. 6 shows, in contrast to the galaxies, the dark matter is much more weakly clustered at $`z=3`$ than at $`z=0`$, indicating that galaxies were strongly biased at birth.

## 4 Conclusions

Unlike most computational problems in many areas of science, the cosmological problem is blessed with known, well-specified initial conditions. Within a general class of models, it is possible to calculate the properties of primordial perturbations in the cosmic energy density generated by quantum processes during an early inflationary epoch. In a wide family of inflationary models, these perturbations are adiabatic, scale-invariant and have Gaussian-distributed Fourier amplitudes. The model also requires an assumption about the nature of the dark matter and the possibilities have now been narrowed down to non-baryonic candidates of which cold dark matter particles seem the most promising. An empirical test of the initial conditions for the formation of structure predicted by the model is provided by the cosmic microwave background radiation. The tiny temperature fluctuations it exhibits have exactly the properties expected in the model. Furthermore, the CMB data can be used to fix some of the key model parameters such as $`\mathrm{\Omega }`$ and $`\mathrm{\Omega }_b`$, while these data, combined with other recent datasets such as the 2dFGRS, allow the determination of many of the remaining parameters such as $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and $`h`$. It this specificity of the cosmological problem that has turned simulations into the primary tool for connecting cosmological theory to astronomical observations.

In addition to well-specified initial conditions, the cosmological dark matter problem has the advantage that the only physical interaction that is important is gravity. The problem can thus be posed as a gravitational N-body problem and approached using the many sophisticated techniques that have been developed over the past two decades to tackle this problem. Although on small scales there remain a number of unresolved issues, it is fair to say that on scales larger than a few megaparsecs, the distribution of dark matter in CDM models is essentially understood. The inner structure of dark matter halos, on the other hand, is still a matter of debate and the mass function of dark matter halos has only been reliably established by simulations down to masses of order $`10^{11}\mathrm{M}_{}`$. Resolving these outstanding issues is certainly within reach, but this will require carefully designed simulations and large amounts of computing power.

The frontier of the subject at present lies in simulations of the formation, evolution and structure of galaxies. This problem requires first of all a treatment of gas dynamics in a cosmological context and a number of techniques, relying on direct simulations or on semi-analytical approximations, are being explored. There are quite a few different approaches to cosmological gasdynamics, but it is reassuring that they all give similar results in the simplest relevant problem, the evolution of non-radiative gas during the formation of a galaxy cluster. No detailed comparisons exist yet for the more complicated case in which the gas is allowed to cool, but at least one of the gasdynamic simulation techniques, SPH, gives quite similar results to a simple semi-analytic approach. Realistic models of galaxy formation, however, will require much more than a correct treatment of cooling gas. Such models will necessarily have to include a plethora of astrophysical phenomena such as star formation, feedback, metal enrichment, etc. The huge disparity between the submegaparsec scales on which these processes operate and the gigaparsec scale of the large-scale structure makes it impossible to contemplate a comprehensive ab initio calculation. The way forward is clearly through a hybrid approach combining direct simulation of processes operating on a limited range of scales with a phenomenological treatment of the others. There is currently a great deal of activity in the phenomenology of galaxy formation.

In spite of the uncertainties that remain, all the indications are that our Universe is well described by a model in which

* the overall geometry is flat;
* the dominant dynamical components are cold dark matter ($`30\%`$) and dark energy ($`70\%`$) with baryons playing very much a supporting role ($`4\%`$);
* the initial conditions are quantum fluctuations in the primordial energy density generated during inflation and
* structure has grown primarily as a result of the gravitational instability experienced by mass fluctuations in an expanding universe.

A skeptic is entitled to feel that the current paradigm is odd, to say the least. Not only is there a need to invoke vast amounts of as yet undetected non-baryonic cold dark matter, but there is also the need to account for the dominant presence of a dark energy whose very existence is a mystery within conventional models of fundamental physics. Odd as it may seem, however, this model accounts remarkably well for a large and diverse collection of empirical facts that span 13 billion years of evolution.

## 5 Acknowledgements

I am grateful to my collaborators for their contribution to the work reviewed here, especially Carlton Baugh, Andrew Benson, Shaun Cole, Adrian Jenkins, Cedric Lacey, Peder Norberg, John Peacock, Will Percival, and Simon White.

## 6 References

Adelberger K. L., Steidel C. C., Giavalisco M., Dickinson M., Pettini M., Kellogg M. (1998), ApJ, 505, 18

Bardeen, J. M., Bond, J. R., Kaiser, N. and Szalay, A. S. (1986), Ap.J., 304, 15

Barnes, J., Hut, P. (1986), Nature, 324, 446

Baugh C. M. (1996), MNRAS, 280, 26

Baugh C. M., Cole S., Frenk C. S. (1996b), MNRAS, 283, 1361

Baugh C. M., Cole S., Frenk C. S. Lacey C. G., (1998), ApJ, 498, 504

Benson A. J., Cole S., Frenk C. S., Baugh C. M., Lacey C. G. (2000), MNRAS, 311, 79

Benson, A. J., Frenk, C. S., Baugh, C. M., Cole, S., Lacey, C. G. (2001a), MNRAS, 327, 1041

Benson, A. J., Pearce, F. R., Frenk, C. S., Baugh, C. M. & Jenkins, A. R. (2001b), MNRAS, 320, 261

Berlind, A. & Weinberg, D. (2002), ApJ, in press (astro-ph/0109001)

de Bernardis, P. etal (2000), Nature 404, 995

Borgani, S. etal (2001), ApJ, 561, 13

Cen, R. (1992), ApJ Suppl 78, 341

Clowe, D., Luppino, G. A., Kaiser, N., Gioia, I. M. (2000), ApJ, 539, 540

Cole, S., Hatton, S., Weinberg, D. H., Frenk, C. S. (1998), MNRAS, 300, 945

Cole, S., Lacey, C. G., Baugh, C. M., Frenk, C.S. (2000), MNRAS, 319, 168

Coles, P. (1993) MNRAS, 262, 1065

Colless, M. etal (the 2dFGRS team) (2001), MNRAS, 328, 1039

Couchman, H. M. P., Thomas, P. A., Pearce, F. R. (1995), ApJ, 452, 797

Dave, R., Hernquist, L., Katz, N. & Weinberg, D. (1999), in Proceedings of Rencontres Internationales de l’IGRAP, Clustering at High Redshift, Marseille 1999, astro-ph/9910221

Davis, M., Efstathiou, G., Frenk, C. S. & White, S. D. M. (1985), ApJ., 292, 371

Efstathiou, G., Davis, M., Frenk, C. S. & White, S. D. M. (1985), ApJ Suppl, 57, 241

Efstathiou, G. etal (The 2dFGRS team) (2002), MNRAS, in press, (astro-ph/0109152)

Efstathiou, G., Sutherland, W. J., Maddox, S. J. (1990), Nature, 348, 705

Eke, V. R., Cole, S. & Frenk, C. S. (1996), MNRAS, 282, 263

Eke, V. R., Cole, S. Frenk, C. S. & Henry, J. P. (1998), MNRAS, 298, 1145

Evrard, A. E. (1997), MNRAS, 292, 289.

Evrard, A.E., MacFarland, T., Couchman, H.M.P., Colberg, J.M., Yoshida, N., White, S.D.M., Jenkins, A., Frenk, C.S., Pearce, F.R., Efstathiou, G., Peacock, J. & Thomas, P., (The Virgo Consortium) (2002), ApJ, in press.

Fischer, P. etal (2000), AJ 120, 1198

Frenk, C. S. etal (1999), ApJ, 525, 554

Freedman, W. L. etal (2001), ApJ, 553, 47

Gingold, R. A. & Monaghan, J. J. (1977), MNRAS 181, 375

Gnedin, N. Y. (1995), ApJ Suppl, 97, 231

Guth, A. H. (1981), Phys Rev D, 23, 347

Hanany, S. etal (2000), ApJ Lett, 545, 5

Helly, J., Cole, S. M., Frenk, C. S., Baugh, C., Benson, A. & Lacey, C. G. (2002), submitted to MNRAS (astro-ph/0202485)

Jenkins, A., Frenk, C. S., Pearce, F. R., Thomas, P., Colberg, J., White, S. D. M., Couchman, H., Peacock, J., Efstathiou, G. & Nelson, A. (1998), ApJ, 499, 20

Jenkins A., Frenk C. S., White S. D. M., Colberg J. M., Cole S., Evrard A. E., Yoshida N. (2001), MNRAS, 321,372

Kaiser, N. (1987), MNRAS, 227, 1

Kauffmann, G., Colberg, J., Diaferio, A., White, S. D. M. (1999a), MNRAS, 303, 188

Kauffmann, G., Colberg, J., Diaferio, A., White, S. D. M. (1999b), MNRAS, 307, 529

Kauffmann G., White S. D. M., Guiderdoni B. (1993), MNRAS, 264, 201

Klypin, A., Gottlober, S., Kravstov, A., Khokhlov, A.M. (1999), ApJ, 516, 530

Leitch, E.M. etal (2002), ApJ, 568, 28L

Lifshitz, E. M. (1946), J Phys, 10, 116

Loveday, J., Peterson, B., Efstathiou, G. & Maddox, S. (1992), ApJ, 319, 338

McKay, T. etal (2001) ApJ, in press (astro-ph/0108013)

Moore, B., Ghigna, S., Governato, F., Lake, G., Quinn, T., Stadel, J., Tozzi, P. (1999), ApJ, 524, 19

Navarro, J, F., Power, C., Jenkins, A. R., Frenk, C. S. & White, S. D. M. (2002), in preparation

Norberg, P., Baugh, C. M., Hawkings, E., Maddox, S., Peacock, J,. Cole, S., Frenk, C. S., etal (The 2dFGRS team) (2001a), MNRAS, 328, 64

Norberg, P., Cole, S., Baugh, C. M., Frenk, C. S. etal (The 2dFGRS team) (2001b), submitted to MNRAS, (Astro-ph/0111011)

Peacock, J. A. & Smith, R. (2000) 318, 1144

Peacock, J. etal (The 2dFGRS team) (2001), Nature, 410, 169

Pearce, F. R., Jenkins, A. R., Frenk, C. S., Colberg, J., White, S. D. M., Thomas, P., Couchman, H. M. P., Peacock, J. A., Efstathiou, G. (the Virgo Consortium) (1999), ApJ, 521, 992

Pearce, F. R., Jenkins, A. R., Frenk, C. S., White, S. D. M., Thomas, P., Couchman, H. M. P., Peacock, J. A., Efstathiou, G. (the Virgo Consortium) (2001), MNRAS, 326, 649

Peebles, P. J. E. (1970), AJ, 75, 13

Peebles, P. J. E. (1980), The large scale structure of the Universe. (Princeton: Princeton University Press).

Pen, U-L. (1998), ApJ Suppl, 115, 19

Percival, W. etal (the 2dFGRS team) 2001, MNRAS, 327, 1297

Perlmutter, S. etal (The Supernova Cosmology Project) (1999), ApJ. 517, 565

Pierpaoli, E., Scott, D. & White, M. (2001), MNRAS 325, 77

Riess, A. etal (1998), AJ, 116, 1009

Seljak, U. (2000), MNRAS, 318, 203

Smoot, G. etal (1992), Apj 396, 1

Somerville R. S., Primack J. R. (1999), MNRAS, 310, 1087

Somerville, Rachel S., Primack, Joel R., Faber, S. M. (2001), MNRAS, 320, 504

Springel, V., Yoshida, N. & White, S. D. M. (2001), New Astronomy, 6, 79

Stadel, J. (2000), PhD thesis, University of Washington, Seattle

Starobinskii, A. (1982), Phys. Lett. 117B, 175

Steidel, C. C., Giavalisco, M., Pettini, M., Dickinson, M., & Adelberger, K. L. (1996), ApJ, 462, L17

Tytler, D., O’Meara, J. M., Suzuki, N. & Lubin, D. (2000), Physica Scripta, 85, 12.

Van Waerbeke, L. etal (2001), A&A, 374, 757.

Viana, P. T. P. & Liddle, A. (1996), MNRAS, 281, 323

Wilson, G., Kaiser, N., Luppino, G., & Cowie, L. L. (2001), ApJ, 555, 572

White, S. D. M. & Frenk, C. S. (1991), ApJ, 379, 52

White, S. D. M., Frenk, C. S. & Davis, M. (1983) ApJ 274, L1.

White, S. D. M., Navarro, J. F., Evrard, A. E., Frenk, C. S. (1993), Nature, 366, 429.

York, D.G, etal (the SDSS collaboration) (2000) AJ, 120, 1579

Yoshida, N., Stoehr, F., Springel, V. & White, S. D. M. (2002), submitted to MNRAS (Astro-ph/0202341)

Zucca, E. etal (1997), AA, 326, 477.


format: text

text:
# 1 Introduction

## 1 Introduction

The existence of plane-fronted parallel wave solutions in general relativity has been known since a long time -. These solutions are not only interesting by themselves, but also because similar solutions accompany massless fields of lower spin, such as electro-magnetic, scalar and spinor fields , -.

Gravitational $`pp`$-waves can be described by space-time metrics

$$g_{\mu \nu }dx^\mu dx^\nu =dudvK(u,x,y)du^2+dx^2+dy^2,$$
(1)
where the light-cone co-ordinates $`(u,v)`$, transverse to the $`x`$-$`y`$-plane, are related to time- and longitudinal co-ordinates $`(t,z)`$ by:

$$u=ctz,v=ct+z.$$
(2)
There are similar solutions with the roles of $`v`$ and $`u`$ interchanged, which we will not discuss explicitly here. With the metric (1), the only non-vanishing connection co-efficients are:

$$\mathrm{\Gamma }_{uu}^v=K_{,u},\mathrm{\Gamma }_{uu}^x=\frac{1}{2}\mathrm{\Gamma }_{xu}^v=\frac{1}{2}K_{,x},\mathrm{\Gamma }_{uu}^y=\frac{1}{2}\mathrm{\Gamma }_{yu}^v=\frac{1}{2}K_{,y}.$$
(3)
The free gravitational $`pp`$-waves are solutions of the free Einstein equations, which for the metrics (1) simplify to

$$\mathrm{\Delta }_{trans}K=K_{,xx}+K_{,yy}=0.$$
(4)
The simplest non-trivial one is

$$K(u,x,y)=\frac{1}{2}\kappa _+(u)(x^2y^2)+\kappa _\times (u)xy.$$
(5)
Here $`\kappa _{(+,\times )}(u)`$ are the amplitudes of the two different polarization states, as appropriate to quadrupole-type waves propagating at the speed of light.

Metrics of the type (1) are also found in the solution of the coupled Maxwell-Einstein, massless Klein-Gordon-Einstein and massless Dirac-Einstein equations. As a physically interesting example, consider the gravitational field of a light wave propagating along the $`z`$-axis. The light-wave is characterized completely by a vector potential 1-form

$$𝐀=\mathrm{sin}ku(a_xdx+a_ydy),$$
(6)
with $`(a_x,a_y)`$ the components of the constant transverse polarization vector. The corresponding electric and magnetic fields are

$$𝐄(u)=\omega 𝐚\mathrm{cos}ku,𝐁(u)=𝐤\times 𝐚\mathrm{cos}ku,$$
(7)
where $`\omega =kc`$. These fields provide a solution of the coupled Maxwell-Einstein with a metric of the form (1), with

$$\begin{array}{ccc}K(u,x,y)\hfill & =\hfill & \frac{2\pi \epsilon _0G}{c^4}\left(𝐄^2(u)+c^2𝐁^2(u)\right)\left(x^2+y^2\right)\hfill \\ & & \\ & =\hfill & \frac{2\pi \epsilon _0G}{c^2}k^2𝐚^2(1+\mathrm{cos}2ku)\left(x^2+y^2\right).\hfill \end{array}$$
(8)
The motion of a particle of mass $`m`$ and charge $`q`$ propagating in gravitational and electro-magnetic background fields is governed, in the limit of neglecting radiation reaction terms, by the covariant Lorentz force

$$\ddot{x}^\mu +\mathrm{\Gamma }_{\lambda \nu }^\mu \dot{x}^\lambda \dot{x}^\nu =\frac{q}{m}F_\nu ^\mu \dot{x}^\nu .$$
(9)
Here the overdot denotes a derivative w.r.t. proper time, defined in the usual way from the universal conserved value of the total four-velocity:

$$u_\mu ^2=g_{\mu \nu }\frac{dx^\mu }{d\tau }\frac{dx^\mu }{d\tau }=c^2.$$
(10)
We apply these equations to study the motion of a particle in the background of the waves fields (1,6,7). Recalling the results of ref., the light-cone co-ordinate $`u`$ can be used as a proper-time variable via the relation

$$\dot{u}=\gamma =\text{const.}u=\gamma \tau .$$
(11)
The expression (10) for the total four-velocity then leads to the following relation between proper time and laboratory time:

$$\frac{dt}{d\tau }=\gamma \frac{dt}{du}=\sqrt{\frac{1\gamma ^2K/c^2}{1𝐯^2/c^2}},$$
(12)
showing that there is both a gravitational and a kinematic time dilation. Substitution back into eq.(10) allows one to cast it into a true conservation law of the form

$$hK+\frac{1𝐯^2/c^2}{(1v_z/c)^2}=\frac{c^2}{\gamma ^2}.$$
(13)
Two independent equations of motion remain, which we take to be those for the tranverse co-ordinates $`𝝃=(x,y)`$. The Lorentz-Einstein equation (9) for these components becomes

$$\frac{d^2𝝃}{du^2}=\frac{2\pi \epsilon _0G}{c^2}k^2𝐚^2\left(1+\mathrm{cos}2ku\right)𝝃\frac{q}{m\gamma }k𝐚\mathrm{cos}ku.$$
(14)
We can define the components w.r.t. the polarization of the light wave as

$$𝝃=𝝃_{}+𝝃_{},𝝃_{}=\frac{𝝃𝐚}{|𝐚|^2}𝐚,𝝃_{}=\frac{𝝃\times 𝐚}{|𝐚|}.$$
(15)
Under this decomposition eq.(14) splits into a homogeneous and an inhomogeneous equation

$$\begin{array}{c}\frac{d^2𝝃_{}}{du^2}+\frac{2\pi \epsilon _0G}{c^2}k^2𝐚^2\left(1+\mathrm{cos}2ku\right)𝝃_{}=0,\hfill \\ \\ \frac{d^2𝝃_{}}{du^2}+\frac{2\pi \epsilon _0G}{c^2}k^2𝐚^2\left(1+\mathrm{cos}2ku\right)𝝃_{}=\frac{q}{m\gamma }k𝐚\mathrm{cos}ku.\hfill \end{array}$$
(16)
In refs.  it was observed, that in the static limit $`k0`$ taken such that $`E_0=cB_0=ck|𝐚|=`$ constant, one obtains a simple harmonic motion with frequency $`\nu `$ (in Hz) given by:

$$2\pi \nu =\sqrt{\frac{4\pi \epsilon _0G}{c^2}}E_0=0.3\times 10^{18}E_0\text{(V/m)}.$$
(17)
In this paper we discuss the motion in the non-static case $`k0`$.

## 2 Transverse motion

We first focus on the case $`q=0`$, describing a neutral massive test particle in an electro-magnetic and gravitational $`pp`$-wave. According to (16), the transverse coordinates $`(x,y)`$ all obey the same homogeneous Mathieu equation. Besides, we choose to restrict our study to the case where the only contribution to the metric co-efficient $`K`$ is the electromagnetic one, given by (8). In other words, no free gravitational $`pp`$-waves are taken to be present: $`\kappa _{+,\times }(u)`$=0. It then follows that the motion of the particle is invariant under the exchange of the tranverse coordinates and rotations in the transverse plane. The equation for $`\xi =(x,y)`$ now reads (with $`a^2=𝐚^2`$):

$$\frac{d^2\xi }{du^2}+\frac{2\pi \epsilon _0G}{c^2}k^2a^2\left(1+\mathrm{cos}2ku\right)\xi =0,$$
(18)
A solution procedure for eq.(18) in the weak-coupling regime was presented in ref. . Here we discuss the construction of the full solution in terms of a generalized Fourier series expansion. Because of the particular form of the equation (18), we can construct a series in terms of either even or odd multiples of the fundamental wave factor $`k`$. In this section we consider the generalized odd solutions, which can be parametrized as

$$\xi (u)=e^{\sigma ku}\underset{n=0}{\overset{\mathrm{}}{}}\left(\xi _n\mathrm{cos}(2n+1)ku+\eta _n\mathrm{sin}(2n+1)ku\right),$$
(19)
with real Fourier components $`(\xi _n,\eta _n)`$. Note that for $`\sigma ^20`$ (and taking the real part) one obtains solutions which are purely of trigonometric type; for $`\sigma ^2>0`$ the solutions describe oscillations enhanced by parametric resonance. The other type of solutions having even Fourier part is discussed in the appendix. The substitution (19) converts the Mathieu equation into an infinite-dimensional linear matrix equation

$$\left(\begin{array}{cccc}A_0& B& 0& \mathrm{}\\ B& A_1& B& \mathrm{}\\ 0& B& A_2& \mathrm{}\\ .& .& .& .\end{array}\right)\left(\begin{array}{c}\rho _0\\ \rho _1\\ \rho _2\\ .\end{array}\right)=0,$$
(20)
where the blocks $`A_m`$ and $`B`$ are $`2\times 2`$ matrices given by

$$\begin{array}{ccc}\hfill A_0& =& \left(\begin{array}{cc}\sigma ^21+\frac{3\kappa ^2}{4}& 2\sigma \\ 2\sigma & \sigma ^21+\frac{\kappa ^2}{4}\end{array}\right),\hfill \\ & & \\ \hfill A_m& =& \left(\begin{array}{cc}\sigma ^2(2m+1)^2+\frac{\kappa ^2}{2}& 2(2m+1)\sigma \\ 2(2m+1)\sigma & \sigma ^2(2m+1)^2+\frac{\kappa ^2}{2}\end{array}\right),(m1),\hfill \end{array}$$
(21)
and

$$B=\frac{\kappa ^2}{4}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\rho _m=\left(\begin{array}{c}\xi _m\\ \eta _m\end{array}\right).$$
(22)
A solution of the system (20) exists provided the linear matrix operator has zero-modes. To determine their existence we cast the equation into upper-triangular form:

$$\left(\begin{array}{cccc}a_0& B& 0& \mathrm{}\\ 0& a_1& B& \mathrm{}\\ 0& 0& a_2& \mathrm{}\\ .& .& .& .\end{array}\right)\left(\begin{array}{c}\overline{\rho }_0\\ \overline{\rho }_1\\ \overline{\rho }_2\\ .\end{array}\right)=0.$$
(23)
The $`2\times 2`$ matrices $`a_m`$ can be expressed in terms of continuing fractions as

$$a_m=A_mB^2\left[A_{m+1}B^2\left[A_{m+2}\mathrm{}\right]^1\right]^1,$$
(24)
and the basis of vectors has been transformed by $`\overline{\rho }_0=\rho _0`$, and

$$\overline{\rho }_m=\rho _m+a_m^1B\rho _{m1},(m1).$$
(25)
This transformation is singular if and only if any of the $`a_m`$ $`(m1)`$ has zero-modes, which is precisely the condition for the Mathieu equation to have non-trivial solutions.

Zero-modes of $`a_0`$ represent the lowest-frequency solutions of the Mathieu equation; they exhibit parametric resonance for $`\sigma >0`$. The existence of zer-modes implies that the determinant of $`a_0`$ vanishes. This determinant can be computed in successive approximations by truncating the infinite fraction (24). The zeroth-order approximation is

$$a_0=A_0\left|a_0^{(0)}\right|=(\sigma ^2+1)^2+\kappa ^2(\sigma ^21)+\frac{3\kappa ^4}{16}.$$
(26)
The first- and higher-order approximations are of the form

$$a_0=\left(\begin{array}{cc}\sigma ^21+\frac{3\kappa ^2}{4}+\kappa ^4\mathrm{\Delta }_{11}& 2\sigma +\kappa ^4\mathrm{\Delta }_{12}\\ (2\sigma +\kappa ^4\mathrm{\Delta }_{12})& \sigma ^21+\frac{\kappa ^2}{4}+\kappa ^4\mathrm{\Delta }_{11}\end{array}\right)+𝒪(\kappa ^6),$$
(27)
where the quantities $`\mathrm{\Delta }_{ij}(\sigma )`$ are rational functions of $`\sigma ^2`$; the determinant then is

$$|a_0|=(\sigma ^2+1)^2+\kappa ^2(\sigma ^21)+\kappa ^4\left(\frac{3}{16}\frac{\mathrm{\Delta }}{8}\right)+𝒪(\kappa ^6).$$
(28)
In first order $`a_0^{(1)}=A_0B^2A_1^1`$, with

$$\mathrm{\Delta }_{11}^{(1)}=\frac{1}{48}\frac{9\sigma ^2}{27+6\sigma ^2+\sigma ^4/3},\mathrm{\Delta }_{12}^{(1)}=\frac{1}{8}\frac{1}{27+6\sigma ^2+\sigma ^4/3},$$
(29)
and

$$\mathrm{\Delta }^{(1)}=\frac{\sigma ^422\sigma ^2+9}{\sigma ^4+18\sigma ^2+81}.$$
(30)
Similarly, in second order with

$$a_0=\left(A_0A_1A_2B^2(A_0+A_2)\right)\left(A_1A_2B^2\right)^1,$$
(31)
one finds up to terms $`𝒪(\kappa ^6)`$:

$$\begin{array}{ccc}\hfill 16(\sigma ^215)\mathrm{\Delta }_{11}^{(2)}& =& \frac{5}{8}+\frac{3}{8}\frac{(\sigma ^2+25)(\sigma ^494\sigma ^2+225)}{(\sigma ^494\sigma ^2+225)^2+256\sigma ^2(\sigma ^215)^2},\hfill \\ & & \\ \hfill \mathrm{\Delta }_{12}^{(2)}& =& \frac{3}{8}\frac{(\sigma ^2+25)^2}{(\sigma ^494\sigma ^2+225)^2+256\sigma ^2(\sigma ^215)^2}.\hfill \end{array}$$
(32)
Then in the determinant

$$\mathrm{\Delta }^{(2)}=\frac{(\sigma ^21)[(\sigma ^225)(\sigma ^494\sigma ^2+225)+160\sigma ^2(\sigma ^215)]}{(\sigma ^494\sigma ^2+225)^2+256\sigma ^2(\sigma ^215)^2}.$$
(33)
Observe that in both first and second order approximation

$$\underset{\sigma ^20}{lim}\mathrm{\Delta }^{(k)}=\frac{1}{9},\underset{\sigma ^2\mathrm{}}{lim}\mathrm{\Delta }^{(k)}=1,k=(1,2).$$
(34)
From the above analysis it follows, that up to terms of order $`\kappa ^6`$ a solution for $`\sigma ^2`$ exists for which the determinant vanishes: $`|a_0|=0`$, as given by the expression (28). This condition has approximate solutions

$$\sigma ^2=\left(1+\frac{\kappa ^2}{2}\right)\pm \sqrt{2\kappa ^2\left(1+\frac{\kappa ^2}{32}(1+2\mathrm{\Delta })\right)}.$$
(35)
Positive non-zero values of $`\sigma ^2`$ exist only for the positive square root, provided

$$2\kappa ^2\left(1+\frac{\kappa ^2}{32}(1+2\mathrm{\Delta })\right)>\left(1+\frac{\kappa ^2}{2}\right)^2.$$
(36)
or

$$\left(1\frac{\kappa ^2}{2}+\frac{\kappa ^2}{4}\sqrt{1+2\mathrm{\Delta }}\right)\left(1\frac{\kappa ^2}{2}\frac{\kappa ^2}{4}\sqrt{1+2\mathrm{\Delta }}\right)<0.$$
(37)
Hence $`\kappa ^2`$ must be in the range

$$\frac{1}{2+\sqrt{1+2\mathrm{\Delta }}}<\frac{\kappa ^2}{4}<\frac{1}{2\sqrt{1+2\mathrm{\Delta }}}.$$
(38)
Purely periodic solutions of the Mathieu equation exist for a discrete set of values of $`\kappa `$ for which $`\sigma =0`$; for these values the inequality (36) turns into and equality. For all other values outside the range (38) the solutions for $`\sigma `$ are imaginary; the series (19) then is purely trigonometric, but for non-rational $`\sigma `$ the result is not periodic.

From the inequality (38) it follows, that the smallest possible value of $`\kappa ^2`$ for which parametric resonance occurs is of the order unity: $`\kappa ^21`$. As $`\kappa ^2`$ represents the ratio of field intensity and frequency, it follows that one needs very strong fields at very low frequencies $`(2\pi \nu =ck)`$:

$$\kappa ^2=\left(\frac{0.3\times 10^{18}E_0(\text{V/m})}{2\pi \nu (\text{Hz})}\right)^21.$$
(39)
An upper limit to the fields strength of the order of $`10^{19}`$ V/m arises because of pair creation processes which destroy the electric field at higher intenties . Therefore we obtain an upper limit of the order of $`1`$ Hz for the frequencies which can give rise to parametric resonance. For lower fields and higher frequencies one finds standard periodic behavior described by solutions with $`\sigma =0`$ or imaginary.

An important characteristic of the motion of a chargeless particle in the background of the fields (1,6,7) is that a particle initially at the origin in the transverse plane before the arrival of the wave, remains at the origin at all later times. This allows one to translate directly the geodesic motion of a test particle in terms of the relative rate of acceleration of this particle w.r.t. another one located in the origin: $`𝝃`$ = 0.

For charged particles this situation changes. If $`q0`$ and we take the transverse co-ordinate $`\xi =\xi _{}`$ parallel to the electric field, the equation of transverse motion is:

$$\frac{d^2\xi }{du^2}+\frac{\kappa ^2k^2}{2}\left(1+\mathrm{cos}2ku\right)\xi =\epsilon \mathrm{cos}ku,$$
(40)
with

$$\epsilon =\frac{qka}{m\gamma }=\frac{qE_0}{mc\gamma }.$$
(41)
In the inhomogeneous case the absence of an expontential term in the driving force implies that any special solution of the inhomogeneous equation can have an expansion (19) only with $`\sigma =0`$; for these special solutions the inhomogeneous Mathieu equation becomes

$$\left(\begin{array}{cccc}a_0& B& 0& \mathrm{}\\ 0& a_1& B& \mathrm{}\\ 0& 0& a_2& \mathrm{}\\ .& .& .& .\end{array}\right)_{\sigma =0}\left(\begin{array}{c}\overline{\rho }_0\\ \overline{\rho }_1\\ \overline{\rho }_2\\ .\end{array}\right)=\left(\begin{array}{c}\theta \\ 0\\ 0\\ .\end{array}\right),\theta =\frac{\epsilon }{k^2}\left[\begin{array}{c}1\\ 0\end{array}\right].$$
(42)
Here we are interested especially in obtaining a particular solution of the previous type: $`\overline{\rho }_00`$, $`\overline{\rho }_n=0`$, $`n1`$ . The non-zero component $`\overline{\rho }_0`$ is a solution of the inhomogeneous equation

$$a_0\overline{\rho }_0=\theta .$$
(43)
Now $`a_0`$ is of the form (27), with $`\sigma =0`$:

$$a_0=\left(\begin{array}{cc}1+\frac{3\kappa ^2}{4}+\mathrm{\Delta }_{11}(0)\kappa ^4& 0\\ 0& 1+\frac{\kappa ^2}{4}+\mathrm{\Delta }_{11}(0)\kappa ^4\end{array}\right)+𝒪(\kappa ^6),$$
(44)
where beyond the lowest-order approximation we find $`\mathrm{\Delta }_{11}(\sigma =0)=\frac{1}{144}`$. As a result a special solution of the inhomogeneous equation is given by

$$\overline{\rho }_0=\frac{\epsilon }{k^2(1\frac{3\kappa ^2}{4}\mathrm{\Delta }_{11}(0)\kappa ^4+\mathrm{})}\left[\begin{array}{c}1\\ 0\end{array}\right],\overline{\rho }_n=0,n1.$$
(45)
The general solution of this system is the sum of a particular solution plus a solution of the homogeneous equation, which can have $`\sigma 0`$. An important difference with the previous case is that a particle initially at the origin does not remain at the origin, but is accelerated by the electric field. Thus the study of the relative acceleration between particles on different world lines becomes more complicated; for a general discussion we refer to ref. .

## 3 QED corrections

As the magnitude of the electric field in the region of parametric resonance is necessarily extremely high, QED effects such as vacuum polarization can in principle not be ignored. However, as we now show they do not affect the classical wave solutions we have studied in the previous paragraphs.

The non-linear effects of quantum fluctuations on the propagation of light in flat (Minkowski) space-time were first studied by Euler and Heisenberg . The effects are summerized by a well-known effective lagrangean incorporating the corrections due to the electromagnetic higher-order terms; its generally covariant form reads:

$$=\epsilon _0c^2\sqrt{g}\left(\frac{1}{4}F^2+\beta \left[(F^2)^2+\frac{7}{4}\varphi ^2\right]\right)$$
(46)
where:

$$\begin{array}{c}F^2=g^{\mu \kappa }g^{\nu \lambda }F_{\mu \nu }F_{\kappa \lambda },\hfill \\ \\ \varphi =F_{\mu \nu }\stackrel{~}{F}^{\mu \nu },\hfill \end{array}$$
(47)
and $`\stackrel{~}{F}_{\mu \nu }`$ is the dual electromagnetic tensor defined by:

$$\stackrel{~}{F}^{\mu \nu }=\frac{1}{2\sqrt{g}}\epsilon ^{\mu \nu \rho \sigma }F_{\rho \sigma }.$$
(48)
The coupling constant $`\beta `$ is given by

$$\beta =\frac{\alpha ^2\mathrm{}^3\epsilon _0}{90m^4c^4}.$$
(49)
We first derive the Maxwell equations from the Euler-Lagrange equations taking $`A_\nu `$ and $`_\nu A_\mu `$ as the variables:

$$\frac{}{A_\nu }_\mu \frac{}{(_\mu A_\nu )}=0.$$
(50)
Combining (46) and (50), we get the following equations:

$$\begin{array}{c}_\mu \left[(18\beta F^2)F^{\mu \nu }14\beta \varphi \stackrel{~}{F}^{\mu \nu }\right]=0.\hfill \end{array}$$
(51)
The energy-momentum tensor is obtained from the lagrangean (46) by:

$$\begin{array}{ccc}\hfill T_{\mu \nu }& =& \frac{2}{\sqrt{g}}\frac{}{g^{\mu \nu }}\hfill \\ & & \\ & =& \epsilon _0c^2\left[(18\beta F^2)F_{\mu \lambda }F_{\nu }^{}{}_{}{}^{\lambda }\frac{1}{4}g_{\mu \nu }(F^24\beta (F^2)^2)\frac{7\beta }{4}g_{\mu \nu }\varphi ^2\right].\hfill \end{array}$$
(52)
With the metric (1) the Ricci tensor $`R_{\mu \nu }`$ has only one non-zero component $`R_{uu}`$, the Einstein field equations reduce to:

$$R_{uu}=\frac{1}{2}\mathrm{\Delta }_{trans}K=8\pi GT_{uu},\text{all other}T_{\mu \nu }=0.$$
(53)
In our particular case of a plane-wave vector potential of the form (6), the only non-vanishing components of the electromagnetic tensor $`F_{\mu \nu }`$ are:

$$F_{ui}=F_{iu}=a_ik\mathrm{cos}ku,i=(x,y).$$
(54)
It follows that $`F^2=\varphi =0`$, and the Maxwell-Einstein field equations reduce to the original form:

$$\begin{array}{c}_\mu F^{\mu \nu }=0,\hfill \\ \\ T_{uu}=\epsilon _0c^2F_{u\lambda }F_{u}^{}{}_{}{}^{\lambda }.\hfill \end{array}$$
(55)
Thus the quantum corrections do not modify the gravito-electro-magnetic plane wave solution.

Appendix

The solutions of the Mathieu equation are classified as having an odd or even periodic part. Eq.(19) represents the general solution of odd type. The general solution of even type can be expanded as

$$\xi (u)=e^{\sigma ku}\left(\frac{\mu _0}{2}+\underset{n=1}{\overset{\mathrm{}}{}}\left[\mu _n\mathrm{cos}2nku+\nu _n\mathrm{sin}2nku\right]\right).$$
(56)
After multiplication with $`e^{\sigma ku}`$ the homogeneous Mathieu equation becomes

$$\begin{array}{ccc}0\hfill & =\hfill & \left(\frac{1}{2}\sigma ^2+\frac{1}{4}\kappa ^2\right)\mu _0+\frac{\kappa ^2}{4}\mu _1\hfill \\ & & \\ & & +\underset{n=1}{\overset{\mathrm{}}{}}\left[\left(\sigma ^24n^2+\frac{\kappa ^2}{2}\right)\mu _n+4n\sigma \nu _n+\frac{\kappa ^2}{4}\left(\mu _{n1}+\mu _{n+1}\right)\right]\mathrm{cos}2nku\hfill \\ & & \\ & & +\underset{n=1}{\overset{\mathrm{}}{}}\left[\left(\sigma ^24n^2+\frac{\kappa ^2}{2}\right)\nu _n4n\sigma \mu _n\right]\mathrm{sin}2nku\hfill \\ & & \\ & & +\frac{\kappa ^2}{4}\left(\nu _2\mathrm{sin}2ku+\underset{n=2}{\overset{\mathrm{}}{}}(\nu _{n1}+\nu _{n+1})\mathrm{sin}2nku\right).\hfill \end{array}$$
(57)
It follows, that

$$\mu _0=\frac{\kappa ^2\mu _1}{2\sigma ^2+\kappa ^2},$$
(58)
and

$$\left(\begin{array}{cccc}D_1& B& 0& \mathrm{}\\ B& D_2& B& \mathrm{}\\ 0& B& D_3& \mathrm{}\\ .& .& .& .\end{array}\right)\left(\begin{array}{c}\zeta _1\\ \zeta _1\\ \zeta _2\\ .\end{array}\right)=0,$$
(59)
where $`B`$ is the $`2\times 2`$-matrix (22), whilst

$$\begin{array}{c}D_1=\left(\begin{array}{cc}\sigma ^24+\frac{\kappa ^2}{2}\frac{\kappa ^4}{4(2\sigma ^2+\kappa ^2)}& 4\sigma \\ 4\sigma & \sigma ^24+\frac{\kappa ^2}{2}\end{array}\right),\hfill \\ \\ D_n=\left(\begin{array}{cc}\sigma ^24n^2+\frac{\kappa ^2}{2}& 4n\sigma \\ 4n\sigma & \sigma ^24n^2+\frac{\kappa ^2}{2}\end{array}\right),n2,\hfill \end{array}$$
(60)
and

$$\zeta _n=\left(\begin{array}{c}\mu _n\\ \nu _n\end{array}\right).$$
(61)
To solve this equation, we bring the infinite-dimensional matrix equation (59) into upper triangular form:

$$\left(\begin{array}{cccc}D_1& B& 0& \mathrm{}\\ B& D_2& B& \mathrm{}\\ 0& B& D_3& \mathrm{}\\ .& .& .& .\end{array}\right)=\left(\begin{array}{cccc}d_1& B& 0& \mathrm{}\\ 0& d_2& B& \mathrm{}\\ 0& 0& d_3& \mathrm{}\\ .& .& .& .\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& \mathrm{}\\ c_1& 1& 0& \mathrm{}\\ 0& c_2& 1& \mathrm{}\\ .& .& .& .\end{array}\right),$$
(62)
with $`d_n`$ and $`c_n`$ the solutions of the equations

$$d_n+Bc_n=D_n,d_{n+1}c_n=BB^1d_n+(B^1d_{n+1})^1=B^1D_n.$$
(63)
The solution for $`d_n`$ takes the form of a continuing fraction

$$d_n=D_nB^2\left[D_{n+1}B^2\left[D_{n+2}\mathrm{}\right]^1\right]^1,$$
(64)
where we have used the fact that $`B`$ is proportional to the unit matrix, hence commutes with all other matrices. If we define

$$\left(\begin{array}{c}\overline{\zeta }_1\\ \overline{\zeta }_2\\ \overline{\zeta }_3\\ .\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& \mathrm{}\\ c_1& 1& 0& \mathrm{}\\ 0& c_2& 1& \mathrm{}\\ .& .& .& .\end{array}\right)\left(\begin{array}{c}\zeta _1\\ \zeta _2\\ \zeta _3\\ .\end{array}\right),$$
(65)
equivalent with $`\overline{\zeta }_1=\zeta _1`$, and

$$\overline{\zeta }_n=\zeta _n+c_{n1}\zeta _{n1},n2.$$
(66)
Now the equation

$$\left(\begin{array}{cccc}d_1& B& 0& \mathrm{}\\ 0& d_2& B& \mathrm{}\\ 0& 0& d_3& \mathrm{}\\ .& .& .& .\end{array}\right)\left(\begin{array}{c}\overline{\zeta }_1\\ \overline{\zeta }_2\\ \overline{\zeta }_3\\ .\end{array}\right)=0,$$
(67)
has solutions with $`\overline{\zeta }_N0`$, and

$$d_k\overline{\zeta }_k=B\overline{\zeta }_{k+1},1kN1\text{and}\overline{\zeta }_n=0,n>N.$$
(68)
The first one of these is the one given by $`d_1\overline{\zeta }_1=0`$, $`\overline{\zeta }_n=0`$, $`n2`$.

Approximations. Solutions of the equation $`d_1\overline{\zeta }_1=0`$ exist if $`detd_1=0`$. To solve this condition, we need to compute the infinite fraction (64). We can do this by successive approximations:

$$\begin{array}{ccc}d_1^{(0)}\hfill & =\hfill & D_1,\hfill \\ & & \\ d_1^{(1)}\hfill & =\hfill & D_1B^2D_2^1=\left(D_1D_2B^2\right)D_2^1,\hfill \\ & & \\ d_1^{(2)}\hfill & =\hfill & \left(D_1D_2D_3B^2D_1B^2D_3\right)\left(D_2D_3B^2\right)^1,\hfill \\ & & \\ \mathrm{}\hfill & & \end{array}$$
(69)
To zeroth order approximation, we can show that $`detD_1`$ = 0 has solutions for $`8<\kappa ^2<16`$:

$$\left|D_1\right|=\frac{\kappa ^2}{8}\left(\kappa ^28\right)\left(\kappa ^216\right)+2\sigma ^2\left(\sigma ^4+\left(8+\frac{3\kappa ^2}{2}\right)\sigma ^2+16+\frac{5\kappa ^2}{8}\right).$$
(70)
It suffices to note, that the first term is negative in the range $`\kappa ^2(8,16)`$, whilst for any $`\kappa ^2`$ the second term can take all positive values between $`(0,\mathrm{})`$.