Organizations that support science (astronomy) such as federal agencies, research centers, observatories, academic institutions, societies, etc. employ advisory committees and boards as a mechanism for reviewing their activities and giving advice on practices, policies and future directions. As with any scientific endeavor, there is concern over complementing these committees with enough members who have as broad a range of expertise and understanding as possible, so that bias is mitigated. However, for a number of reasons (logistical, practical, financial, etc.), committees can also not be infinitely large and thus trade-offs must be made. It is often recognized that conflicts of interest must be acknowledged within these committees, but what is not often recognized it the potential for unmitigated biases and "group think" that can be introduced as part of these committees. In this white paper, we recommend that advisory committees that collect community input, (e.g., the Decadal Survey review committee), also collect, compile and review input demographic data before finalizing reports, (e.g., the final 2020 Decadal Survey Report). A summary of these data should be released alongside the final survey report. This information would enable the committee to understand potential "blind spots" and biases of the data collection phase and inform future data collections of any barriers that affect the omission of perspectives from various demographics.

The purpose of this white paper is to provide guidance to funding agencies, leaders in the discipline, and its constituent departments about strategies for (1) improving access to advanced education for people from populations that have long been underrepresented and (2) improving the climates of departments where students enroll. The twin goals of improving access to increase diversity and improving climate to enhance inclusiveness are mutually reinforcing, and they are both predicated on a fundamental problem of inequality in participation. This white paper has been endorsed by the Board of Trustees of the AAS.

Calculating the long term solution of ordinary differential equations, such as those of the $N$-body problem, is central to understanding a wide range of dynamics in astrophysics, from galaxy formation to planetary chaos. Because generally no analytic solution exists to these equations, researchers rely on numerical methods which are prone to various errors. In an effort to mitigate these errors, powerful symplectic integrators have been employed. But symplectic integrators can be severely limited because they are not compatible with adaptive stepping and thus they have difficulty accommodating changing time and length scales. A promising alternative is time-reversible integration, which can handle adaptive time stepping, but the errors due to time-reversible integration in astrophysics are less understood. The goal of this work is to study analytically and numerically the errors caused by time-reversible integration, with and without adaptive stepping. We derive the modified differential equations of these integrators to perform the error analysis. As an example, we consider the trapezoidal rule, a reversible non-symplectic integrator, and show it gives secular energy error increase for a pendulum problem and for a H\'{e}non---Heiles orbit. We conclude that using reversible integration does not guarantee good energy conservation and that, when possible, use of symplectic integrators is favored. We also show that time-symmetry and time-reversibility are properties that are distinct for an integrator.

We present a new symplectic integrator designed for collisional gravitational $N$-body problems which makes use of Kepler solvers. The integrator is also reversible and conserves 9 integrals of motion of the $N$-body problem to machine precision. The integrator is second order, but the order can easily be increased by the method of \citeauthor{yos90}. We use fixed time step in all tests studied in this paper to ensure preservation of symplecticity. We study small $N$ collisional problems and perform comparisons with typically used integrators. In particular, we find comparable or better performance when compared to the 4th order Hermite method and much better performance than adaptive time step symplectic integrators introduced previously. We find better performance compared to SAKURA, a non-symplectic, non-time-reversible integrator based on a different two-body decomposition of the $N$-body problem. The integrator is a promising tool in collisional gravitational dynamics.

We investigate the observational consequences of the quintessence field rolling to and oscillating near a minimum in its potential, "if" it happens close to the present epoch (z<0.2). We show that in a class of models, the oscillations lead to a rapid growth of the field fluctuations and the gravitational potential on subhorizon scales. The growth in the gravitational potential occurs on timescales << H^(1). This effect is present even when the quintessence parameters are chosen to reproduce an expansion history consistent with observations. For linearized fluctuations, we find that although the gravitational potential power spectrum is enhanced in a scale-dependent manner, the shape of the dark matter/galaxy power spectrum is not significantly affected. We find that the best constraints on such a transition in the quintessence field is provided via the integrated Sachs-Wolfe (ISW) effect in the CMB temperature power spectrum. Going beyond the linearized regime, the quintessence field can fragment into large, localized, long lived excitations (oscillons) with sizes comparable to galaxy clusters; this fragmentation could provide additional observational constraints. Two quoted "signatures" of modified gravity are a scale-dependent growth of the gravitational potential and a difference between the matter power spectrum inferred from measurements of lensing and galaxy clustering. Here, both effects are achieved by a minimally coupled scalar field in general relativity with a canonical kinetic term.

The metric of a perturbed Robertson-Walker spacetime is characterized by three functions: a scale-factor giving the expansion history and two potentials which generalize the single potential of Newtonian gravity. The Newtonian potential induces peculiar velocities and, from these, the growth of matter fluctuations. Massless particles respond equally to the Newtonian potential and to a curvature potential. The difference of the two potentials, called the gravitational slip, is predicted to be very small in general relativity but can be substantial in modified gravity theories. The two potentials can be measured, and gravity tested on cosmological scales, by combining weak gravitational lensing or the Integrated Sachs-Wolfe effect with galaxy peculiar velocities or clustering.

N-body simulations have revealed a wealth of information about dark matter halos however their results are largely empirical. Using analytic means, we attempt to shed light on simulation results by generalizing the self-similar secondary infall model to include tidal torque. In this first of two papers, we describe our halo formation model and compare our results to empirical mass profiles inspired by N-body simulations. Each halo is determined by four parameters. One parameter sets the mass scale and the other three define how particles within a mass shell are torqued throughout evolution. We choose torque parameters motivated by tidal torque theory and N-body simulations and analytically calculate the structure of the halo in different radial regimes. We find that angular momentum plays an important role in determining the density profile at small radii. For cosmological initial conditions, the density profile on small scales is set by the time rate of change of the angular momentum of particles as well as the halo mass. On intermediate scales, however, $\rho\propto r^{-2}$, while $\rho\propto r^{-3}$ close to the virial radius.

Using a generalized self-similar secondary infall model, which accounts for tidal torques acting on the halo, we analyze the velocity profiles of halos in order to gain intuition for N-body simulation results. We analytically calculate the asymptotic behavior of the internal radial and tangential kinetic energy profiles in different radial regimes. We then numerically compute the velocity anisotropy and pseudo-phase-space density profiles and compare them to recent N-body simulations. For cosmological initial conditions, we find both numerically and analytically that the anisotropy profile asymptotes at small radii to a constant set by model parameters. It rises on intermediate scales as the velocity dispersion becomes more radially dominated and then drops off at radii larger than the virial radius where the radial velocity dispersion vanishes in our model. The pseudo-phase-space density is universal on intermediate and large scales. However, its asymptotic slope on small scales depends on the halo mass and on how mass shells are torqued after turnaround. The results largely confirm N-body simulations but show some differences that are likely due to our assumption of a one-dimensional phase space manifold.

Indirect detection of high-energy particles from dark matter interactions is a promising avenue for learning more about dark matter, but is hampered by the frequent coincidence of high-energy astrophysical sources of such particles with putative high-density regions of dark matter. We calculate the boost factor and gamma-ray flux from dark matter associated with two shell-like caustics of luminous tidal debris recently discovered around the Andromeda galaxy, under the assumption that dark matter is its own supersymmetric antiparticle. These shell features could be a good candidate for indirect detection of dark matter via gamma rays because they are located far from the primary confusion sources at the galaxy's center, and because the shapes of the shells indicate that most of the mass has piled up near apocenter. Using a numerical estimator specifically calibrated to estimate densities in N-body representations with sharp features and a previously determined N-body model of the shells, we find that the largest boost factors do occur in the shells but are only a few percent. We also find that the gamma-ray flux is an order of magnitude too low to be detected with Fermi for likely dark matter parameters, and about 2 orders of magnitude less than the signal that would have come from the dwarf galaxy that produces the shells in the N-body model. We further show that the radial density profiles and relative radial spacing of the shells, in either dark or luminous matter, is relatively insensitive to the details of the potential of the host galaxy but depends in a predictable way on the velocity dispersion of the progenitor galaxy.

We develop a theoretical framework for describing the hierarchical structure of the phase space of cold dark matter haloes, due to gravitationally bound substructures. Because it includes the full hierarchy of the cold dark matter initial conditions and is hence complementary to the halo model, the stable clustering hypothesis is applied for the first time here to the small-scale phase space structure. As an application, we show that the particle dark matter annihilation signal could be up to two orders of magnitude larger than that of the smooth halo within the Galactic virial radius. The local boost is inversely proportional to the smooth halo density, and thus is O(1) within the solar radius, which could translate into interesting signatures for dark matter direct detection experiments: The temporal correlation of dark matter detection can change by a factor of 2 in the span of 10 years, while there will be significant correlations in the velocity space of dark matter particles. This can introduce O(1) uncertainty in the direction of local dark matter wind, which was believed to be a benchmark of directional dark matter searches or the annual modulation signal.

Most of the mass content of dark matter haloes is expected to be in the form of tidal debris. The density of debris is not constant, but rather can grow due to formation of caustics at the apocenters and pericenters of the orbit, or decay as a result of phase mixing. In the phase space, the debris assemble in a hierarchy which is truncated by the primordial temperature of dark matter. Understanding this phase structure can be of significant importance for the interpretation of many astrophysical observations and in particular dark matter detection experiments. With this purpose in mind, we develop a general theoretical framework to describe the hierarchical structure of the phase space of cold dark matter haloes. We do not make any assumption of spherical symmetry and/or smooth and continuous accretion. Instead, working with correlation functions in the action-angle space, we can fully account for the hierarchical structure (predicting a two-point correlation function ~ (\Delta J)^{-1.6} in the action space), as well as the primordial discreteness of the phase space. As an application, we estimate the boost to the dark matter annihilation signal due to the structure of the phase space within virial radius: the boost due to the hierarchical tidal debris is of order unity, whereas the primordial discreteness of the phase structure can boost the total annihilation signal by up to an order of magnitude. The latter is dominated by the regions beyond 20% of the virial radius, and is largest for the recently formed haloes with the least degree of phase mixing.

When thin accretion disks around black holes are perturbed, the main restoring force is gravity. If gas pressure, magnetic stresses, and radiation pressure are neglected, the disk remains thin as long as orbits do not intersect. Intersections would result in pressure forces which limit the growth of perturbations. We find that a discrete set of perturbations is possible for which orbits remain non-intersecting for arbitrarily long times. These modes define a discrete set of frequencies. We classify all long-lived perturbations for arbitrary potentials and show how their mode frequencies are related to pattern speeds computed from the azimuthal and epicyclic frequencies. We show that modes are concentrated near radii where the pattern speed has vanishing radial derivative. We explore these modes around Kerr black holes as a possible explanation for the high-frequency quasi-periodic oscillations of black hole binaries such as GRO J1655-40. The long-lived modes are shown to coincide with diskoseismic waves in the limit of small sound speed. While the waves have long lifetime, they have the wrong frequencies to explain the pairs of high-frequency quasi-periodic oscillations observed in black hole binaries.

The acceleration of the universe can be explained either through dark energy or through the modification of gravity on large scales. In this paper we investigate modified gravity models and compare their observable predictions with dark energy models. Modifications of general relativity are expected to be scale-independent on super-horizon scales and scale-dependent on sub-horizon scales. For scale-independent modifications, utilizing the conservation of the curvature scalar and a parameterized post-Newtonian formulation of cosmological perturbations, we derive results for large scale structure growth, weak gravitational lensing, and cosmic microwave background anisotropy. For scale-dependent modifications, inspired by recent $f(R)$ theories we introduce a parameterization for the gravitational coupling $G$ and the post-Newtonian parameter $\gamma$. These parameterizations provide a convenient formalism for testing general relativity. However, we find that if dark energy is generalized to include both entropy and shear stress perturbations, and the dynamics of dark energy is unknown a priori, then modified gravity cannot in general be distinguished from dark energy using cosmological linear perturbations.

This paper describes a fourth-order integration algorithm for the gravitational N-body problem based on discrete Lagrangian mechanics. When used with shared timesteps, the algorithm is momentum conserving and symplectic. We generalize the algorithm to handle individual time steps; this introduces fifth-order errors in angular momentum conservation and symplecticity. We show that using adaptive block power of two timesteps does not increase the error in symplecticity. In contrast to other high-order, symplectic, individual timestep, momentum-preserving algorithms, the algorithm takes only forward timesteps. We compare a code integrating an N-body system using the algorithm with a direct-summation force calculation to standard stellar cluster simulation codes. We find that our algorithm has about 1.5 orders of magnitude better symplecticity and momentum conservation errors than standard algorithms for equivalent numbers of force evaluations and equivalent energy conservation errors.

Weakly interacting massive particles are part of the lepton-photon plasma in the early universe until kinetic decoupling, after which time the particles behave like a collisionless gas with nonzero temperature. The Boltzmann equation for WIMP-lepton collisions is reduced to a Fokker-Planck equation for the evolution of the WIMP distribution including scalar density perturbations. This equation and the Einstein and fluid equations for the plasma are solved numerically including the acoustic oscillations of the plasma before and during kinetic decoupling, the frictional damping occurring during kinetic decoupling, and the free-streaming damping occurring afterwards and throughout the radiation-dominated era. An excellent approximation reduces the solution to quadratures for the cold dark matter density and velocity perturbations. The subsequent evolution is followed through electron pair annihilation and the radiation-matter transition; analytic solutions are provided for both large and small scales. For a 100 GeV WIMP with bino-type interactions, kinetic decoupling occurs at a temperature $T_d=23$ MeV. The transfer function in the matter-dominated era leads to an abundance of small cold dark matter halos; with a smooth window function the Press-Schechter mass distribution is $dn/d\ln M\propto M^{-1/3}$ for $M<10^{-4} (T_d/$10 MeV)$^{-3}$ M$_\odot$.

Treating the motion of a dust particle suspended in a liquid as a random walk, Einstein in 1905 derived an equation describing the diffusion of the particle's probability distribution in configuration space. Fokker and Planck extended this work to describe the velocity distribution of the particles. Their equation and its solutions have been applied to many problems in nature starting with the motion of Brownian particles in a liquid. Chandrasekhar derived the Fokker-Planck equation for stars and showed that long-range gravitational encounters provide a drag force, dynamical friction, which is important in the evolution of star clusters and the formation of galaxies. In certain circumstances, Fokker-Planck evolution also describes the evolution of dark (invisible) matter in the universe. In the early universe, the thermal decoupling of weakly interacting massive particles from the plasma of relativistic leptons and photons is governed by Fokker-Planck evolution. The resulting dissipation imprints a minimum length scale for cosmic density fluctuations. Still later, these density fluctuations produce stochastic gravitational forces on the dark matter as it begins to cluster under gravity. The latter example provides an exact derivation of the Fokker-Planck equation without the usual assumption of Markovian dynamics.

The strongest evidence for dark energy comes presently from geometric techniques such as the supernova distance-redshift relation. By combining the measured expansion history with the Friedmann equation one determines the energy density and its time evolution, hence the equation of state of dark energy. Because these methods rely on the Friedmann equation which has not been independently tested it is desirable to find alternative methods that work for both general relativity and other theories of gravity. Assuming that sufficiently large patches of a perturbed Robertson-Walker spacetime evolve like separate Robertson-Walker universes, that shear stress is unimportant on large scales and that energy and momentum are locally conserved, we derive several relations between long-wavelength metric and matter perturbations. These relations include generalizations of the initial-value constraints of general relativity. For a class of theories including general relativity we reduce the long-wavelength metric, density, and velocity potential perturbations to quadratures including curvature perturbations, entropy perturbations, and the effects of nonzero background curvature. When combined with the expansion history measured geometrically, the long-wavelength solution provide a test that may distinguish modified gravity from other explanations of dark energy.

We present a parallel implementation of the particle-particle/particle-mesh (P3M) algorithm for distributed memory clusters. The GRACOS (GRAvitational COSmology) code uses a hybrid method for both computation and domain decomposition. Long-range forces are computed using a Fourier transform gravity solver on a regular mesh; the mesh is distributed across parallel processes using a static one-dimensional slab domain decomposition. Short-range forces are computed by direct summation of close pairs; particles are distributed using a dynamic domain decomposition based on a space-filling Hilbert curve. A nearly-optimal method was devised to dynamically repartition the particle distribution so as to maintain load balance even for extremely inhomogeneous mass distributions. Tests using $800^3$ simulations on a 40-processor beowulf cluster showed good load balance and scalability up to 80 processes. We discuss the limits on scalability imposed by communication and extreme clustering and suggest how they may be removed by extending our algorithm to include adaptive mesh refinement.

A novel covariant formalism for the treatment of the transfer and Compton scattering of partially polarized light is presented. This was initially developed to aid in the computation of relativistic corrections to the polarization generated by the Sunyaev-Zeldovich effect (demonstrated in a companion paper), but it is of more general utility. In this approach, the polarization state of a light beam is described by a tensor constructed from the time average of quadratic products of the electric field components in a local observer frame. This leads naturally to a covariant description which is ideal for calculations involving the boosting of polarized light beams between Lorentz frames, and is more flexible than the traditional Stokes parameter approach in which a separate set of polarization basis vectors is required for each photon. The covariant kinetic equation for Compton scattering of partially polarized light by relativistic electrons is obtained in the tensor formalism by a heuristic semi-classical line of reasoning. The kinetic equation is derived first in the electron rest frame in the Thomson limit, and then is generalized to account for electron recoil and allow for scattering from an arbitrary distribution of electrons.

We apply the coherency tensor formalism to the calculation of the spectral distortions imprinted in the intensity and polarization of the cosmic microwave background radiation due to the kinematic and thermal Sunyaev-Zeldovich effects (SZE). We obtain the first relativistic corrections to the intensity produced by the kinematic and thermal SZE, and the first correction to the polarization magnitude due to electron thermal motion.

We present a kinetic theory for the evolution of the phase-space distribution of dark matter particles in galaxy halos in the presence of a cosmological spectrum of fluctuations. This theory introduces a new way to model the formation and evolution of halos, which traditionally have been investigated by analytic gravitational infall models or numerical N-body methods. Unlike the collisionless Boltzmann equation, our kinetic equation contains nonzero terms on the right-hand side arising from stochastic fluctuations in the gravitational potential due to substructures in the dark matter mass distribution. Using statistics for constrained Gaussian random fields in standard cosmological models, we show that our kinetic equation to second-order in perturbation theory is of the Fokker-Planck form, with one scattering term representing drift and the other representing diffusion in velocity-space. The drift is radial, and the drift and diffusion coefficients depend only on positions and not velocities; our relaxation process in the quasilinear regime is therefore different from the standard two-body relaxation. We provide explicit expressions relating these coefficients to the linear power spectrum of mass fluctuation and present results for the currently favored cold dark matter model with a nonzero cosmological constant. Solutions to this kinetic equation will provide a complete description of the cold dark matter spatial and velocity distributions for the average halo during the early phases of galaxy halo formation.

In at least two black hole binary systems, the Rossi X-Ray Timing Explorer has detected high frequency quasi-periodic oscillations (HFQPOs) with a 2:3 frequency commensurability. We propose a simple hot spot model to explain the positions, amplitudes, and widths of the HFQPO peaks. Using the exact geodesic equations for the Kerr metric, we calculate the trajectories of massive test particles, which are treated as isotropic, monochromatic emitters in their rest frames. By varying the hot spot parameters, we are able to explain the different features observed in ``Type A'' and ``Type B'' QPOs from XTE J1550-564. In the context of this model, the observed power spectra allow us to infer values for the black hole mass and angular momentum, and also constrain the parameters of the model.

Observations from the Rossi X-Ray Timing Explorer have shown the existence of high frequency quasi-periodic oscillations (HFQPOs) in the X-ray flux from accreting black hole binary systems. In at least two systems, these HFQPOs come in pairs with a 2:3 frequency commensurability. We propose a simple ``hot spot'' model to explain the position and amplitude of the HFQPO peaks. Using the exact geodesic equations for the Kerr metric, we calculate the trajectories of massive test particles, which are treated as isotropic, monochromatic emitters in their rest frames. Photons are traced from the accretion disk to a distant observer to produce time- and frequency-dependent images of the orbiting hot spot and background disk. The power spectrum of the X-ray light curve consists of multiple peaks at integral combinations of the black hole coordinate frequencies. In particular, if the radial frequency is one-third of the azimuthal frequency (as is the case near the inner-most stable circular orbit), beat frequencies appear in the power spectrum at two-thirds and four-thirds of the fundamental azimuthal orbital frequency, in agreement with observations. We also model the effects of shearing the hot spot in the disk, producing an arc of emission that also follows a geodesic orbit, as well as the effects of non-planar orbits that experience Lens-Thirring precession around the black hole axis. By varying the arc length, we are able to explain the different QPO features observed in ``Type A'' and ``Type B'' X-ray outbursts from XTE J1550-564. In the context of this model, the observed power spectra allow us to infer values for the black hole mass and angular momentum, and also constrain the parameters of the model, such as the hot spot size and luminosity.

We show that the linear dynamics of cosmological perturbations can be described by coupled wave equations, allowing their efficient numerical and, in certain limits, analytical integration directly in position space. The linear evolution of any perturbation can then be analyzed with the Green's function method. Prior to hydrogen recombination, assuming tight coupling between photons and baryons, neglecting neutrino perturbations, and taking isentropic (adiabatic) initial conditions, the obtained Green's functions for all metric, density, and velocity perturbations vanish beyond the acoustic horizon. A localized primordial cosmological perturbation expands as an acoustic wave of photon-baryon density perturbation with narrow spikes at its acoustic wavefronts. These spikes provide one of the main contributions to the cosmic microwave background radiation anisotropy on all experimentally accessible scales. The gravitational interaction between cold dark matter and baryons causes a dip in the observed temperature of the radiation at the center of the initial perturbation. We first model the radiation by a perfect fluid and then extend our analysis to account for finite photon mean free path. The resulting diffusive corrections smear the sharp features in the photon and baryon density Green's functions over the scale of Silk damping.

We suggest that the cosmic microwave background (CMB) temperature correlation function C(theta) as a function of angle provides a direct connection between experimental data and the fundamental cosmological quantities. The evolution of inhomogeneities in the prerecombination universe is studied using their Green's functions in position space. We find that a primordial adiabatic point perturbation propagates as a sharp-edged spherical acoustic wave. Density singularities at its wavefronts create a feature in the CMB correlation function distinguished by a dip at theta ~ 1.2 deg. Characteristics of the feature are sensitive to the values of cosmological parameters, in particular to the total and the baryon densities.

This paper describes the generation of initial conditions for numerical simulations in cosmology with multiple levels of resolution, or multiscale simulations. We present the theory of adaptive mesh refinement of Gaussian random fields followed by the implementation and testing of a computer code package performing this refinement called GRAFIC2. This package is available to the computational cosmology community at http://arcturus.mit.edu/grafic/ or by email from the author.

These lecture notes discuss several topics in the physics of cosmic structure formation starting from the evolution of small-amplitude fluctuations in the radiation-dominated era. The topics include relativistic cosmological perturbation theory with the scalar-vector-tensor decomposition, the evolution of adiabatic and isentropic initial fluctuations, microwave background anisotropy, spatial and angular power spectra, the cold dark matter linear transfer function, Press-Schechter theory, and a brief introduction to numerical simulation methods.

The power spectrum P(k)\propto k^n with n=-2 is close to the shape of the measured galaxy spectrum on small scales. Unfortunately this spectrum has proven rather difficult to simulate. Further, 2-dimensional simulations have suggested a breakdown of self-similar scaling for spectra with n<-1 due to divergent contributions from the coupling of long wave modes. This paper is the second (numerical) part of our investigation into the nonlinear gravitational clustering of scale-free spectra. Using high-resolution N-body simulations we find that the n=-2 power spectrum, as well as trajectories of the amplitude and phase of Fourier modes, display self-similar scaling. The evolution of the phase shift does show a different scaling at late times, but this was shown analytically to arise from the purely kinematical effect of bulk flows. Thus our analytical and N-body results verify that self-similarity in gravitational clustering holds for -3<n<1. The N-body spectrum is also compared with analytic fitting formulae, which are found to slightly underestimate the power in the nonlinear regime. The asymptotic shape of the spectrum at high-k is a power law with the same slope as predicted by the stable clustering hypothesis.

We present an efficient method to compute CMB anisotropies in non-flat universes. First we derive the Boltzmann equation for cosmic microwave background temperature and polarization fluctuations produced by scalar perturbations in a general Robertson-Walker universe. We then apply the integral method to solve this equation, writing temperature and polarization anisotropies as a time integral over a geometrical term and a source term. The geometrical terms can be written using ultra-spherical Bessel functions, which depend on curvature. These cannot be precomputed in advance as in flat space. Instead we solve directly their differential equation for selected values of the multipoles. The resulting computational time is comparable to the flat space case and improves over previous methods by 2-3 orders of magnitude. This allows one to compute highly accurate CMB temperature and polarization spectra, matter transfer functions and their CMB normalizations for any cosmological model, thereby avoiding the need to use various approximate fitting formulae that exist in the literature.

Any viable cosmological model must produce enough structure at early epochs to explain the amount of gas associated with high-redshift damped Ly$\alpha$ systems. We study the evolution of damped Ly$\alpha$ systems at redshifts $z\ge 2$ in cold dark matter (CDM) and cold+hot dark matter (CDM+HDM) models using both N-body and hydrodynamic simulations. Our approach incorporates the effects of gas dynamics, and we find that all earlier estimates which assumed that all the baryons in dark matter halos would contribute to damped Ly$\alpha$ absorption have overestimated the column density distribution $f(N)$ and the fraction of neutral dense gas $\Omega_g$ in damped Ly$\alpha$ systems. The differences are driven by ionization of hydrogen in the outskirts of galactic halos and by gaseous dissipation near the halo centers, and they tend to exacerbate the problem of late galaxy formation in CDM+HDM models. We only include systems up to the highest observed column density $N\sim 10^{21.8}$ cm$^{-2}$ in the estimation of $\Omega_g$ for a fair comparison with data. If the observed $f(N)$ and $\Omega_g$ inferred from a small number of confirmed and candidate absorbers are robust, the amount of gas in damped Ly$\alpha$ systems at high redshifts in the $\Omega_\nu=0.2$ CDM+HDM model falls well below the observations.

The EFAR project is a study of 736 candidate early-type galaxies in 84 clusters lying in two regions towards Hercules-Corona Borealis and Perseus-Cetus at distances $cz \approx 6000-15000$ km/s. In this paper we describe a new method of galaxy photometry adopted to derive the photometric parameters of the EFAR galaxies. The algorithm fits the circularized surface brightness profiles as the sum of two seeing-convolved components, an $R^{1/4}$ and an exponential law. This approach allows us to fit the large variety of luminosity profiles displayed by the EFAR galaxies homogeneously and to derive (for at least a subset of these) bulge and disk parameters. Multiple exposures of the same objects are optimally combined and an optional sky-fitting procedure has been developed to correct for sky subtraction errors. Extensive Monte Carlo simulations are analyzed to test the performance of the algorithm and estimate the size of random and {\it systematic} errors. Random errors are small, provided that the global signal-to-noise ratio of the fitted profiles is larger than $\approx 300$. Systematic errors can result from 1) errors in the sky subtraction, 2) the limited radial extent of the fitted profiles, 3) the lack of resolution due to seeing convolution and pixel sampling, 4) the use of circularized profiles for very flattened objects seen edge-on and 5) a poor match of the fitting functions to the object profiles. Large systematic errors are generated by the widely used simple $R^{1/4}$ law to fit luminosity profiles when a disk component, as small as 20% of the total light, is present.

We present a new method to add long wavelength power to an evolved $N$-body simulation, making use of the Zel'dovich (1970) approximation to change positions and velocities of particles. We describe the theoretical framework of our technique and apply it to a P$^3$M cosmological simulation performed on a cube of $100$ Mpc on a side, obtaining a new ``simulation'' of $800$ Mpc on a side. We study the effect of the power added by long waves by mean of several statistics of the density and velocity field, and suggest possible applications of our method to the study of the large-scale structure of the universe.

The EFAR project is a study of 736 candidate elliptical galaxies in 84 clusters lying in two regions towards Hercules-Corona Borealis and Perseus-Pisces-Cetus at distances cz = 6000-15000 km/s. In this paper (the first of a series) we present an introduction to the EFAR project and describe in detail the selection of the clusters and galaxies in our sample. Fundamental data for the galaxies and clusters are given, including accurate new positions for each galaxy and redshifts for each cluster. The galaxy selection functions are determined using diameters measured from Schmidt sky survey images for 2185 galaxies in the cluster fields. Future papers in this series will present the spectroscopic and photometric observations of this sample, investigate the properties of the fundamental plane for ellipticals, and determine the large-scale peculiar velocity fields in these two regions of the universe.

Building a self-gravitating hydrodynamic code as a combination of a hydrodynamic solver and a gravity solver is discussed. We show that straightforward combining those two solvers generally leads to a code that does not conserve energy locally, and instead a special Consistency Condition ought to be satisfied. A particular example of combining Softened Lagrangian Hydrodynamics (SLH) with a P3M gravity solver is used to demonstrate the effect of the Consistency Condition for a self-gravitating hydrodynamic code. The need to supplement the SLH method with the P3M gravity solver arose because the Moving Mesh Gravity solver, used in conjunction with the SLH method previously, was found to produce inaccurated results. We also show that most existing cosmological hydrodynamic codes implicitly satisfy the Consistency Condition.

We investigate three different local approximations for nonlinear gravitational instability in the framework of cosmological Lagrangian fluid dynamics of cold dust. They include the Zel'dovich approximation (ZA), the ``non-magnetic'' approximation of Bertschinger \& Jain (1994, NMA), and a new ``local tidal'' approximation (LTA). The LTA is exact for any perturbations whose gravitational and velocity equipotentials have the same constant shape with time, including spherical, cylindrical, and plane-parallel perturbations. We tested all three local approximations with the collapse of a homogeneous triaxial ellipsoid, for which an exact solution exists for an ellipsoid embedded in empty space and an excellent approximation is known in the cosmological context. We find that the LTA is significantly more accurate in general than the ZA and the NMA. Like the ZA, but unlike the NMA, the LTA generically leads to pancake collapse. For a randomly chosen mass element in an Einstein-de Sitter universe, assuming a Gaussian random field of initial density fluctuations, the LTA predicts that at least 78\% of initially underdense regions collapse owing to nonlinear effects of shear and tides.

We develop an algorithm for setting up initial Gaussian random density and velocity fields containing one or more peaks or dips, in an arbitrary cosmological scenario. The intention is to generate appropriate initial conditions for cosmological N-body simulations that focus on the evolution of the progenitors of the present-day galaxies and clusters. The procedure is an application of the direct and accurate prescription of Hoffman & Ribak (1991) for generating constrained random fields. For each peak a total of 21 physical characteristics can be specified, including its scale, position, density Hessian, velocity, and velocity gradient. The velocity (or, equivalently, gravity) field constrants are based on a generalization of the formalism developed by Bardeen et al. (1986). The resulting density field is sculpted such that it induces the desired amount of net gravitational and tidal forces. We provide a detailed mathematical presentation of the formalism. Afterwards we provide analytical estimates of the likelihood of the imposed constraints. Amongst others, it is shown that the tidal field has a strong tendency to align itself along the principal axes of the mass tensor. The method is illustrated by means of some concrete examples. In addition to the illustration of constraint-field correlation functions and how they add up to the mean fields, followed by illustrations of the variance characteristics of field realizations, we concentrate in particular on the consequences of imposing gravitational field constraints (or, equivalent in the linear regime for growing mode fluctuations, peculiar velocity field constraints).

COSMICS is a package of fortran programs useful for computing transfer functions and microwave background anisotropy for cosmological models, and for generating gaussian random initial conditions for nonlinear structure formation simulations of such models. Four programs are provided: {\bf linger\_con} and {\bf linger\_syn} integrate the linearized equations of general relativity, matter, and radiation in conformal Newtonian and synchronous gauge, respectively; {\bf deltat} integrates the photon transfer functions computed by the linger codes to produce photon anisotropy power spectra; and {\bf grafic} tabulates normalized matter power spectra and produces constrained or unconstrained samples of the matter density field. Version 1.0 of COSMICS is available at http://arcturus.mit.edu/cosmics/ . The current release gives fortran-77 programs that run on workstations and vectorized supercomputers. Unix makefiles are included that make it simple to build and test the package. A future release will include portable parallel versions of the linger codes using standard message-passing libraries.

This paper presents a systematic treatment of the linear theory of scalar gravitational perturbations in the synchronous gauge and the conformal Newtonian (or longitudinal) gauge. It differs from others in the literature in that we give, in both gauges, a complete discussion of all particle species that are relevant to any flat cold dark matter (CDM), hot dark matter (HDM), or CDM+HDM models (including a possible cosmological constant). The particles considered include CDM, baryons, photons, massless neutrinos, and massive neutrinos (an HDM candidate), where the CDM and baryons are treated as fluids while a detailed phase-space description is given to the photons and neutrinos. Particular care is applied to the massive neutrino component, which has been either ignored or approximated crudely in previous works. Isentropic initial conditions on super-horizon scales are derived. The coupled, linearized Boltzmann, Einstein and fluid equations that govern the evolution of the metric and density perturbations are then solved numerically in both gauges for the standard CDM model and two CDM+HDM models with neutrino mass densities $\onu=0.2$ and 0.3, assuming a scale-invariant, adiabatic spectrum of primordial fluctuations. We also give the full details of the cosmic microwave background anisotropy, and present the first accurate calculations of the angular power spectra in the two CDM+HDM models including photon polarization, higher neutrino multipole moments, and helium recombination. The numerical programs for both gauges are available at http://arcturus.mit.edu/cosmics/ .

We have developed a code which links scalar-mode fluctuations in the early universe with those observable at the present time by integrating the coupled, linearized, Einstein, Boltzmann, and fluid equations in a perturbed flat Robertson-Walker spacetime. The results are useful both for calculations of the cosmic microwave background anisotropy and the linear power spectrum of matter fluctuations. This paper introduces the serial and parallel codes, presents timing results, and gives sample output in graphical and animated form. This is a latex version containing most of an HTML-format technical paper submitted for Supercomputing '95. The preferred method of viewing is to point your WWW client to http://arcturus.mit.edu/SC95/

This is a set of lecture notes basd on the lectures on cosmological dynamics given by E. Bertschinger at Les Houches in August 1993. The contents include elementary mechanics in cosmology, Eulerian and Lagrangian fluid dynamics, hot dark matter, and relativistic cosmological perturbation theory. To typeset the notes one must first obtain the style files and figures in file figmac.uu sent separately. Place all the files in one directory, run latex three times, and run dvips or equivalent to produce a postscript file.

We present a computational algorithm for computing short range forces between particles. The algorithm has two distinguishing features. First, it is optimized for multi-processor computers, and will use as many processors as are available. Second, it is optimized for inhomogeneous, dynamic particle distributions; for any distribution the computational load is distributed evenly to all processors, and the communication time is less than 15\% of the total run time. In this talk we present our new algorithm. We developed the program for a grand-challenge problem in cosmology, simulation of the formation of large-scale structure in the universe. This simulation, run on the Thinking Machines Corporation CM-5, uses the particle-particle/particle-mesh (PPPM) \cite{hock_east} algorithm. The particle-particle phase is computed using the algorithm we describe in this paper. We discuss this and other applications.

The gravitational evolution of scale free initial spectra $P(k)\propto k^n$ in an Einstein-de Sitter universe is widely believed to be self-similar for $-3<n<4$. However, for $-3<n<-1$ the existence of self-similar scaling has not been adequately demonstrated. Here we investigate the possible breaking of self-similar scaling due to the nonlinear contributions of long wave modes. For $n<-1$ the nonlinear terms in the Fourier space fluid equations contain terms that diverge due to contributions from wavenumber $k\to 0$ (the long wave limit). To assess the possible dynamical effects of this divergence the limit of long wave contributions is investigated in detail using two different analytical approaches. Perturbative contributions to the power spectrum are examined. It is shown that for $n<-1$ there are divergent contributions at all orders. However, at every order the leading order divergent terms cancel out exactly. This does not rule out the existence of a weaker but nevertheless divergent net contribution. The second approach consists of a non-perturbative approximation, developed to study the nonlinear effects of long wave mode coupling. A solution for the phase shift of the Fourier space density is obtained which is divergent for $n<-1$. A kinematical interpretation of the divergence of the phase shift, related to the translational motion induced by the large-scale bulk velocity, is given. Our analysis indicates that the amplitude of the density is {\it not} affected by the divergent terms. Thus both analytical approaches lead to the conclusion that the self-similar scaling of physically relevant measures of the growth of density perturbations is preserved.

We examine high-resolution gravitational N-body simulations of the $\Omega=1$ cold dark matter (CDM) model in order to determine whether there is any normalization of the initial density fluctuation spectrum that yields acceptable results for galaxy clustering and velocities. Dense dark matter halos in the evolved mass distribution are identified with luminous galaxies; the most massive halos are also considered as sites for galaxy groups, with a range of possibilities explored for the group mass to light ratios. We verify the earlier conclusions of White et al. (1987) for the low amplitude (high bias) CDM model --- the galaxy correlation function is marginally acceptable but that there are too many galaxies. We also show that the peak biasing method does not accurately reproduce the results obtained using dense halos identified in the simulations themselves. The COBE anisotropy implies a higher normalization, resulting in problems with excessive pairwise galaxy velocity dispersion unless a strong velocity bias is present. Although we confirm the strong velocity bias of halos reported by Couchman \& Carlberg (1992), we show that the galaxy motions are still too large on small scales. We find no amplitude for which the CDM model can reconcile simultaneously the galaxy correlation function, the low pairwise velocity dispersion, and the richness distribution of groups and clusters. With the normalization implied by COBE, the CDM spectrum has too much power on small scales if $\Omega=1$.

We use numerical simulations of critically-closed cold dark matter (CDM) models to study the effects of numerical resolution on observable quantities. We study simulations with up to $256^3$ particles using the particle-mesh (PM) method and with up to $144^3$ particles using the adaptive particle-particle --particle-mesh (P$^3$M) method. Comparisons of galaxy halo distributions are made among the various simulations. We also compare distributions with observations and we explore methods for identifying halos, including a new algorithm that finds all particles within closed contours of the smoothed density field surrounding a peak. The simulated halos show more substructure than predicted by the Press-Schechter theory. We are able to rule out all $\Omega=1$ CDM models for linear amplitude $\sigma_8\gsim 0.5$ because the simulations produce too many massive halos compared with the observations. The simulations also produce too many low mass halos. The distribution of halos characterized by their circular velocities for the P$^3$M simulations is in reasonable agreement with the observations for $150\kms\lsim V_{\rm circ} \lsim 350\kms$.}}

The abundance of galactic systems at high redshifts can impose a strong constraint on the cold+hot dark matter (CDM+HDM) models. The hot component reduces the excessive small-scale power in the COBE-normalized CDM model but also delays the epoch of galaxy formation. We present results from the first numerical simulations that have enough dynamic range to address accurately the issue of high-redshift halo abundances in CDM+HDM models. Equivalent high-resolution particle-particle/particle-mesh $N$-body simulations are performed for spatially flat models with $\Omega_\nu =0.3$ and 0.2 (with $H_0=50$ km s$^{-1}$ Mpc$^{-1}$ and $\Omega_b=0.05$). We study the constraints placed on the models by the high-redshift quasar space density and by the mass fraction in neutral dense gas associated with damped Ly$\alpha$ systems. We find that even with optimistic assumptions, the much-studied $\Omega_\nu=0.3$ model does not produce enough massive halos to account for the observed abundance of quasars at $z>4$. The model passes this test if $\Omega_\nu$ is decreased to 0.2. Both models do not produce enough high column-density halos to account for the amount of gas in damped Ly$\alpha$ systems at $z\go 3$: the $\Omega_\nu=0.3$ model falls short by a factor $\sim$80; the $\Omega_\nu=0.2$ model by a factor $\sim$3. We conclude that only CDM+HDM models with $\Omega_\nu\lo 0.2$ can match observations at high redshift, implying an upper bound of 4.7 eV on the most massive light neutrino (presumably the $\tau$).

We derive the evolution equations for the electric and magnetic parts of the Weyl tensor for cold dust from both general relativity and Newtonian gravity. In a locally inertial frame at rest in the fluid frame, the Newtonian equations agree with those of general relativity. We give explicit expressions for the electric and magnetic parts of the Weyl tensor in the Newtonian limit. In general, the magnetic part does not vanish, implying that the Lagrangian evolution of the fluid is not purely local.

We solve the nonlinear evolution of pressureless, irrotational density fluctuations in a perturbed Robertson-Walker spacetime using a new Lagrangian method based on the velocity gradient and gravity gradient tensors. Borrowing results from general relativity, we obtain a set of Newtonian ordinary differential equations for these quantities following a given mass element. Using these Lagrangian fluid equations we prove the following results: (1) The spherical tophat perturbation, having zero shear, is the slowest configuration to collapse for a given initial density and growth rate. (2) Initial density maxima are not generally the sites where collapse first occurs. (3) Initially underdense regions may undergo collapse if the shear is not too small. If the magnetic part of the Weyl tensor vanishes, the nonlinear evolution is described purely locally by our equations; this condition holds for spherical, cylindrical, and planar perturbations and may be a good approximation in other circumstances. Assuming the vanishing of the magnetic part of the Weyl tensor, we compute the exact nonlinear gravitational evolution of cold matter. We find that 56\% of initially underdense regions collapse in an Einstein-de Sitter universe for a homogeneous and isotropic random field. We also show that, given this assumption, the final stage of collapse is generically two-dimensional, leading to strongly prolate filaments rather than Zel'dovich pancakes. While this result may explain the prevalence of filamentary collapses in N-body simulations, it is not true in general, suggesting that the magnetic part of the Weyl tensor need not vanish in the Newtonian limit.

We present a systematic treatment of the linear theory of scalar gravitational perturbations in the synchronous gauge and the conformal Newtonian (or longitudinal) gauge. We first derive the transformation law relating the two gauges. We then write down in parallel in both gauges the coupled, linearized Boltzmann, Einstein and fluid equations that govern the evolution of the metric perturbations and the density fluctuations of the particle species. The particle species considered include cold dark matter (CDM), baryons, photons, massless neutrinos, and massive neutrinos (a hot dark matter or HDM candidate), where the CDM and baryon components are treated as fluids while a detailed phase-space description is given to the photons and neutrinos. The linear evolution equations presented are applicable to any $\Omega=1$ model with CDM or a mixture of CDM and HDM. Isentropic initial conditions on super-horizon scales are derived. The equations are solved numerically in both gauges for a CDM+HDM model with $\Omega_{\rm cold}=0.65,$ $\Omega_{\rm hot}=0.3$, and $\Omega_{\rm baryon}=0.05$. We discuss the evolution of the metric and the density perturbations and compare their different behaviors outside the horizon in the two gauges. In a companion paper we integrate the geodesic equations for the neutrino particles in the perturbed conformal Newtonian background metric computed here. The purpose is to obtain an accurate sampling of the neutrino phase space for the HDM initial conditions in $N$-body simulations of the CDM+HDM models.

This article reviews the prevailing paradigm for how galaxies and larger structures formed in the universe: gravitational instability. Basic observational facts are summarized to motivate the standard cosmological framework underlying most detailed investigations of structure formation. The observed universe approaches spatial uniformity on scales larger than about $10^{26}$ cm. On these scales gravitational dynamics is almost linear and therefore relatively easy to relate to observations of large-scale structure. On smaller scales cosmic structure is complicated not only by nonlinear gravitational clustering but also by nonlinear nongravitational gas dynamical processes. The complexity of these phenomena makes galaxy formation one of the grand challenge problems of the physical sciences. No fully satisfactory theory can presently account in detail for the observed cosmic structure. However, as this article summarizes, significant progress has been made during the last few years.

The Eulerian cosmological fluid equations are used to study the nonlinear mode coupling of density fluctuations. We evaluate the second-order power spectrum including all four-point contributions. In the weakly nonlinear regime we find that the dominant nonlinear contribution for realistic cosmological spectra is made by the coupling of long-wave modes and is well estimated by second order perturbation theory. For a linear spectrum like that of the cold dark matter model, second order effects cause a significant enhancement of the high $k$ part of the spectrum and a slight suppression at low $k$ near the peak of the spectrum. Our perturbative results agree well in the quasilinear regime with the nonlinear spectrum from high-resolution N-body simulations. We find that due to the long-wave mode coupling, characteristic nonlinear masses grow less slowly in time (i.e., are larger at higher redshifts) than would be estimated using the linear power spectrum. For the cold dark matter model at $(1+z)=(20,10,5,2)$ the nonlinear mass is about $(180,8,2.5,1.6)$ times (respectively) larger than a linear extrapolation would indicate, if the condition rms $\delta\rho/\rho =1$ is used to define the nonlinear scale. At high redshift the Press-Schechter mass distribution significantly underestimates the abundance of high-mass objects for the cold dark matter model. Although the quantitative results depend on the definition of the nonlinear scale, these basic consequences hold for any initial spectrum whose post-recombination spectral index $n$ decreases sufficiently rapidly with increasing $k$, a feature which arises quite generally during the transition from a radiation- to matter-dominated universe.

We use the POTENT reconstruction of the mass density field in the nearby universe to estimate the amplitude of the density fluctuation power spectrum for various cosmological models. We find sigma_8\Omega_m^{0.6}= 1.3^{+0.4}_{-0.3}, almost independently of the power spectrum. This value agrees well with the COBE normalization for the standard CDM model, while some alternative models predict an excessive amplitude compared with COBE. Flat low \Omega_m models and tilted models with spectral index n<0.8 are particularly discordant.

We study the quasilinear evolution of the one-point probability density functions (PDFs) of the smoothed density and velocity fields in a cosmological gravitating system beginning with Gaussian initial fluctuations. Our analytic results are based on the Zel'dovich approximation and laminar flow. A numerical analysis extends the results into the multistreaming regime using the smoothed fields of a CDM N-body simulation. We find that the PDF of velocity, both Lagrangian and Eulerian, remains Gaussian under the laminar Zel'dovich approximation, and it is almost indistinguishable from Gaussian in the simulations. The PDF of mass density deviates from a normal distribution early in the quasilinear regime and it develops a shape remarkably similar to a lognormal distribution with one parameter, the \rms density fluctuation $\sigma$. Applying these results to currently available data we find that the PDFs of the velocity and density fields, as recovered by the \pot\ procedure from observed velocities assuming $\Omega=1$, or as deduced from a redshift survey of \iras\ galaxies assuming that galaxies trace mass, are consistent with Gaussian initial fluctuations.

We describe an efficient Particle-Mesh algorithm for the Connection Machine CM-5. Our particular method parallelizes well and the computation time per time step decreases as the particles become more clustered. We achieve floating-point computation rates of 4--5 MFlops/sec/processing node and total operations (the sum of floating-point and integer arithmetic plus communications) of 5--10 MOps/sec/processing node. The rates scale almost linearly from 32 to 256 processors. Although some of what we discuss is specific to the CM-5, many aspects (e.g., the computation of the force on a mesh) are generic to all implementations, and other aspects (e.g., the algorithm for assignment of the density to the mesh) are useful on any parallel computer.

The cosmological fluid equations are used to study the nonlinear mode coupling of density fluctuations. We find that for realistic cosmological spectra there is a significant contribution to the nonlinear evolution on scales of interest to large-scale structure from the long-wave part of the initial spectrum. A consequence of this mode coupling is that at high redshift, $z$, the nonlinear scale [defined by $\sigma(z)=1$] can be significantly larger than a linear extrapolation would indicate. For the standard CDM spectrum with a $\sigma_8=1$ normalization the mass corresponding to the nonlinear scale at $z=20, 10, 5$ is about $100, 10, 3$ times (respectively) larger than the linearly extrapolated value. We also investigate the possibility of divergent contributions to the density field from long-wave modes if the spectral index of the power spectrum $\ n<-1$. Using an approximate non-perturbative approach we find that for $n>-3\,$ the divergent contribution appears only in the phase. This can be related to the large-scale bulk velocity, and clarifies previous results from N-body simulations.

This paper reviews the essential physics of gravitational instability in a Robertson-Walker background spacetime. Three approaches are presented in a pedagogical manner, based on (1) the Eulerian fluid equations, (2) the Lagrangian description of trajectories, and (3) the Lagrangian fluid equations. Linear and nonlinear limits are discussed for each case. Shear and tides are shown to play a key role in nonlinear gravitational instability. The Lagrangian fluid approach is used to show that several widely held beliefs about gravitational instability are false. The following collapse theorem is proven: for a given initial density fluctuation and growth rate, the spherical tophat perturbation collapses more slowly than any other configuration. We also show that density maxima are not the first points to collapse and that underdense regions may collapse if their initial shear is sufficiently high. The Lagrangian fluid approach leads to an almost closed set of local evolution equations for individual mass elements. The magnetic part of the Weyl tensor, which may be present even in the nonrelativistic (Newtonian) limit, may prevent a purely local description. However, neglecting the magnetic part of the Weyl tensor, we obtain predictions for high-redshift collapse that are in good agreement with a high-resolution cold dark matter N-body simulation.

We have used maximum-likelihood estimation to determine the quadrupole amplitude $Q_{\rm rms-PS}$ and the spectral index $n$ of the density fluctuation power spectrum at recombination from the \cobe\ DMR data. We find a strong correlation between the two parameters of the form $Q_{\rm rms-PS}=(15.7\pm 2.6)\exp[0.46(1-n)]\mk$ for fixed $n$. Our result is slightly smaller than and has a smaller statistical uncertainty than the 1992 estimate of Smoot et al.

This paper presents a general-relativistic N-body technique for evolving the phase space distribution of massive neutrinos in linear perturbation theory. The method provides a much more accurate sampling of the neutrino phase space for the HDM initial conditions of N-body simulations in a cold+hot dark matter universe than previous work. Instead of directly sampling the phase space at the end of the linear era, we first compute the evolution of the metric perturbations by numerically integrating the coupled, linearized Einstein, Boltzmann, and fluid equations for all particle species. We then sample the phase space shortly after neutrino decoupling at redshift z=10^9 when the distribution is Fermi-Dirac. To follow the trajectory of each neutrino, we subsequently integrate the geodesic equations for each neutrino in the perturbed background spacetime from z=10^9 to z=13.55, using the linearized metric found in the previous calculation to eliminate discreteness noise. The positions and momenta resulting from this integration represent a fair sample of the full neutrino phase space and can be used as HDM initial conditions for N-body simulations of nonlinear structure evolution in this model. A total of 21 million neutrino particles are used in a 100 Mpc box, with Omega_cdm=0.65, Omega_hdm=0.30, Omega_baryon=0.05, and Hubble constant H_0=50. We find that correlations develop in the neutrino densities and momenta which are absent when only the zeroth-order Fermi-Dirac distribution is considered.