phy-4660/feedback/papers/0504097/main.tex

%auto-ignore


\bibliographystyle{naturemag}

\onecolumn

\section*{\ \newline \Large Simulating the joint evolution of quasars, galaxies\vspace*{0.1cm}\newline
    and their large-scale distribution}

%\baselineskip16pt

\noindent{\sffamily Volker~Springel$^{1}$, %
Simon~D.~M.~White$^{1}$, %
Adrian~Jenkins$^{2}$, %
Carlos~S.~Frenk$^{2}$, \newline%
Naoki~Yoshida$^{3}$, %
Liang~Gao$^{1}$, %
Julio~Navarro$^{4}$, % 
Robert~Thacker$^{5}$, %
Darren~Croton$^{1}$, \newline%
John~Helly$^{2}$, %
John~A.~Peacock$^{6}$, %
Shaun~Cole$^{2}$, %
Peter~Thomas$^{7}$, % 
Hugh~Couchman$^{5}$, \newline% 
August~Evrard$^{8}$, % 
J\"org~Colberg$^{9}$ \& %
Frazer~Pearce$^{10}$}
\ \\

\noindent%
{\footnotesize\it%
$^{1}${Max-Planck-Institute for Astrophysics, Karl-Schwarzschild-Str. 1, 85740 Garching, Germany}\\
$^{2}${Inst. for Computational Cosmology, Dep. of Physics, Univ. of
  Durham, South Road, Durham  DH1 3LE, UK}\\
$^{3}${Department of Physics, Nagoya University, Chikusa-ku, Nagoya
  464-8602, Japan}\\
$^{4}${Dep. of Physics \& Astron., University of
    Victoria, Victoria, BC, V8P 5C2, Canada}\\
$^{5}${Dep. of Physics \& Astron., McMaster Univ., 1280 Main
  St. West, Hamilton, Ontario, L8S 4M1, Canada}\\
$^{6}${Institute of Astronomy, University of Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK}\\
$^{7}${Dep. of Physics \& Astron., University of
    Sussex, Falmer, Brighton BN1 9QH, UK}\\
$^{8}${Dep. of Physics \& Astron., Univ. of Michigan,
    Ann Arbor, MI 48109-1120, USA}\\
$^{9}${Dep. of Physics \& Astron., Univ. of Pittsburgh,
    3941 O'Hara Street, Pittsburgh PA 15260, USA}\\
$^{10}${Physics and Astronomy Department, Univ. of Nottingham,
    Nottingham NG7 2RD, UK}\\
}



\baselineskip26pt % On-and-half-space the manuscript.
\setlength{\parskip}{12pt}
\setlength{\parindent}{22pt}%




\noindent{\bf The cold dark matter model has become the leading
  theoretical paradigm for the formation of structure in the Universe.
  Together with the theory of cosmic inflation, this model makes a
  clear prediction for the initial conditions for structure formation
  and predicts that structures grow hierarchically through
  gravitational instability.  Testing this model requires that the
  precise measurements delivered by galaxy surveys can be compared to
  robust and equally precise theoretical calculations. Here we present
  a novel framework for the quantitative physical interpretation of
  such surveys. This combines the largest simulation of the growth of
  dark matter structure ever carried out with new techniques for
  following the formation and evolution of the visible components.  We
  show that baryon-induced features in the initial conditions of the
  Universe are reflected in distorted form in the low-redshift galaxy
  distribution, an effect that can be used to constrain the nature of
  dark energy with next generation surveys.}




Recent large surveys such as the 2 degree Field Galaxy Redshift Survey
(2dFGRS) and the Sloan Digital Sky Survey (SDSS) have characterised
much more accurately than ever before not only the spatial clustering,
but also the physical properties of low-redshift galaxies.  Major
ongoing campaigns exploit the new generation of 8m-class telescopes
and the Hubble Space Telescope to acquire data of comparable quality
at high redshift. Other surveys target the weak image shear caused by
gravitational lensing to extract precise measurements of the
distribution of dark matter around galaxies and galaxy clusters. The
principal goals of all these surveys are to shed light on how galaxies
form, to test the current paradigm for the growth of cosmic structure,
and to search for signatures which may clarify the nature of dark
matter and dark energy.  These goals can be achieved only if the
accurate measurements delivered by the surveys can be compared to
robust and equally precise theoretical predictions.  Two problems have
so far precluded such predictions: (i) accurate estimates of
clustering require simulations of extreme dynamic range, encompassing
volumes large enough to contain representative populations of rare
objects (like rich galaxy clusters or quasars), yet resolving the
formation of individual low luminosity galaxies; (ii) critical aspects
of galaxy formation physics are uncertain and beyond the reach of
direct simulation (for example, the structure of the interstellar
medium, its consequences for star formation and for the generation of
galactic winds, the ejection and mixing of heavy elements, AGN feeding
and feedback effects \ldots) -- these must be treated by
phenomenological models whose form and parameters are adjusted by
trial and error as part of the overall data-modelling process. We have
developed a framework which combines very large computer simulations
of structure formation with post-hoc modelling of galaxy formation
physics to offer a practical solution to these two entwined problems.
  
During the past two decades, the cold dark matter (CDM) model,
augmented with a dark energy field (which may take the form of a
cosmological constant `$\Lambda$'), has developed into the standard
theoretical paradigm for galaxy formation.  It assumes that structure
grew from weak density fluctuations present in the otherwise
homogeneous and rapidly expanding early universe.  These fluctuations
are amplified by gravity, eventually turning into the rich structure
that we see around us today.  Confidence in the validity of this model
has been boosted by recent observations.  Measurements of the cosmic
microwave background (CMB) by the WMAP satellite\cite{Bennett2003}
were combined with the 2dFGRS to confirm the central tenets of the
model and to allow an accurate determination of the geometry and
matter content of the Universe about $380\,000$ years after the Big
Bang\cite{Spergel2003}.  The data suggest that the early density
fluctuations were a Gaussian random field, as predicted by
inflationary theory, and that the current energy density is dominated
by some form of dark energy. This analysis is supported by the
apparent acceleration of the current cosmic expansion inferred from
studies of distant supernovae\cite{Riess1998,Perlmutter1999}, as well
as by the low matter density derived from the baryon fraction of
clusters\cite{White1993}.

While the initial, linear growth of density perturbations can be
calculated analytically, the collapse of fluctuations and the
subsequent hierarchical build-up of structure is a highly nonlinear
process which is only accessible through direct numerical
simulation\cite{Davis1985}.  The dominant mass component, the cold
dark matter, is assumed to be made of elementary particles that
currently interact only gravitationally, so the collisionless dark
matter fluid can be represented by a set of discrete point
particles. This representation as an N-body system is a coarse
approximation whose fidelity improves as the number of particles in
the simulation increases.  The high-resolution simulation described
here -- dubbed the {\it Millennium Simulation} because of its size --
was carried out by the Virgo Consortium, a collaboration of British,
German, Canadian, and US astrophysicists. It follows $N= 2160^3\simeq
1.0078\times 10^{10}$ particles from redshift $z=127$ to the present
in a cubic region $500\,h^{-1}{\rm Mpc}$ on a side, where $1+z$ is the
expansion factor of the Universe relative to the present and $h$ is
Hubble's constant in units of $100\,{\rm km\,s^{-1}Mpc^{-1}}$. With
ten times as many particles as the previous largest computations of
this kind\cite{Colberg2000,Evrard2002,Wambsganss2004} (see
Supplementary Information), it offers substantially improved spatial
and time resolution within a large cosmological volume. Combining this
simulation with new techniques for following the formation and
evolution of galaxies, we predict the positions, velocities and
intrinsic properties of all galaxies brighter than the Small
Magellanic Cloud throughout volumes comparable to the largest current
surveys.  Crucially, this also allows us to establish evolutionary
links between objects observed at different epochs.  For example, we
demonstrate that galaxies with supermassive central black holes can
plausibly form early enough in the standard cold dark matter cosmology
to host the first known quasars, and that these end up at the centres
of rich galaxy clusters today.


\begin{figure*}
\noindent\hspace*{-0.5cm}%
\resizebox{17.0cm}{!}{\includegraphics{fig1.eps}} %
\caption{The dark matter density field on various scales.  Each
  individual image shows the projected dark matter density field in a
  slab of thickness $15\,h^{-1}{\rm Mpc}$ (sliced from the periodic
  simulation volume at an angle chosen to avoid replicating structures
  in the lower two images), colour-coded by density and local dark
  matter velocity dispersion.  The zoom sequence displays consecutive
  enlargements by factors of four, centred on one of the many galaxy
  cluster halos present in the simulation.}
\label{FigDMDist}
\end{figure*}


\subsubsection*{Dark matter halos and galaxies}
The mass distribution in a $\Lambda$CDM universe has a complex
topology, often described as a ``cosmic web'' \cite{Bond1996}. This is
visible in full splendour in Fig.~\ref{FigDMDist} (see also the
corresponding Supplementary Video). The zoomed out panel at the bottom
of the figure reveals a tight network of cold dark matter clusters and
filaments of characteristic size $\sim 100\,h^{-1} {\rm Mpc}$. On
larger scales, there is little discernible structure and the
distribution appears homogeneous and isotropic.  Subsequent images
zoom in by factors of four onto the region surrounding one of the many
rich galaxy clusters.  The final image reveals several hundred dark
matter substructures, resolved as independent, gravitationally bound
objects orbiting within the cluster halo. These substructures are the
remnants of dark matter halos that fell into the cluster at earlier
times.


\begin{figure}
\hspace*{-1.0cm}%
\resizebox{16.0cm}{!}{\includegraphics{fig2.eps}}
\caption{ \baselineskip20pt Differential halo number density as a
function of mass and epoch.  The function $n(M,z)$ gives the comoving
number density of halos less massive than $M$. We plot it as the halo
multiplicity function $M^2\rho^{-1}\,{\rm d}n/{\rm d}M$, where $\rho$
is the mean density of the universe. Groups of particles were found
using a friends-of-friends algorithm\cite{Davis1985} with linking
length equal to 0.2 of the mean particle separation. The fraction of
mass bound to halos of more than 20 particles (vertical dotted line)
grows from $6.42\times 10^{-4}$ at $z=10.07$ to 0.496 at $z=0$. Solid
lines are predictions from an analytic fitting function proposed in
previous work\cite{Jenkins2001}, while the dashed lines give the
Press-Schechter model\cite{Press1974} at $z=10.07$ and $z=0$.
\label{FigMassFunc} }
\end{figure}


The space density of dark matter halos at various epochs in the
simulation is shown in Fig.~\ref{FigMassFunc}. At the present day,
there are about 18 million halos above a detection threshold of 20
particles; 49.6\% of all particles are included in these halos.  These
statistics provide the most precise determination to date of the mass
function of cold dark matter halos\cite{Jenkins2001,Reed2003}. In the
range that is well sampled in our simulation ($z \le 12$, $M\ge
1.7\times 10^{10}\,h^{-1}{\rm M}_\odot$), our results are remarkably
well described by the analytic formula proposed by Jenkins et
al.\cite{Jenkins2001} from fits to previous simulations. Theoretical
models based on an ellipsoidal excursion set
formulation\cite{Sheth2002} give a less accurate, but still reasonable
match. However, the commonly used Press-Schechter
formula\cite{Press1974} underpredicts the high-mass end of the mass
function by up to an order of magnitude. Previous studies of the
abundance of rare objects, such as luminous quasars or clusters, based
on this formula may contain large errors\cite{Efstathiou1988}. We
return below to the important question of the abundance of quasars at
early times.

To track the formation of galaxies and quasars in the simulation, we
implement a semi-analytic model to follow gas, star and supermassive
black hole processes within the merger history trees of dark matter
halos and their substructures (see Supplementary Information).  The
trees contain a total of about 800 million nodes, each corresponding
to a dark matter subhalo and its associated galaxies. This methodology
allows us to test, during postprocessing, many different
phenomenological treatments of gas cooling, star formation, AGN
growth, feedback, chemical enrichment, etc. Here, we use an update of
models described in\cite{Springel2001b,Kauffmann2000} which are
similar in spirit to previous semi-analytic
models\cite{WhiteFrenk1991,Kauffmann1993,Cole1994,Baugh1996,Sommerville1999,Kauffmann1999};
the modelling assumptions and parameters are adjusted by trial and
error in order to fit the observed properties of low redshift
galaxies, primarily their joint luminosity-colour distribution and
their distributions of morphology, gas content and central black hole
mass.  Our use of a high-resolution simulation, particularly our
ability to track the evolution of dark matter substructures, removes
much of the uncertainty of the more traditional semi-analytic
approaches based on Monte-Carlo realizations of merger trees.  Our
technique provides accurate positions and peculiar velocities for all
the model galaxies.  It also enables us to follow the evolutionary
history of individual objects and thus to investigate the relationship
between populations seen at different epochs. It is the ability to
establish such evolutionary connections that makes this kind of
modelling so powerful for interpreting observational data.

\subsubsection*{The fate of the first quasars}

Quasars are among the most luminous objects in the Universe and can be
detected at huge cosmological distances.  Their luminosity is thought
to be powered by accretion onto a central, supermassive black
hole. Bright quasars have now been discovered as far back as redshift
$z=6.43$ (ref.~\cite{Fan2003}), and are believed to harbour central
black holes of mass a billion times that of the sun. At redshift
$z\sim 6$, their comoving space density is estimated to be $\sim (2.2
\pm 0.73)\times 10^{-9}\,h^3{\rm Mpc}^{-3}$ (ref.~\cite{Fan2004}).
Whether such extreme rare objects can form at all in a $\Lambda$CDM
cosmology is an open question.

A volume the size of the Millennium Simulation should contain, on
average, just under one quasar at the above space density.  Just what
sort of object should be associated with these ``first quasars'' is,
however, a matter of debate. In the local universe, it appears that
every bright galaxy hosts a supermassive black hole and there is a
remarkably good correlation between the mass of the central black hole
and the stellar mass or velocity dispersion of the bulge of the host
galaxy\cite{Tremaine2002}.  It would therefore seem natural to assume
that at any epoch, the brightest quasars are always hosted by the
largest galaxies. In our simulation, `large galaxies' can be
identified in various ways, for example, according to their dark
matter halo mass, stellar mass, or instantaneous star formation rate.
We have identified the 10 `largest' objects defined in these three
ways at redshift $z=6.2$. It turns out that these criteria all select
essentially the same objects: the 8 largest galaxies by halo mass are
identical to the 8 largest galaxies by stellar mass, only the ranking
differs.  Somewhat larger differences are present when galaxies are
selected by star formation rate, but the 4 first-ranked galaxies are
still amongst the 8 identified according to the other 2 criteria.



\begin{figure}
\vspace*{-1.0cm}\hspace*{-0.3cm}%
\resizebox{8.2cm}{!}{\includegraphics{fig3a.eps}} %
\resizebox{8.2cm}{!}{\includegraphics{fig3b.eps}}\vspace*{0.05cm}\\%
\hspace*{-0.3cm}%
\resizebox{8.2cm}{!}{\includegraphics{fig3c.eps}} %
\resizebox{8.2cm}{!}{\includegraphics{fig3d.eps}}\\%
\caption{Environment of a `first quasar candidate' at high and low redshifts.
  The two panels on the left show the projected dark matter
  distribution in a cube of comoving sidelength $10\,h^{-1}{\rm Mpc}$,
  colour-coded according to density and local dark matter velocity
  dispersion.  The panels on the right show the galaxies of the
  semi-analytic model overlayed on a gray-scale image of the dark
  matter density. The volume of the sphere representing each galaxy is
  proportional to its stellar mass, and the chosen colours encode the
  restframe stellar $B-V$ colour index. While at $z=6.2$ (top) all
  galaxies appear blue due to ongoing star formation, many of the
  galaxies that have fallen into the rich cluster at $z=0$ (bottom)
  have turned red.
\label{FigFirstQuasar}}
\end{figure}


In Figure~\ref{FigFirstQuasar}, we illustrate the environment of a
``first quasar'' candidate in our simulation at $z=6.2$.  The object
lies on one of the most prominent dark matter filaments and is
surrounded by a large number of other, much fainter galaxies. It has a
stellar mass of $6.8\times 10^{10}\,h^{-1}{\rm M}_\odot$, the largest
in the entire simulation at $z=6.2$, a dark matter virial mass of
$3.9\times 10^{12}\,h^{-1}{\rm M}_\odot$, and a star formation rate of
$235\, {\rm M_\odot yr^{-1}}$.  In the local universe central black
hole masses are typically $\sim 1/1000$ of the bulge stellar
mass\cite{Merrit2001}, but in the model we test here these massive
early galaxies have black hole masses in the range $10^8 - 10^9{\rm
M}_\odot$, significantly larger than low redshift galaxies of similar
stellar mass. To attain the observed luminosities, they must convert
infalling mass to radiated energy with a somewhat higher efficiency
than the $\sim 0.1\,c^2$ expected for accretion onto a {\em
non-spinning} black hole.


Within our simulation we can readily address fundamental questions
such as: ``Where are the descendants of the early quasars today?'', or
``What were their progenitors?''. By tracking the merging history
trees of the host halos, we find that all our quasar candidates end up
today as central galaxies in rich clusters. For example, the object
depicted in Fig.~\ref{FigFirstQuasar} lies, today, at the centre of
the ninth most massive cluster in the volume, of mass
$M=1.46\times10^{15}\,h^{-1}{\rm M}_\odot$. The candidate with the
largest virial mass at $z=6.2$ (which has stellar mass $4.7\times
10^{10}\,h^{-1}{\rm M}_\odot$, virial mass $4.85\times
10^{12}\,h^{-1}{\rm M}_\odot$, and star formation rate $218\, {\rm
M_\odot yr^{-1}}$) ends up in the second most massive cluster, of mass
$3.39\times10^{15}\,h^{-1}{\rm M}_\odot$.  Following the merging tree
backwards in time, we can trace our quasar candidate back to redshift
$z=16.7$, when its host halo had a mass of only $1.8\times
10^{10}\,h^{-1}{\rm M}_\odot$. At this epoch, it is one of just 18
objects that we identify as collapsed systems with $\ge 20$
particles. These results confirm the view that rich galaxy clusters
are rather special places.  Not only are they the largest virialised
structures today, they also lie in the regions where the first
structures developed at high redshift. Thus, the best place to search
for the oldest stars in the Universe or for the descendants of the
first supermassive black holes is at the centres of present-day rich
galaxy clusters.


\subsubsection*{The clustering evolution of dark matter and galaxies} 

The combination of a large-volume, high-resolution N-body simulation
with realistic modelling of galaxies enables us to make precise
theoretical predictions for the clustering of galaxies as a function
of redshift and intrinsic galaxy properties. These can be compared
directly with existing and planned surveys.  The 2-point correlation
function of our model galaxies at redshift $z=0$ is plotted in
Fig.~\ref{FigClustering} and is compared with a recent measurement
from the 2dFGRS\cite{Hawkins2003}.  The prediction is remarkably close
to a power-law, confirming with much higher precision the results of
earlier semi-analytic\cite{Kauffmann1999,Benson2000} and
hydrodynamic\cite{Weinberg2004} simulations. This precision will allow
interpretation of the small, but measurable deviations from a pure
power-law found in the most recent data\cite{Padilla2003,Zehavi2004}.
The simple power-law form contrasts with the more complex behaviour
exhibited by the dark matter correlation function but is really no
more than a coincidence.  Correlation functions for galaxy samples
with different selection criteria or at different redshifts do not, in
general, follow power-laws.


\begin{figure}
\vspace*{-0.0cm}\ \\%
\resizebox{14.0cm}{!}{\includegraphics{fig4.eps}}\\
\caption{Galaxy 2-point correlation function at the present epoch.
  Red symbols (with vanishingly small Poisson error-bars) show
  measurements for model galaxies brighter than $M_K = -23$.  Data for
  the large spectroscopic redshift survey 2dFGRS\cite{Hawkins2003} are
  shown as blue diamonds.  The SDSS\cite{Zehavi2002} and
  APM\cite{Padilla2003} surveys give similar results. Both, for the
  observational data and for the simulated galaxies, the correlation
  function is very close to a power-law for $r\le 20\, h^{-1}{\rm
  Mpc}$. By contrast the correlation function for the dark matter
  (dashed line) deviates strongly from a power-law.
\label{FigClustering}}
\end{figure}

 
\begin{figure}
\hspace*{-1.0cm}\ \resizebox{17cm}{!}{\includegraphics{fig5.eps}}
\caption{Galaxy clustering as a function of luminosity and colour.  In
  the panel on the left, we show the 2-point correlation function of
  our galaxy catalogue at $z=0$ split by luminosity in the bJ-band
  (symbols). Brighter galaxies are more strongly clustered, in
  quantitative agreement with observations\cite{Norberg2001} (dashed
  lines).  Splitting galaxies according to colour (right panel), we
  find that red galaxies are more strongly clustered with a steeper
  correlation slope than blue
  galaxies. Observations\cite{Madgwick2003} (dashed lines) show a
  similar trend, although the difference in clustering amplitude is
  smaller than in this particular semi-analytic model.
\label{FigClusteringSubsamples}}
\end{figure}



Although our semi-analytic model was not tuned to match observations
of galaxy clustering, in not only produces the excellent overall
agreement shown in Fig.~\ref{FigClustering}, but also reproduces the
observed dependence of clustering on magnitude and colour in the
2dFGRS and SDSS\cite{Norberg2001,Zehavi2002,Madgwick2003}, as shown in
Figure~\ref{FigClusteringSubsamples}.  The agreement is particularly
good for the dependence of clustering on luminosity. The colour
dependence of the slope is matched precisely, but the amplitude
difference is greater in our model than is
observed\cite{Madgwick2003}.  Note that our predictions for galaxy
correlations split by colour deviate substantially from
power-laws. Such predictions can be easily tested against survey data
in order to clarify the physical processes responsible for the
observed difference.

 
In contrast to the near power-law behaviour of galaxy correlations on
small scales, the large-scale clustering pattern may show interesting
structure.  Coherent oscillations in the primordial plasma give rise
to the well-known acoustic peaks in the
CMB\cite{deBernardis2000,Mauskopf2000,Spergel2003} and also leave an
imprint in the linear power spectrum of the dark matter. Detection of
these ``baryon wiggles'' would not only provide a beautiful
consistency check for the cosmological paradigm, but could also have
important practical applications. The characteristic scale of the
wiggles provides a ``standard ruler'' which may be used to constrain
the equation of state of the dark energy\cite{Blake2003}. A critical
question when designing future surveys is whether these baryon wiggles
are present and are detectable in the {\em galaxy} distribution,
particularly at high redshift.



On large scales and at early times, the mode amplitudes of the {\it
dark matter} power spectrum grow linearly, roughly in proportion to
the cosmological expansion factor. Nonlinear evolution accelerates the
growth on small scales when the dimensionless power $\Delta^2(k) = k^3
P(k)/(2\pi^2)$ approaches unity; this regime can only be studied
accurately using numerical simulations. In the Millennium Simulation,
we are able to determine the nonlinear power spectrum over a larger
range of scales than was possible in earlier work\cite{Jenkins1998},
almost five orders of magnitude in wavenumber $k$.

At the present day, the acoustic oscillations in the matter power
spectrum are expected to fall in the transition region between linear
and nonlinear scales.  In Fig.~\ref{FigWiggles}, we examine the matter
power spectrum in our simulation in the region of the
oscillations. Dividing by the smooth power spectrum of a $\Lambda$CDM
model with no baryons\cite{Bardeen1986} highlights the baryonic
features in the initial power spectrum of the simulation, although
there is substantial scatter due to the small number of large-scale
modes.  Since linear growth preserves the relative mode amplitudes, we
can approximately correct for this scatter by scaling the measured
power in each bin by a multiplicative factor based on the initial
difference between the actual bin power and the mean power expected in
our $\Lambda$CDM model.  This makes the effects of nonlinear evolution
on the baryon oscillations more clearly visible.  As
Fig.~\ref{FigWiggles} shows, nonlinear evolution not only accelerates
growth but also reduces the baryon oscillations: scales near peaks
grow slightly more slowly than scales near troughs. This is a
consequence of the mode-mode coupling characteristic of nonlinear
growth. In spite of these effects, the first two ``acoustic peaks''
(at $k\sim 0.07$ and $k\sim 0.13\,h\,{\rm Mpc}^{-1}$, respectively) in
the dark matter distribution do survive in distorted form until the
present day (see the lower right panel of Fig.~\ref{FigWiggles}).


\begin{figure}
\begin{center}
\vspace*{-1.6cm}\hspace*{-1.0cm}%
\resizebox{15.8cm}{!}{\includegraphics{fig6.eps}}\vspace*{-1.0cm}%
\end{center}
\caption{ Power spectra of the dark matter and galaxy distributions in
  the baryon oscillation region. All measurements have been divided by
  a linearly evolved, CDM-only power spectrum\cite{Bardeen1986}. Red
  circles show the dark matter, and green squares the galaxies. Blue
  symbols give the actual realization of the initial fluctuations in
  our simulation, which scatters around the mean input power (black
  lines) due to the finite number of modes.  Since linear growth
  preserves relative mode amplitudes, we correct the power in each bin
  to the expected input power and apply these scaling factors at all
  other times.  At $z=3.06$, galaxies with stellar mass above
  $5.83\times 10^9\,h^{-1}{\rm M}_\odot$ and space-density of $8\times
  10^{-3}\,h^{3}{\rm Mpc}^{-3}$ were selected. Their large-scale
  density field is biased by a factor $b=2.7$ with respect to the dark
  matter (the galaxy measurement has been divided by $b^2$). At $z=0$,
  galaxies brighter than $M_B = -17$ and a space density higher by a
  factor $\sim 7.2$ were selected. They exhibit a slight antibias,
  $b=0.92$.  Corresponding numbers for $z=0.98$ are $M_B = -19$ and
  $b=1.15$.
\label{FigWiggles}}
\end{figure}


Are the baryon wiggles also present in the galaxy distribution?
Fig.~\ref{FigWiggles} shows that the answer to this important question
is `yes'. The $z=0$ panel shows the power spectrum for all model
galaxies brighter than $M_B = -17$. On the largest scales, the galaxy
power spectrum has the same shape as that of the dark matter, but with
slightly lower amplitude corresponding to an ``antibias'' of
8\%. Samples of brighter galaxies show less antibias while for the
brightest galaxies, the bias becomes slightly positive. The figure
also shows measurements of the power spectrum of luminous galaxies at
redshifts $z=0.98$ and $z=3.06$. Galaxies at $z=0.98$ were selected to
have a magnitude $M_B<-19$ in the restframe, whereas galaxies at
$z=3.06$ were selected to have stellar mass larger than $5.83\times
10^9\,h^{-1}{\rm M}_\odot$, corresponding to a space density of
$8\times 10^{-3}\,h^{3}{\rm Mpc}^{-3}$, similar to that of the
Lyman-break galaxies observed at $z\sim 3$\cite{Adelberger1998}.
Signatures of the first two acoustic peaks are clearly visible at both
redshifts, even though the density field of the $z=3$ galaxies is much
more strongly biased with respect to the dark matter (by a factor
$b=2.7$) than at low redshift.  Selecting galaxies by their star
formation rate rather than their stellar mass (above $10.6\,{\rm
M_\odot yr^{-1}}$ for an equal space density at $z=3$) produces very
similar results.

Our analysis demonstrates conclusively that baryon wiggles should
indeed be present in the galaxy distribution out to redshift
$z=3$. This has been assumed but not justified in recent proposals to
use evolution of the large-scale galaxy distribution to constrain the
nature of the dark energy.  To establish whether the baryon
oscillations can be measured in practice with the requisite accuracy
will require detailed modelling of the selection criteria of an actual
survey and a thorough understanding of the systematic effects that
will inevitably be present in real data. These issues can only be
properly addressed by means of specially designed mock catalogues
constructed from realistic simulations. We plan to construct suitable
mock catalogues from the Millennium Simulation and make them publicly
available. Our provisional conclusion, however, is that the next
generation of galaxy surveys offers excellent prospects for
constraining the equation of state of the dark energy.


N-body simulations of CDM universes are now of such size and quality
that realistic modelling of galaxy formation in volumes matched to
modern surveys has become possible. Detailed studies of galaxy and AGN
evolution exploiting the unique dataset of the Millennium Simulation
therefore enable stringent new tests of the theory of hierarchical
galaxy formation.  Using the simulation we demonstrated that quasars
can plausibly form sufficiently early in a $\Lambda$CDM universe to be
compatible with observation, that their progenitors were already
massive by $z \sim 16$, and that their $z=0$ descendents lie at the
centres of cD galaxies in rich galaxy clusters. Interesting tests of
our predictions will become possible if observations of the black hole
demographics can be extended to high redshift, allowing, for example,
a measurement of the evolution of the relationship between
supermassive black hole masses and the velocity dispersion of their
host stellar bulges.

We have also demonstrated that a power-law galaxy autocorrelation
function can arise naturally in a $\Lambda$CDM universe, but that this
suggestively simple behaviour is merely a coincidence. Galaxy surveys
will soon reach sufficient statistical power to measure precise
deviations from power-laws for galaxy subsamples, and we expect that
comparisons of the kind we have illustrated will lead to tight
constraints on the physical processes included in the galaxy formation
modelling.  Finally, we have demonstrated for the first time that the
baryon-induced oscillations recently detected in the CMB power
spectrum should survive in distorted form not only in the nonlinear
dark matter power spectrum at low redshift, but also in the power
spectra of realistically selected galaxy samples at $0<z<3$. Present
galaxy surveys are marginally able to detect the baryonic features at
low redshifts\cite{Cole2005,Eisenstein2005}. If future surveys improve
on this and reach sufficient volume and galaxy density also at high
redshift, then precision measurements of galaxy clustering will shed
light on one of the most puzzling components of the universe, the
elusive dark energy field.



\vspace*{1cm}

\section*{Methods}


The Millennium Simulation was carried out with a specially customised
version of the {\small GADGET2}~(Ref. \cite{Springel2001})~code, using
the ``TreePM'' method\cite{Xu1995} for evaluating gravitational
forces. This is a combination of a hierarchical multipole expansion,
or ``tree'' algorithm\cite{Barnes1986}, and a classical, Fourier
transform particle-mesh method\cite{Hockney1981}. The calculation was
performed on 512 processors of an IBM p690 parallel computer at the
Computing Centre of the Max-Planck Society in Garching, Germany. It
utilised almost all the 1~TB of physically distributed memory
available.  It required about $350\,000$ processor hours of CPU time,
or 28 days of wall-clock time.  The mean sustained floating point
performance (as measured by hardware counters) was about 0.2~TFlops,
so the total number of floating point operations carried out was of
order $5\times 10^{17}$.

The cosmological parameters of our $\Lambda$CDM-simulation are:
$\Omega_{\rm m}= \Omega_{\rm dm}+\Omega_{\rm b}=0.25$, $\Omega_{\rm
b}=0.045$, $h=0.73$, $\Omega_\Lambda=0.75$, $n=1$, and
$\sigma_8=0.9$. Here $\Omega_{\rm m}$ denotes the total matter density
in units of the critical density for closure, $\rho_{\rm crit}=3
H_0^2/(8\pi G)$. Similarly, $\Omega_{\rm b}$ and $\Omega_\Lambda$
denote the densities of baryons and dark energy at the present
day. The Hubble constant is parameterised as $H_0 = 100\, h\, {\rm
km\, s^{-1} Mpc^{-1}}$, while $\sigma_8$ is the {\em rms} linear mass
fluctuation within a sphere of radius $8\, h^{-1}{\rm Mpc}$
extrapolated to $z=0$.  Our adopted parameter values are consistent
with a combined analysis of the 2dFGRS\cite{Colless2001} and first
year WMAP data\cite{Spergel2003}.


The simulation volume is a periodic box of size $500\,h^{-1}{\rm Mpc}$
and individual particles have a mass of $8.6\times 10^8\,h^{-1}{\rm
M}_{\odot}$.  This volume is large enough to include interesting rare
objects, but still small enough that the halos of all luminous
galaxies brighter than $0.1 L_\star$ are resolved with at least 100
particles. At the present day, the richest clusters of galaxies
contain about 3 million particles. The gravitational force law is
softened isotropically on a comoving scale of $5\,h^{-1}{\rm kpc}$
(Plummer-equivalent), which may be taken as the spatial resolution
limit of the calculation. Thus, our simulation achieves a dynamic
range of $10^5$ in 3D, and this resolution is available everywhere in
the simulation volume.

Initial conditions were laid down by perturbing a homogeneous,
`glass-like', particle distribution\cite{White1996} with a realization
of a Gaussian random field with the $\Lambda$CDM linear power spectrum
as given by the Boltzmann code {\small CMBFAST}\cite{Seljak1996}. The
displacement field in Fourier space was constructed using the
Zel'dovich approximation, with the amplitude of each random phase mode
drawn from a Rayleigh distribution. The simulation started at redshift
$z=127$ and was evolved to the present using a leapfrog integration
scheme with individual and adaptive timesteps, with up to $11\,000$
timesteps for individual particles.  We stored the full particle data
at 64 output times, each of size 300 GB, giving a raw data volume of
nearly 20 TB.  This allowed the construction of finely resolved
hierarchical merging trees for tens of millions of halos and for the
subhalos that survive within them. A galaxy catalogue for the full
simulation, typically containing $\sim2\times 10^6$ galaxies at $z=0$
together with their full histories, can then be built for any desired
semi-analytic model in a few hours on a high-end workstation.



The semi-analytic model itself can be viewed as a simplified
simulation of the galaxy formation process, where the star formation
and its regulation by feedback processes is parameterised in terms of
simple analytic physical models. These models take the form of
differential equations for the time evolution of the galaxies that
populate each hierarchical merging tree.  In brief, these equations
describe radiative cooling of gas, star formation, growth of
supermassive black holes, feedback processes by supernovae and AGN,
and effects due to a reionising UV background. In addition, the
morphological transformation of galaxies and the process of metal
enrichment are modelled as well.  To make direct contact with
observational data, we apply modern population synthesis models to
predict spectra and magnitudes for the stellar light emitted by
galaxies, also including simplified models for dust obscuration. In
this way we can match the passbands commonly used in observations.

The basic elements of galaxy formation modelling follow previous
studies\cite{WhiteFrenk1991,Kauffmann1993,Cole1994,Baugh1996,Sommerville1999,Kauffmann1999,Springel2001b}
(see also Supplementary Information), but we also use novel approaches
in a number of areas. Of substantial importance is our tracking of
dark matter substructure. This we carry out consistently and with
unprecedented resolution throughout our large cosmological volume,
allowing an accurate determination of the orbits of galaxies within
larger structures, as well as robust estimates of the survival time of
structures infalling into larger objects. Also, we use dark matter
substructure properties, like angular momentum or density profile, to
directly determine sizes of galactic disks and their rotation curves.
Secondly, we employ a novel model for the build-up of a population of
supermassive black holes in the universe. To this end we extend the
quasar model developed in previous work\cite{Kauffmann2000} with a
`radio mode', which describes the feedback activity of central AGN in
groups and clusters of galaxies. While largely unimportant for the
cumulative growth of the total black hole mass density in the
universe, our results show that the radio mode becomes important at
low redshift, where it has a strong impact on cluster cooling
flows. As a result, it reduces the brightness of central cluster
galaxies, an effect that shapes the bright end of the galaxy
luminosity function, bringing our predictions into good agreement with
observation.



\bibliography{main}


\paragraph*{Supplementary Information} accompanies the paper on 
{\bf www.nature.com/nature}.

\vspace*{-0.5cm}\paragraph*{Acknowledgements} We would like to thank the anonymous
referees who helped to improve the paper substantially. The
computations reported here were performed at the {\em Rechenzentrum der
Max-Planck-Gesellschaft} in Garching, Germany.


\vspace*{-0.5cm}\paragraph*{Competing interests} The authors declare that they have no
competing financial interests.


\vspace*{-0.5cm}\paragraph*{Correspondence} and requests for materials should
be addressed to V.S.~(email: vspringel@mpa-garching.mpg.de).