text:
# The Two-Dimensional Bose-Einstein Condensate

## I Introduction

In a recent paper gorlitz , one of the groups that pioneered the formation and detection of Bose-Einstein condensation (BEC) in harmonically trapped atomic gases Wieman ; Ketterle ; Hulet  reports the creation of (pseudo-) two-dimensional condensates. These have been produced by taking a three-dimensional condensate of 23Na atoms and carrying out two independent processes on it: 1) Initially, one of the confining frequencies (that in the $`z`$ direction, $`\omega _z`$, in order to minimize the effects of gravity) is increased until the condensate radius in that dimension is smaller than the healing length associated with the interaction between atoms (taken here to be repulsive and parametrized by the two-body scattering length $`a`$). This is not sufficient to reduce the dimensionality, however, since the atoms, each of which has mass $`m`$, will literally squeeze into the third dimension if there are more than

$$N=\sqrt{\frac{32\mathrm{}}{225ma^2}}\sqrt{\frac{\omega _z^3}{\omega _x^2\omega _y^2}}$$
(1)
of them in the trap. 2) Consequently, the number of atoms in the condensate must be reduced; this is achieved by exposing the condensate to a thermal beam. The reduction in effective dimensionality becomes apparent when the aspect ratio of the condensate, which is independent of $`N`$ in 3D, starts to change as the number of atoms is gradually reduced. The condensates thus produced have a number of atoms that ranges between $`10^4`$ and $`10^5`$.

This constitutes an important experimental contribution to the long-standing debate about the existence of BEC in two dimensions. It is a well-known fact (and a standard textbook exercise) that BEC cannot happen in a 2D homogeneous ideal gas; a rigorous mathematical theorem HT  extends this result to the case where there are interactions between the bosons. When the system is in a harmonic trap, on the other hand, BEC can occur in two dimensions bagklep2  below a temperature $`k_\text{B}T_c=\mathrm{}\omega \sqrt{6N/\pi ^2}`$, but the theorem is once again valid when interactions between the bosons are considered mulLong : while there is a BEC in the 2D system, it occurs at $`T=0`$.

The preceding discussion is valid only in the thermodynamic limit, which in the particular case of an isotropically trapped system consists of making $`N\mathrm{}`$ and $`\omega 0`$ in such a way that $`N\omega ^2`$ remains finite Damle ; mulLong . The question remains whether a phenomenon resembling BEC—that is, the accumulation of a macroscopic number of particles in a single quantum state—occurs or not when the system consists of a *finite* number of particles confined by a trap of finite frequency, as is certainly the case in experimental situations. If there is such a phenomenon, one would like to know more about the process by which this “condensate” is destabilized at finite temperatures as $`N`$ grows.

Some authors kagan87 ; kagan00  have considered, in the finite homogeneous 2D case, the possibility of a BEC. A similar analysis petrov  was considered for a quasi-2D trapped gas; the latter reference finds that the phase fluctuations in the condensate vary with temperature and particle number as

$$(\delta \varphi )^2T\mathrm{log}N$$
(2)
which diverges for finite temperature as $`N\mathrm{}`$, as one would expect. For finite $`N`$, on the other hand, the fluctuations are tempered at very low temperature; since the coherence length, though finite, is still larger than the characteristic length imposed on the system by the trap filshev  or by walls, one can speak of a “quasicondensate.” It is this quasicondensation that we wish to study in this paper. Reference jakkola  reports the observation of a quasicondensate in a homogeneous system of atomic hydrogen adsorbed on 4He; at this point we cannot tell if the condensates reported in LABEL:gorlitz are actually quasicondensates.

One can approach the study of the 2D Bose gas by employing mean-field theory, which in this context refers to the Hartree-Fock-Bogoliubov (HFB) equations and various simplifications thereof that we will review in quantitative detail and classify in the present work. The HFB theory has been remarkably successful in treating the 3D case; though we give a few representative references below, we refer the reader to LABEL:RMPBEC, which includes a comprehensive review and an exhaustive list of references, and will henceforth concentrate on the work that has been carried out in 2D.

The HFB theory assumes from the outset that the system is partially condensed, and proposes separate equations to describe the condensate and the uncondensed (thermal) component. The condensate is described by a macroscopic wavefunction that obeys a generalized Gross-Pitaevskiĭ (GP) equation Gross ; Pita , while the noncondensate consists of a superposition of Bogoliubov quasiparticle and “quasihole” excitations weighted by a Bose distribution. The expression for the noncondensate can be simplified by neglecting the quasihole excitations (the Hartree-Fock scheme) HuseSiggia ; BagKlep , by performing a semiclassical WKB approximation VdB , or by combining both of these GPSTT ; GPS ; HKN . Each one of the above schemes can be further simplified, when the thermodynamic limit is approached, by neglecting the kinetic energy of the condensate: this is the Thomas-Fermi limit GSL ; BaymPeth . Finally one can neglect the interactions between thermal atoms, arriving at the semi-ideal model Stamper ; Minguzzi .

The semi-ideal model has been implemented in 2D Turks ; kim1 , but has been found to yield unphysical results as the interaction strength becomes sizable. The full-blown HFB model does not appear to be much more successful, at least when considered in the semiclassical limit. In a previous paper mulShort  one of us found that, below a certain temperature, the introduction of a condensate in the Thomas-Fermi limit (corresponding to the thermodynamic limit) renders the HFB equations incompatible, with the noncondensate density becoming infinite at every point in space. The singularity occurs at the low end of the energy spectrum, indicating that the condensate is being destabilized by long-wavelength phonons; this interpretation in terms of phase fluctuations had already been proposed for the homogeneous system HT ; Reatto . In this work we also report our inability to find self-consistent semiclassical solutions to the HFB model when a finite trapped system is considered.

Moreover, it was recently discovered bhaduri ; brandon  that it is possible, in the 2D case exclusively, to solve the HFB equations semiclassically at *any* temperature without even having to invoke the presence of a condensate (thus obtaining what we will call an “uncondensed” solution). In other words, it is possible to simply cross out the condensate component and solve for the system to a temperature close to $`T=0`$. The solution thus obtained shows an accumulation of atoms at the center of the trap and yields a bulge in density similar to that caused by the presence of a condensate, even though no state is macroscopically occupied.

These results appear to reinforce the conclusion that BEC cannot happen in the two-dimensional trapped system. However, we are still confronted with the experimental results described above. Furthermore, two different Monte Carlo simulations Heinrichs ; Pearson  show significant concentrations of particles in the lowest energy state for finite $`N`$, though it must be said that this method provides little information about the types of excitations that contribute to the disappearance of the condensate, and that it is difficult to carry out such simulations on very large systems.

The HFB model cannot be discounted at this stage either: when we restrict ourselves to the Hartree-Fock approximation, it is possible to find self-consistent solutions (henceforth referred to as “condensed”) involving a condensate, both for finite systems and in the thermodynamic limit, when using either the discrete set of equations kim2  or the WKB approximation Ibiza  to treat the noncondensate. The Hartree-Fock approximation neglects phononlike excitations, so it is not surprising that it yields solutions. However, one may be able to justify its usefulness.

It is known that, in the infinite homogeneous system, infrared singularities occur but are renormalized by interactions, providing, in essence, a cutoff at a low wavenumber $`k_0\sqrt{nmU}`$ fishhoh ; prokofio ; popovbook , where $`n`$ is the density, $`m`$ the particle mass, and $`U`$ the effective interaction strength. Indeed, it is possible to estimate the Berezinskii-Kosterlitz-Thouless (BKT) transition temperature by simply cutting off the ideal-gas density expression at this $`k_0`$ prokofio . Presumably a similar situation occurs in the trapped case. The Hartree-Fock approach provides a convergent theory by cutting off the singularities at a wavenumber similar to that of more rigorous theories. Whether such an approach gives a reasonable estimate of the BKT transition temperature, the superfluid density, or a quasicondensate density for the finite system will need to await a more rigorous theoretical approach to the interacting 2D trapped gas Stoofnew .

Given the above limitations, we analyze the character of the BEC in the 2D trapped system by solving the coupled equations of the theory. We find that, in the Hartree-Fock scheme, it is possible to find both condensed and uncondensed solutions for the two-dimensional equations. We also calculate the free energy corresponding to each one and find that the condensed solution has a lower free energy at all temperatures, which appears to imply that, at least at this level of approximation, the uncondensed solution is unphysical or metastable. The condensed solution will be “preferred” over the uncondensed one, and we assume that the solution represents an approximation to a quasi-condensate.

It is evident from our discussion that our approach to BEC in 2D trapped systems is a preliminary one and that a further analysis that takes fully into account the BKT transition is necessary. A start to answering this need has been presented in LABEL:Stoofnew, and we intend to return to this problem ourselves in the future.

## II The Model

### II.1 The Hartree-Fock-Bogoliubov equations

Throughout this paper we use a dimensionless system of units in which all lengths are scaled by the oscillator length $`x_0\sqrt{\mathrm{}/m\omega }`$ and all energies are expressed in terms of the one-dimensional ground-state energy of the oscillator, $`\mathrm{}\omega /2`$. Dimensionless variables will in general carry a tilde: for example, the total density $`n`$ becomes $`\stackrel{~}{n}x_0^2n`$ and the chemical potential $`\stackrel{~}{\mu }\mu /\frac{1}{2}\mathrm{}\omega `$.

The HFB equations RMPBEC ; GPS ; Fetter72 ; Griffin96  result from assuming that 1) the (repulsive) interactions between atoms consist exclusively of two-body low-energy collisions that can be described by a delta-function pseudopotential HuangYang  of strength $`g`$ (related to its dimensionless counterpart $`\stackrel{~}{\gamma }`$ through $`g\frac{1}{2}\mathrm{}\omega x_0^2\stackrel{~}{\gamma }`$), that 2) the many-body field operator $`\mathrm{\Psi }`$ can be decomposed via

$$\mathrm{\Psi }=\mathrm{\Psi }+\stackrel{~}{\psi }\stackrel{~}{\mathrm{\Phi }}+\stackrel{~}{\psi },$$
(3)
where the ensemble average $`\stackrel{~}{\mathrm{\Phi }}`$ is a real macroscopic wavefunction that describes the condensate (reflecting the imposition of macroscopic long-range order on the condensed system), and that 3) products of noncondensate operators can be simplified using a finite-temperature version of Wick’s theorem FetWal . If we insert (3) into the many-body grand-canonical Hamiltonian

$$\stackrel{~}{H}=d^2x\mathrm{\Psi }^{}(\stackrel{~}{\mathrm{\Lambda }}+\frac{\stackrel{~}{\gamma }}{2}\mathrm{\Psi }^{}\mathrm{\Psi }^{}\mathrm{\Psi })\mathrm{\Psi },$$
(4)
where $`\stackrel{~}{\mathrm{\Lambda }}=\stackrel{~}{}^2+x^2\stackrel{~}{\mu }`$, and neglect anomalous averages via the Popov approximation GPS , we obtain an expression that can be diagonalized and yields an infinite set of coupled differential equations.

On one hand, the macroscopic wavefunction mentioned above is the square root of the dimensionless condensate density $`n_0`$ and obeys the generalized Gross-Pitaevskiĭ equation,

$$\stackrel{~}{\mathrm{\Lambda }}\stackrel{~}{\mathrm{\Phi }}+\stackrel{~}{\gamma }(\stackrel{~}{n}_0+2\stackrel{~}{n}^{})\stackrel{~}{\mathrm{\Phi }}=0,$$
(5)
where $`\stackrel{~}{n}^{}\stackrel{~}{n}\stackrel{~}{n}_0`$ is the noncondensate density. The factor of $`2`$ in (5) and hereafter is a consequence of the direct (Hartree) and exchange (Fock) terms being identical, which follows from the fact that the delta-function interaction that we are considering has zero range Griffin96 . (The term involving the condensate does not include that factor; this is a consequence of the restricted grand-canonical ensemble that we are using. See the end of the next section.)

The noncondensate, in turn, is described by an infinite number of pairs of functions that obey Fetter72

$$\left(\begin{array}{cc}\stackrel{~}{\mathrm{\Lambda }}+2\stackrel{~}{\gamma }\stackrel{~}{n}& \stackrel{~}{\gamma }\stackrel{~}{n}_0\\ \stackrel{~}{\gamma }\stackrel{~}{n}_0& \stackrel{~}{\mathrm{\Lambda }}+2\stackrel{~}{\gamma }\stackrel{~}{n}\end{array}\right)\left(\begin{array}{cc}u_j& \\ v_j& \end{array}\right)=\stackrel{~}{ϵ}_j\left(\begin{array}{cc}u_j& \\ v_j& \end{array}\right)$$
(6)
and which generate the noncondensate density via

$$\stackrel{~}{n}^{}(𝐱)=\underset{j}{}\left(\left(\left|u_j\right|^2+\left|v_j\right|^2\right)f_j+\left|v_j\right|^2\right),$$
(7)
where we have introduced the Bose-Einstein distribution factor $`f_j(e^{\stackrel{~}{ϵ}_j/t}1)^1`$, which appears when the free energy of the system is minimized, and the dimensionless temperature $`tk_\text{B}T/\frac{1}{2}\mathrm{}\omega `$. The last term of (7) describes the zero-temperature depletion of the condensate, which accounts for less than 1% of the particles and is therefore negligible GPS . This self-consistent set of equations is closed, and the chemical potential found, by imposing that the total number of particles remain fixed:

$$N=N_0+N^{}=d^2x\stackrel{~}{n}(𝐱)=d^2x\left(\stackrel{~}{n}_0+\stackrel{~}{n}^{}\right).$$
(8)
At temperatures such that $`k_\text{B}T\mathrm{}\omega `$, one can use the semiclassical (WKB) approximation VdB  that results from taking $`u_ju(𝐱)e^{i\varphi }`$ and $`v_jv(𝐱)e^{i\varphi }`$. The phase common to both defines a quasiparticle momentum through $`𝜿=\stackrel{~}{\mathbf{}}\varphi `$; $`u`$$`v`$, and $`𝜿`$ vary sufficiently slowly with $`𝐱`$ that their spatial derivatives can be neglected GPS . The distribution factor is now $`f_jf(𝐱,𝜿)`$, infinite sums are transformed into momentum integrals—which in two dimensions can be solved in closed form mulShort —and the equations become algebraic, yielding the Bogoliubov energy spectrum GPS

$`\stackrel{~}{ϵ}_{\text{HFB}}(𝜿,𝐱)=\sqrt{(\kappa ^2+x^2+2\stackrel{~}{\gamma }\stackrel{~}{n}\stackrel{~}{\mu })^2\stackrel{~}{\gamma }^2\stackrel{~}{n}_0^2}`$ (9)
and the following integral expression for the noncondensate density:

$`\stackrel{~}{n}^{}(𝐱)`$ $`=`$ $`{\displaystyle \frac{t}{4\pi }}{\displaystyle _{((x^2+2\stackrel{~}{\gamma }\stackrel{~}{n}\stackrel{~}{\mu })^2\stackrel{~}{\gamma }^2\stackrel{~}{n}_0^2)^{1/2}/t}^{\mathrm{}}}{\displaystyle \frac{d\xi }{e^\xi 1}}`$ (10)
$`=`$ $`{\displaystyle \frac{t}{4\pi }}\mathrm{log}\left(1e^{(\sqrt{(x^2+2\stackrel{~}{\gamma }\stackrel{~}{n}\stackrel{~}{\mu })^2\stackrel{~}{\gamma }^2\stackrel{~}{n}_0^2})/t}\right).`$
At high enough temperatures there is no condensate, and the density on the left-hand side of (10) is just the total density. We thus have to solve a single self-consistent equation,

$$\stackrel{~}{n}^{}(𝐱)\stackrel{~}{n}(𝐱)=\frac{t}{4\pi }\mathrm{log}\left(1e^{(x^2+2\stackrel{~}{\gamma }\stackrel{~}{n}\stackrel{~}{\mu })/t}\right);$$
(11)
the chemical potential is once again calculated by requiring $`N`$ to be fixed. Equation (11) has been found by the authors of LABEL:bhaduri to be soluble at *all temperatures*, a result that we confirm in the present work. Reference mulShort  had mistakenly concluded that it was impossible to solve (11) below a certain temperature.

The Hartree-Fock approximation GSL ; GPS ; HuseSiggia  amounts to neglecting the $`v_j`$ in (6) and, when combined with the semiclassical treatment, results in the disappearance of the last term within the square root of Eq. (9). The energy spectrum now becomes

$$\stackrel{~}{ϵ}_{\text{HF}}\stackrel{~}{ϵ}=\kappa ^2+x^2+2\stackrel{~}{\gamma }\stackrel{~}{n}\stackrel{~}{\mu }$$
(12)
and the noncondensate density turns into

$$\stackrel{~}{n}^{}(𝐱)=\frac{t}{4\pi }\mathrm{log}\left(1e^{(x^2+2\stackrel{~}{\gamma }\stackrel{~}{n}\stackrel{~}{\mu })/t}\right),$$
(13)
an expression that differs from (11) in that the left-hand side corresponds only to the noncondensate density. The three-dimensional version of this equation, coupled with (5) and (8), has been frequently used in the literature to study the three-dimensional gas. Reference HKN , for example, exhibits a detailed comparison of its predictions to those of Monte Carlo simulations and finds excellent agreement between the two.

In two dimensions, the Hartree-Fock equations have been solved without resorting to the WKB approximation kim2 . The authors of this reference succeeded in finding self-consistent density profiles and used them to study the temperature dependence of the condensate fraction. Our results Ibiza , which do take advantage of the semiclassical approximation, agree quite well with theirs (see Fig. 1).

### II.2 The Thermodynamic Limit

When $`N`$ is large we can neglect the kinetic energy of the system, which in 2D can be shown to be of order $`1/N`$ mulShort , and obtain the Thomas-Fermi approximation BaymPeth

$$\stackrel{~}{\gamma }\stackrel{~}{n}_0=(\stackrel{~}{\mu }x^22\stackrel{~}{\gamma }\stackrel{~}{n}^{})\mathrm{\Theta }(\stackrel{~}{\mu }x^22\stackrel{~}{\gamma }\stackrel{~}{n}^{}),$$
(14)
where $`\mathrm{\Theta }(x)`$ is the Heaviside step function, introduced to ensure that the density profile is everywhere real and positive. Also, since the thermodynamic limit requires that $`\omega 0`$, the WKB approximation becomes rigorous and can be used with confidence. Thus we can insert expression (14) in the Bogoliubov energy spectrum (9) and show that the latter reduces to

$$\stackrel{~}{ϵ}_{\text{HFB}}(𝜿,𝐱)\kappa \sqrt{\kappa ^2+2\stackrel{~}{\gamma }\stackrel{~}{n}_0}\kappa \sqrt{2\stackrel{~}{\gamma }\stackrel{~}{n}_0}$$
(15)
for small quasimomenta. The fact that it is linear clearly shows us that, in this approximation, the low end of the energy spectrum corresponds to phononlike quasiparticles. Now, if we introduce (14) into the noncondensate density (10), we can see that the argument of the logarithm vanishes at all temperatures, making the density diverge at every point in space. The first line of (10) shows that the divergence in the integral is caused at its lower limit; this restates the conclusion arrived at in LABEL:mulShort: low-energy phonons destabilize the condensate in the two-dimensional thermodynamic limit when the noncondensate quasiparticles obey the Bogoliubov spectrum.

In the Hartree-Fock approximation, on the other hand, the energy spectrum tends in this limit to $`\stackrel{~}{ϵ}\kappa ^2+\stackrel{~}{\gamma }\stackrel{~}{n}_0`$ and predicts single-particle excitations whose minimum energy is $`\stackrel{~}{\gamma }\stackrel{~}{n}_0`$; this can be interpreted equivalently by assigning a minimum value $`\kappa _c^2=\stackrel{~}{\gamma }\stackrel{~}{n}_0`$ for the excitation quasimomentum. This cutoff is consistent with those proposed in the past fishhoh ; kagan87 ; kagan00 ; prokofio ; popovbook  and removes the infrared singularity in the HFB equations.

This momentum cutoff is robust enough that it enables one to carry out Hartree-Fock calculations even in the thermodynamic limit: in fact, as was first found in LABEL:GSL, the introduction of this limit actually simplifies the calculations, and it is possible to find self-consistent solutions by simultaneously treating the noncondensate in the Hartree-Fock approximation and the condensate in the Thomas-Fermi limit kim2 ; Ibiza . This model cannot provide realistic density profiles at every point in space, since the condensate density (14) has a discontinuous derivative at its edge, but predicts quite reasonable results outside of this region, as can be seen in Fig. 1.

### II.3 The Free Energy

We have seen that in 2D the mean-field-theory BEC equations admit solutions both with and without a condensate. The unphysical solution should be that with the highest free energy, since equilibrium at finite temperatures occurs when the grand potential attains a minimum; in fact, the Bose-Einstein distribution factor in Eq. (7) comes from minimizing this quantity HuseSiggia ; deGennes , which in our dimensionless units adopts the form

$$\mathrm{\Omega }=(U\mu NTS)/\frac{1}{2}\mathrm{}\omega d^2x(\mathrm{{\rm Y}}\stackrel{~}{\mu }\stackrel{~}{n}t\mathrm{\Sigma }).$$
(16)
At this point we have to keep in mind that the grand-canonical free energy is a function of $`\stackrel{~}{\mu }`$, not of $`N`$; in order to make a meaningful comparison of these energies at the same $`N`$, then, we have to compare the Helmholtz free energies, given by $`\stackrel{~}{A}=\mathrm{\Omega }+\stackrel{~}{\mu }N`$. The expressions given below have all been derived in the grand-canonical ensemble and thus contain the chemical potential; we will eliminate this dependence on $`\stackrel{~}{\mu }`$ by adding $`\stackrel{~}{\mu }N`$ to the expressions we obtain.

In the Hartree-Fock approximation, the grand-canonical energy density of the system is given by GPS

$$\mathrm{{\rm Y}}\stackrel{~}{\mu }\stackrel{~}{n}=\stackrel{~}{\mathrm{\Phi }}(\stackrel{~}{\mathrm{\Lambda }}+\frac{\stackrel{~}{\gamma }}{2}\stackrel{~}{n}_0)\stackrel{~}{\mathrm{\Phi }}+\frac{1}{(2\pi )^2}\frac{d^2\kappa \stackrel{~}{ϵ}}{e^{\stackrel{~}{ϵ}/t}1}\stackrel{~}{\gamma }\stackrel{~}{n}^2;$$
(17)
the first term corresponds to the condensate energy, while the second one is the sum, weighted by the Bose-Einstein distribution, of the energies of the excited states; this last expression includes an extra term $`\stackrel{~}{\gamma }\stackrel{~}{n}^2`$ that has to be subtracted explicitly, as is usually the case in Hartree-Fock calculations FetWal . We can simplify (17) further by invoking Eq. (5):

$$\mathrm{{\rm Y}}\stackrel{~}{\mu }\stackrel{~}{n}=\stackrel{~}{\gamma }\stackrel{~}{n}^2+\frac{\stackrel{~}{\gamma }}{2}\stackrel{~}{n}_0^2+\frac{1}{(2\pi )^2}\frac{d^2\kappa \stackrel{~}{ϵ}}{e^{\stackrel{~}{ϵ}/t}1}.$$
(18)
The entropy of the system can be found from the combinatorial expression LanLif ; FetWal

$$S=k_\text{B}\underset{i}{}(f_i\mathrm{log}f_i(f_i+1)\mathrm{log}(f_i+1)),$$
(19)
which in the WKB and Hartree-Fock approximations yields the entropy density GPS

$$t\mathrm{\Sigma }=\frac{1}{(2\pi )^2}\frac{d^2\kappa \stackrel{~}{ϵ}}{e^{\stackrel{~}{ϵ}/t}1}\frac{t}{(2\pi )^2}d^2\kappa \mathrm{log}(1e^{\stackrel{~}{ϵ}/t}).$$
(20)
The first term in (20) cancels with the last one in (18), and the other term can be integrated in closed form, yielding

$`\stackrel{~}{A}_\mathrm{c}`$ $`=`$ $`\stackrel{~}{\mu }N{\displaystyle }d^2x(\stackrel{~}{\gamma }(\stackrel{~}{n}^2{\displaystyle \frac{1}{2}}\stackrel{~}{n}_0^2)`$ (21)
$`+{\displaystyle \frac{t^2}{4\pi }}g_2(e^{(x^2+2\stackrel{~}{\gamma }\stackrel{~}{n}\stackrel{~}{\mu })/t}))`$
for the free energy of the condensed solution. (Here we have introduced the Bose-Einstein integral $`g_\sigma (x)_{k=1}^{\mathrm{}}x^k/k^\sigma `$. Note that $`g_1(x)=\mathrm{log}(1x)`$, a function that repeatedly appears above.) As $`T0`$, this expression reduces to the correct value in the homogeneous case BagKlep . In the thermodynamic limit, we can find an expression for the free energy per particle at $`T=0`$ using the Thomas-Fermi density profile:

$$\frac{\stackrel{~}{A}_{\text{TF}}}{N}=\frac{2}{3}\sqrt{\frac{2}{\pi }}\sqrt{N\stackrel{~}{\gamma }}.$$
(22)
The free energy of the uncondensed solution is found by simply crossing out the condensate density in (21):

$$\stackrel{~}{A}_\mathrm{u}=\stackrel{~}{\mu }Nd^2x\left(\stackrel{~}{\gamma }\stackrel{~}{n}^2+\frac{t^2}{4\pi }g_2(e^{(x^2+2\stackrel{~}{\gamma }\stackrel{~}{n}\stackrel{~}{\mu })/t})\right).$$
(23)
Equations (21) and (23), of which more-general versions are derived in LABEL:HoubStoof by a different method HuseSiggia , can be easily seen to become identical at temperatures high enough that the condensate density can be neglected. On the other hand, their low-temperature limits differ, since Eq. (23) tends to a higher value than that attained by (21). This can be traced back to the fact that the grand-canonical ensemble has to be changed in order for it to correctly describe the particle-number fluctuations at low temperatures HoubStoof ; BergemLi . Now, a 2D Bose system has a condensate at least at $`T=0`$, so the two free energies should coincide there; however, we must keep in mind that the semiclassical approximation used to derive Eq. (23) requires that $`k_\text{B}T\mathrm{}\omega `$; thus the method we are using *cannot* describe the appearance of the zero-temperature condensate in the uncondensed solution. At the end of the next section we will study further consequences of this distinction.

## III Numerical Methods and Results

There is a time-honored prescription Griffin96 ; GPS  for finding the self-consistent solution of the Hartree-Fock equations in the presence of a condensate: Initially, we assume that only the condensate is present and solve the Gross-Pitaevskiĭ equation (5) for $`\stackrel{~}{n}^{}=0`$. The wavefunction and eigenvalue that result are fed into the Hartree-Fock expression for the density (13), which, when integrated over all space, yields also a value for the noncondensate fraction $`N^{}`$; one can then readjust the condensate fraction and solve the Gross-Pitaevskiĭ equation that results. The process is then iterated until the chemical potential and the particle fractions stop changing.

By far the most difficult part of this process is the solution of the nonlinear eigenvalue problem (5). Different methods exist in the literature; we have obtained identical results by solving it as an initial value problem edw  and, much more efficiently, by employing the method of spline minimization Krauth ; HKN . This method uses the fact that Eq. (5) is the Euler-Lagrange equation that minimizes the functional

$$J[\stackrel{~}{\mathrm{\Phi }}]=d^2x\left[(\stackrel{~}{\mathbf{}}\stackrel{~}{\mathrm{\Phi }})^2+\stackrel{~}{\mathrm{\Phi }}(x^2+2\stackrel{~}{\gamma }\stackrel{~}{n}(𝐱))\stackrel{~}{\mathrm{\Phi }}+\frac{\stackrel{~}{\gamma }}{2}\stackrel{~}{\mathrm{\Phi }}^4\right].$$
(24)
After setting a small, nonuniform grid of fixed abscissas that represent the coördinate $`x`$, we take the corresponding ordinates, which represent $`\stackrel{~}{\mathrm{\Phi }}`$, as the parameters to be varied until (24) attains its smallest possible value. Since we only have information about the value of the function at a discrete set of points, in order to calculate the necessary derivatives and to integrate we perform a cubic-spline interpolation; the integral is found using ten-point Gauss-Legendre quadrature. The minimization is carried out using the Nelder-Mead method.

We checked our code by comparing its predictions in three dimensions to previously published results. For the case studied in LABEL:HKN, the condensate fractions we found differed from those in the paper by less than one part in $`10^4`$. The code also reproduced previously known results for the ideal gas, including finite-size effects mulLong .

Figure 1, already discussed above, shows one of the solutions that we have found using the Hartree-Fock approximation. The gas has $`N=10^4`$ atoms; the coupling constant $`\stackrel{~}{\gamma }`$ has been chosen so that the system has approximately the same radius as the three-dimensional gas studied in LABEL:HKN, where parameters resembling those of the original JILA trap Wieman  are used. The system is shown at $`T=0.7T_c`$, where $`T_c`$ is the condensation temperature for the ideal gas. We have also found solutions for the Goldman-Silvera-Leggett model, which corresponds to the large-$`N`$ limit of the GP equation. This was done by treating the problem as a simultaneous system of nonlinear equations on a uniform grid and solving it with a least-squares method.

When $`N=10^4`$, as in Fig. 1, it is not possible to find a self-consistent solution for the condensed equations beyond $`T0.8T_c`$, since $`\stackrel{~}{ϵ}_{\text{HF}}`$ becomes negative; this had already been noted in LABEL:HKN for three dimensions, where it was interpreted as a finite-size effect, and occurs at even lower temperatures for 2D. The limitation becomes more severe as $`N`$ increases: for $`N=10^6`$ we cannot find solutions beyond $`T0.5T_c`$. The authors of LABEL:kim2 report predictions at temperatures very close to the transition by using a finite-size correction to the chemical potential nearTc . Our inability to work above certain temperatures might be a consequence of using the semiclassical approximation, though we point out that the Popov approximation is expected to break down close to the transition temperature GPS ; RMPBEC .

We also found self-consistent three-dimensional solutions using the more general Bogoliubov spectrum by applying the same method and using the Hartree-Fock solutions as a starting point for the iteration. We find that these solutions exhibit an enhanced depletion of the condensate, in agreement with those found by other authors nearTc . It was impossible, however, to find this kind of solution in two dimensions, even for systems with $`N`$ as low as $`100`$: after a few iterations, the chemical potential became too large, the Bogoliubov energies became imaginary, and the noncondensate density diverged. When the possibility of phononlike excitations is allowed, then, the condensate is destabilized, just as had been found in the Thomas-Fermi limit mulShort .

The uncondensed case, on the other hand, gives us solutions in these conditions, and to it we now turn. Equation (11) along with the particle-conservation condition (8) for all temperatures is most easily solved by rewriting (11) as

$$Ze^{x^2/t}=2e^{(\pi \stackrel{~}{\gamma })\nu (𝐱)}\mathrm{sinh}\pi \nu (𝐱)$$
(25)
where have introduced the fugacity $`Z=e^{\stackrel{~}{\mu }/t}`$ and $`\nu (𝐱)2\stackrel{~}{n}(𝐱)/t`$. Given $`Z`$, $`t`$, and $`\stackrel{~}{\gamma }`$ it is possible to find $`\nu `$ at every point using a standard root-finding algorithm. One then wants to find the value of $`Z`$ such that the total number of particles is $`N`$. It is better, however, to write Eq. (25) at the origin,

$$Z=2e^{(\pi \stackrel{~}{\gamma })\nu _0}\mathrm{sinh}\pi \nu _0,$$
(26)
eliminate $`Z`$ between (25) and (26), and solve for $`\nu _0=\nu (0)`$, the density at the center of the trap, using the same root finder.

In Fig. 2 we show the (Hartree-Fock) condensed and uncondensed solutions that we have obtained. They are similar in shape and exhibit identical behavior for large $`x`$. The uncondensed solution has a lower value at the origin and predicts a wider radial density profile.

We also calculated the free energy corresponding to each solution; the results are shown in Fig. 3, which shows $`\stackrel{~}{A}/N`$ as a function of temperature for both cases when $`N=10^3`$ and $`N=10^4`$. The free energies appear to coincide at high temperatures; this was to be expected, since Eqs. (21) and (23) become identical at temperatures high enough for the condensate density $`\stackrel{~}{n}_0`$ to be neglected. As for the low-temperature limit, we have already noted that the free energy of the condensed solution tends to the value predicted by the zero-temperature Thomas-Fermi limit, while that of the uncondensed solution tends to a higher value. This value can actually be calculated: it is easy to show bhaduri  that the low-$`T`$ limit of Eq. (25) is

$$2\stackrel{~}{\gamma }n(𝐱)=\stackrel{~}{\mu }x^2,$$
(27)
the Thomas-Fermi limit but with $`\stackrel{~}{\gamma }`$ replaced by $`2\stackrel{~}{\gamma }`$; a larger interaction strength, as we can see from (22), implies a higher free energy. Interestingly, when we compare a given Hartree-Fock solution to an uncondensed solution with half the interaction strength, we find that the density profiles become very similar (and are indistinguishable at $`T=0`$), while the free energies coincide almost exactly for a wide range of temperatures; attractive as this possibility might be, however, it has a serious flaw: the free energies start to differ as the temperature increases, where they should coincide by definition. This tells us that the factor of $`2`$, which turns out to be the same one discussed after Eqs. (5) and (23), has to be retained; the WKB approximate expression for the uncondensed state is valid *only* for $`k_\text{B}T\mathrm{}\omega `$, so dropping the factor of $`2`$ in order to match the $`T=0`$ condensate is invalid. (The authors of LABEL:bhaduri, in fact, omit this factor from their paper, although they address this question in a subsequent publication brandon ). It is in fact this factor that guarantees that the uncondensed solution has a higher free energy than the condensed one at all temperatures; this, despite our inability to find condensed solutions using the Bogoliubov energy spectrum, leads us to conclude that, at least at this level of approximation, the uncondensed solution is unphysical and the two-dimensional finite trapped system will exhibit some sort of condensation at finite temperature.

## IV Conclusion

We have found solutions to the two-dimensional HFB equations in the Hartree-Fock approximation, both for finite numbers of atoms and in the thermodynamic limit; still, when we try to go beyond this scheme and consider the whole Bogoliubov excitation spectrum, the low end of which is described by phonons, we are unable to find self-consistent solutions for finite—even low—values of $`N`$. We have seen that it is possible to describe the system as having no condensate at all, but these solutions, at least for the parameter combinations that we have studied, have a higher free energy than their partly-condensed counterparts; this leads us to conclude that the 2D system will have some kind of condensate at low enough temperatures, insofar as the Hartree-Fock approach can successfully describe the quasicondensate that we expect to find in a finite system.

One could argue that the whole mean-field approach we have adopted is wrong in two dimensions, since a condensation into a single state is being assumed from the start. However, an alternative treatment Shev92  that does not make this assumption ends up with equations identical to (6), so we are left none the wiser. We have also bypassed the fact that the interaction strength $`g`$ is not really constant in 2D, but rather depends logarithmically on the relative momentum Shev91 ; petrov ; this, however, has been found to be of little consequence brandon . Another possibility is that condensation occurs into a band of states, forming a “smeared” or generalized condensate Gir ; Noz ; this alternative situation is not ruled out by the standard proof of Hohenberg’s theorem MulHT .

Monte Carlo simulations *do* seem to predict the presence of a condensate in two dimensions Heinrichs ; Pearson . Furthermore, two-dimensional condensates appear to have been produced in the laboratory gorlitz ; jakkola . Presumably the MC simulations show both the effects of a quasicondensation of a finite system and the BKT transition to the superfluid state. Our attempt here has been to test the possibility of representing these Bose effects by a relatively simple set of HFB or Hartree-Fock equations. While no solutions can be found for the HFB set, even for a finite system, the Hartree-Fock equations do provide a description of a condensed state. We feel there is reason to believe that this description should be a fair representation of the actual situation.

###### Acknowledgements.

Markus Holzmann lent us the program he wrote for LABEL:HKN and thus allowed us to compare our results to those in the reference; we are grateful to him for that, and to P. G. Kevrekidis and E. Navayazdani for illuminating conversations. We also thank the Department of Astronomy at the University of Massachusetts for granting us access to their computational facilities. Finally, we acknowledge the input given to this paper by two anonymous referees.


format: text

text:
# Metallicity Evolution of Damped Lyman Alpha Systems In ΛCDM Cosmology

## 1 Introduction

Since their discovery, it has been recognized that studies of the damped Ly$`\alpha `$ systems help reveal the chemical enrichment history of the universe in neutral gas. These quasar absorption line systems, defined to exhibit $`N(\mathrm{HI})2\times \mathrm{\hspace{0.33em}10}^{20}\mathrm{cm}^2`$, dominate the neutral hydrogen content of the universe up to $`z4`$ (Storrie-Lombardi and Wolfe 2000; Peroux et al. 2002) and exhibit metal-line profiles which yield accurate metallicity measurements (e.g., Prochaska et al. 2001). The first comprehensive metallicity survey of the damped systems (Pettini et al. 1994) revealed a sub-solar metallicity in systems at $`z2`$ with a considerable scatter from system to system. By weighting these individual Zn metallicity measurements by the $`N(\mathrm{HI})`$ value of the system, the authors derived a global mean metallicity $`\mathrm{log}(Z_{HI}/Z_{})1.1`$ dex. This important quantity is independent of the morphology, age, and history of the individual damped Ly$`\alpha `$ systems. Further observations (Pettini et al. 1997) supported the initial work and suggested minimal evolution in $`Z_{HI}`$ from $`z12.5`$ (Pettini 1999).

Prochaska & Wolfe (2000a) extended metallicity surveys of the damped systems to $`z4.5`$ through echelle observations on the Keck I 10m telescope. Their results also favored minimal evolution in the metallicity of HI gas from $`z24`$, although their initial sample suffered from small number statistics. More recently, Prochaska, Gawiser, & Wolfe (2001) and Prochaska & Wolfe (2002) have enlarged the sample to $`N50`$ systems at $`z2`$ and have confirmed the main conclusions from Prochaska & Wolfe (2000a): (1) the $`N(\mathrm{HI})`$-weighted mean metallicity exhibits minimal evolution from $`z1.73.5`$ and a small but significant decline at $`z3.5`$; (2) no individual damped system shows a metallicity below 1/1000 solar; and (3) the unweighted mean metallicity closely tracks the weighted mean over this entire redshift range.

At low redshift $`(z1`$), the situation is far less certain. There are many fewer observations (Boissé et al. 1998; Pettini et al. 2000) and the majority of systems have very low $`N(\mathrm{HI})`$ values such that the weighted mean is entirely dominated by one or two DLAs. Although the current value for $`Z_{HI}`$ at $`z1`$ roughly matches the $`z2`$ values, it is very possible that this is the result of small number statistics (Kulkarni & Fall 2002).

On the theoretical front, various groups have composed models to predict the chemical evolution history of the universe. The majority of these (Pei & Fall 1995; Pei et al. 1999; Malaney & Chaboyer 1996; Edmunds & Phillips 1997) compute the evolution of gas metallicity globally. All such calculations have been constrained to yield $`Z_{HI}=Z_{}`$ at $`z=0`$ and all predict marked evolution in $`Z_{HI}`$ from $`z=4`$ to the present epoch. In comparison with the observations, these scenarios predict more substantial evolution in $`Z_{HI}`$ at high redshift than is currently observed. As our subsequent results would indicate, this global approach does not accurately reflect the unique evolution of the DLA sub population. Somerville, Primack, & Faber (2001) used semi-analytic models to derive metallicity evolution and find reasonable agreement with observations. It is, however, difficult to gauge whether the lack of treatment of the intergalactic medium in the semi-analytic models has substantially affected the metallicity in DLAs.

We have followed a direct approach using cosmological hydrodynamic simulations with galaxy formation. In an earlier paper (Cen & Ostriker 1999) we attempted to compute cosmic chemical evolution using a simulation with a baryonic mass resolution of $`1\times 10^8h^1\mathrm{M}_{}`$ and a spatial resolution of $`200h^1`$kpc comoving. While the agreement with observations found in a variety of environments ranging from rich clusters of galaxies to low column density Lyman alpha clouds was encouraging, we found ourselves unable to adequately address the evolution of metallicity of dense and small systems such as DLAs. Tissera et al. (2001) also used direct (SPH) simulations to study the chemical evolution of DLAs but their simulations suffered substantially from limited mass resolution. Since DLAa are the prime repository of cool gas phase baryons, we thought it important to study them using a new simulation especially designed for this purpose. Our new simulation has a baryonic mass resolution of $`3.5\times 10^5h^1\mathrm{M}_{}`$ and a spatial resolution of $`30h^1`$kpc comoving, and we now investigate the evolution of metallicity of DLAs. The central questions that we would like to address are: Why is the observed evolution of metallicity of DLAs so mild from $`z=40`$, while we know that the star formation primarily occurs during this redshift interval? What is the expected distribution of metallicity among DLAs? How does the evolution of metallicity of DLAs compare to that of the all gas in the universe globally averaged and to stellar metallicity? In what sense (if any) are DLAs progenitors of normal galaxies? A description of the simulation is given in §2. Results are presented defined in §3, followed by conclusions in §4.

We would like to point out at outset that we are interested in “cosmic chemical evolution” of DLAs, which is a convolution of true galactic evolution and observational selection of DLAs. This has a different meaning for investigators who work in galactic studies where the same term means internal galactic chemical evolution with some chosen boundary conditions. Our simulations are admittedly inadequate for explaining any detailed chemical evolution of individual galaxies, for which studies such as that by Mathlin et al. (2001) can make much finer predictions to compare with respective observations. But, as they pointed out, galaxy interaction and merger, and gas infall/outflow (which are known to be important processes in hierarchical structure formation) may have played a major role but were ignored in their treatment. Our simulations, on the other hand, take into account mergers, gas infall/outflow and other galaxy interactions and thus provide a complementary approach to this problem.

## 2 Simulations

The results reported on here are based on a new computation of the evolution of the intergalactic medium in a cold dark matter model with a cosmological constant. The relevant model parameters are: $`\mathrm{\Omega }_0=0.30`$, $`\mathrm{\Omega }_b=0.035`$, $`\mathrm{\Lambda }_0=0.70`$, $`\sigma _8=0.90`$, $`H_0=67`$km/s/Mpc and $`n=1.0`$. The cosmological model is normalized to both the cosmic microwave background (CMB) temperature fluctuations measured by COBE on large scales (e.g., Bunn & White 1997) and the observed abundance of clusters of galaxies in the local universe (e.g., Cen 1998), and it is close to both the concordance model of Ostriker & Steinhardt (1995) and consistent with the recent high redshift supernova results (Reiss et al. 1998) and CMB measurement on intermediate scales (de Bernardis et al. 2000; Balbi et al. 2000). Retrospectively, this model is close to WMAP (Wilkinson Microwave Anisotropy Probe) normalized model (Spergel et al. 2003). The simulation box size is $`L_{box}=25h^1`$Mpc having $`768^3`$ gas cells and $`384^3`$ dark matter particles with the mean baryonic mass in a cell being $`3.3\times 10^5h^1\mathrm{M}_{}`$ and the dark matter particle mass being $`2.0\times 10^7h^1\mathrm{M}_{}`$. The nominal spatial resolution for both gas and dark matter is comoving $`32h^1`$kpc. This mass resolution is considerably better than of most cosmological simulations, but the spatial resolution, while significantly inferior to that obtained in both the SPH and AMR schemes, is, we believe, adequate for the present purpose.

The initial conditions adopted are those for Gaussian processes with the phases of the different waves being random and uncorrelated. The initial condition is generated by the COSMICS software package kindly provided by E. Bertschinger (2001). The simulated box is an unconstrained and truly random realization of the model universe of that size. Numerical methods of the cosmological hydrodynamic code and input physical ingredients have been described in our earlier papers (Cen & Ostriker 1999a,b). Here we give a more detailed description.

The simulation integrates five sets of equations simultaneously: the Euler equations for gas dynamics, rate equations for different hydrogen and helium species at different ionizational states, the Vlasov equation for dynamics of collisionless particles, the Poisson’s equation for obtaining the gravitational potential field and the equation governing the evolution of the intergalactic ionizing radiation field, all in cosmological comoving coordinates. The gasdynamical equations are solved using the TVD shock capturing code (Ryu et al. 1993) on an uniform mesh. The rate equations are treated using sub-cycles within a hydrodynamic time step due to much shorter ionization timescales (i.e., the rate equations are very “stiff”). Dark matter particles are advanced in time using the standard particle-mesh (PM) scheme. The gravitational potential on an uniform mesh is solved using the Fast Fourier Transform (FFT) method.

The radiation field from $`1`$eV to $`100`$keV is followed in detail with allowance for self-consistently produced radiating sources and sinks in the simulation box and for cosmological effects, i.e., radiation transfer for the mean field $`J_\nu `$ is computed with stellar, quasar and bremsstrahlung sources and sinks due to Ly$`\alpha `$ clouds etc. In addition, a local optical depth approximation is adopted to crudely mimic the local shielding effects: each cell is flagged with an hydrogen “optical depth” equal to the product of neutral hydrogen density, hydrogen ionization cross section and the cell size; equivalent ones for neutral helium and singly-ionized helium are also computed. In computing the global sink terms for the radiation field the contribution of each cell is subject to the shielding due to its own “optical depth”. In addition, in computing the local ionization and cooling/heating balance for each cell the same shielding is taken into account to attenuate the external ioinizing radiation field.

The simulation computes for each cell and each timestep detailed cooling and heating processes due to all the principal line and continuum processes for a plasma of primordial composition. Metals ejected from star formation (see below) are followed in detail in a time-dependent, non-equilibrium fashion. Cooling due to metals is computed using a code based on the Raymond-Smith code assuming ionization equilibrium (Cen et al. 1995): at each time step, given the ionizing background radiation field, we compute a lookup table for metal cooling in the temperature-density plane for a gas with solar metallicity, then metal cooling rate for each gas cell is computed using the appropriate entry in that plane multiplied by its metallicity (in solar units). In addition, Compton cooling due to the microwave background radiation field and Comption cooling/heating due to the self-consistently produced (see below) X-ray and high energy background are also included.

We follow star formation using a well defined (heuristic but plausible) prescription used by us in our earlier work (Cen & Ostriker 1992,1993) and similar to that of other investigators (Katz, Hernquist, & Weinberg 1992; Katz, Weinberg, & Hernquist 1996; Steinmetz 1996; Gnedin & Ostriker 1997). A stellar particle of mass $`m_{}=c_{}m_{\mathrm{gas}}\mathrm{\Delta }t/t_{}`$ is created (the same amount is removed from the gas mass in the cell), if the gas in a cell at any time meets the following three conditions simultaneously: (i) flow contracting, (ii) cooling time less than dynamic time, and (iii) Jeans unstable, where $`\mathrm{\Delta }t`$ is the timestep, $`t_{}=\mathrm{max}(t_{\mathrm{dyn}},10^7`$yrs), $`t_{dyn}`$ is the dynamical time of the cell, $`m_{\mathrm{gas}}`$ is the baryonic gas mass in the cell and $`c_{}=0.25`$ is star formation efficiency. In essence, we follow the classic work of Eggen, Lynden-Bell & Sandage (1962) and assume that the dynamical free-fall and galaxy formation timescales are simply related. Each stellar particle has a number of other attributes at birth, including formation time $`t_i`$, initial gas metallicity and the free-fall time in the birth cell $`t_{dyn}`$. The typical mass of a stellar particle in the simulation is about one million solar masses; in other words, these stellar particles are like coeval globular clusters.

Stellar particles are subsequently treated dynamically as collisionless particles, except that feedback from star formation is allowed in three forms: UV ionizing field, supernova kinetic energy, and metal rich gas, all being proportional to the local star formation rate. Supernova energy and metals are ejected into the local gas cells where stellar particles are located. Supernova energy feedback into the intergalactic medium (IGM) is included with an efficiency (in terms of rest-mass energy of total formed stars) of $`e_{SN}=1\times 10^5`$; i.e., for $`m_{}`$ of stars formed $`e_{SN}m_{}c^2`$ amount of energy from supernovae is released to the IGM, where $`c`$ is the speed of light. Two types of ionizing radiation sources are used: one characteristic of star formation regions and the other characteristic of quasars, with efficiencies (i.e., the fraction of rest-mass energy converted into ionizing radiation) of $`e_{UV,}=3\times 10^6`$, and $`e_{UV,Q}=5\times 10^6`$, respectively. We adopt the emission spectrum of massive stars from Scalo (1986) and that of quasars from Edelson & Malkan (1986). In addition, hot, shocked regions (like clusters of galaxies) emit ionizing photons due to bremsstrahlung radiation. The temporal release of the all three feedback components at time $`t`$ has the same form: $`(dt/t_{dyn})[(tt_i)/t_{dyn}]\mathrm{exp}[(tt_i)/t_{dyn}]`$, where $`t_i`$ is the formation time of a stellar particle. The UV component is simply averaged over the box, since the light propagation time across our box is small compared to the timesteps.

It is useful to state clearly how we treat the chemical evolution in our simulations. We did not separately make any adjustments to fit to the observed distributions and evolution of metals, but assumed a specific efficiency of metal formation, an adopted “yield” (Arnett 1996), the percentage of stellar mass that is ejected back into IGM as metals, of $`0.02`$. This is the only adjustable parameter which significantly effects metallicity and, to zero-th order, all quoted metallicity scale with this assumed yield. Metals in the IGM (assuming the standard solar composition) are followed as a separate variable (analogous to the total gas density) with the same hydrocode. In essence, we compute in current simulations the metals enrichment process of the IGM due to type II supernovae. Furthermore, we use Zn to represent the total metals (i.e., O), i.e., we adopt $`[Zn/O]=0`$. While the origin of Zn may be unclear, our justification is purely empirical based on observations. The assumption of $`[Zn/O]=0`$ is consistent with the extant, albeit, small set of DLAs whose oxygen abundances have been detected (Prochaska & Wolfe 2002; Pettini et al. 2002); there is no evidence to the contrary. We convert observed metallicities of DLAs expressed in other species to Zn abundance, calibrated by observed empirical relations for DLAs.

The model reproduces the observed UV background as a function of redshift, the redshift distribution of star formation (“Madau Plot”; Nagamine, Cen & Ostriker 1999), and the galaxy luminosity function (Nagamine et al. 2001b). We find that the computed metallicity distributions over a wide range of environments, including clusters of galaxies, damped Lyman systems, Lyman alpha forest and stars, are in broad agreement with observations (Cen & Ostriker 1999b; Nagamine et al. 2001a,b; Cen et al. 2001), lending us confidence that the computed metal distribution in the intermediate regions under consideration (between Lyman alpha forest and clusters of galaxies) in this high resolution simulation may be a good approximation to the model universe.

We identify a DLA with a system with a line-of-sight (LOS) neutral hydrogen column density higher than $`2\times 10^{20}`$cm-2. This method is similar to what is employed to observationally classify DLAs. In addition, we can easily identify, in three-dimension, host galaxies that are responsible for DLAs. A host galaxy is associated with a DLA by comparing their spatial locations. Galaxies in the simulation box are identified using a scheme described in Nagamine et al. (2000), with each having a set of attributes including mean formation redshift and initial gas metallicity. Each galaxy consists of sub-units each with formation time, mass, metallicity and we compute the luminosity of each galaxy using a population synthesis code (e.g., Bruzual 2000). Each galaxy is also associated with the dark matter mass around it or dark matter halo and virial velocity of the whole system. We note that multiple DLAs may be associated with a single host galaxies, simply because different LOS through a galaxy may be observed as DLAs. The typical impact parameter is comparable to the cell size. Throughout the paper we refer to ”DLA galaxies” as the host galaxies of DLAs.

## 3 Results

### 3.1 Pictures

A visual impression of the gas distribution in our box is presented in Figure (1), which shows the total gas distribution at redshift $`z=3`$. DLAs are typically located in or near the highest density regions denoted in red. More quantitatively, in Figure (2) we show two slices of gas density distribution with a size of $`25\times 25\times 3.125h^3`$Mpc3 each at $`z=0`$. The left slice is typical and the right slice contains the richest cluster in the simulation. We see that rich clusters are hostile environments for DLAs; DLAs tend to reside in filaments which bridge and surround clusters. This is due to the same physical effect described in Blanton et al. (1999). DLAs are positively correlated with high total density but negatively correlated with high environmental temperature (at least at the high temperature end); the specification of one variable - density - does not enable one to predict the likely occurance of collapsed systems. While more quantitative comparisons can not be made, this finding appears to be consistent with the fact that the observed spiral galaxies in clusters have lower HI content (Haynes et al. 1990; Haynes & Giovanelli 1991). In Figure 2 red dots show galaxies having stellar mass $`10^{10.5}\mathrm{M}_{}`$ and we see that roughly 30% of DLAs today are associated with massive galaxies like the Milky Way. We will return to this issue to present more quantitative results later.

### 3.2 Quantitative Results

Let us now turn to quantitative results. Figure (3) shows the metallicity distribution of DLAs as a function of neutral hydrogen column density, $`N_{HI}`$. At each redshift there is a weak trend of higher metallicity correlating with higher column density. Although the variance is large, as indicated by the quartile errorbars, this is, at first sight, inconsistent with observations where the opposite trend is seen. There may be important observational selection effects that have not yet been included in our analyses, which have caused this apparent inconsistency. In particular, dust obscuration may have played a large role. We will explore this possibility here.

To illustrate the dust obscuration effect we show, in each panel in Figure 3, a dashed line with $`\mathrm{log}(N_{HI})+[Zn/H]=20.5`$, intending to approximately separate the region that is expected to be affected by dust obscuration (upper right) and the region that is not (lower left), partly motivated by observations. If DLAs in the upper right side of the dashed lines had obscured QSOs, one would expect to see a large deficit of DLAs there. In fact, there is some observational indication that this effect may be there (Boissé et al. 1998). In Figure 4 we plot the median metallicities (with quartile errorbars) for the “unobscured” DLA population at $`z=23`$. We see that the weak positive correlation between the median metallicity and column density of DLAs seen in Figure 3 is now reversed, with a result that is consistent with the observed trend (Boissé et al. 1998).

As a consistency check we show in Figure 5 the column density distributions of DLAs at three redshifts ($`z=0,2,4`$). We plotted the total DLA population (squares) as well as the “unobscured” DLA population, as defined in Figure 3. Quite surprisingly, the “unobscured” DLA population provides a remarkably good fit to observations shown as the thick solid curve (Storrie-Lombardi & Wolfe 2000). Therefore, our assumption of dust obscuration is at least self-consistent and the model would provide a reasonable fit to observations, if the observed DLAs are a subset of total DLA population that is not obscured by dust. There is some indication that the model may be producing too many DLA at the low column density end even with the dust obscuration effect taken into account. This, however, may be telling us something about the small-scale power in the adopted cosmological model. It may be explained if the adopted LCDM model has too much small scale power, resulting in more numerous dense neutral cores; note that the small-scale power makes a much larger contribution to the DLA population at high redshift. This implication would be consistent with other independent observations including dwarf galaxies (Flores & Primack 1994; Burkert 1995; Klypin et al. 1999; Moore et al. 1999; McGaugh & de Blok 1998; Navarro & Steinmetz 2000). However, it could also be due to limited numerical simulation of the present simulation, which artificially enlarges density concentrations and thus produces more DLAs than there should be. Future high resolution simulations can check this. However, given that the metallicity dependence on column density is very weak, the conclusions drawn with regard to the metallicity of DLAs and its evolution are relatively robust.

Visual examination of Figure 3 may give the reader the impression that there appears to exist a significant number of DLAs at $`N_{HI}>10^{21}`$cm-2 in the model even for the “unobscured” DLA sub population. A more quantitative comparison with observations is now made. The simulation and observations indicate, respectively, that about (30%, 27%) of all DLAs have $`N_{HI}>10^{21}`$cm-2 at $`z=2`$. At $`z=3`$ simulation produces about 26% of DLAs having $`N_{HI}>10^{21}`$cm-2 , while the corresponding observed number is 23%. The simulation is also in good agreement with observations at low redshift ($`z00.5`$) where a high rate of occurance of such high column density DLAs is seen (Rao & Turnshek 2000). These results are consistent with the close agreement seen in Figure 5. At $`z>3.5`$, however, simulation gives 20% of DLA with $`N_{HI}>10^{21}`$cm-2 at $`z=4`$ versus 2.5% for observed DLAs with $`N_{HI}>10^{21}`$cm-2. While the number of observed DLAs at the highest redshift range is still small, this discrepancy at high redshift may be revealing the same alleged flaw of excessive small-scale power in the simulated model.

The neutral gas contained in DLAs is another extensively studied quantity, for its obvious link to stars in that the latter must form out of the former, and is closely related to the column density distribution shown in Figure 5. (But note that the converse statement is not necessarily true; not all DLAs are transformed to galaxies. We find many fall into clusters of galaxies and are added to hot, shocked gas.) It thus seems warranted to present the results from the simulation. The reader is reminded that, unlike the metallicity distributions but like the column density distribution shown in Figure 5, the computed fraction of the neutral gas in DLAs is much more prone to numerical effects as well as variations in the cosmological model. Figure 6 presents the neutral gas in DLAs as a function of redshift, compared to extant observations. Note that both the computed and observed neutral gas in DLAs are subject to the uncertainties on the upper limit of the column density distribution, since it is dominated the highest column density DLAs (given that the slope of the column density distribution is substantially flatter than $`2.0`$). Consequently but unsurprisingly, the computed neutral gas in total DLAs is significantly (about a factor $`1.52.5`$) higher than the observed values at all redshifts, in agreement with the overabundance of high column density DLAs in the simulation seen in Figure 5. But the “unobscured” DLAs provide a reasonable fit to observations, given both observational errors (especially at low redshifts) and numerical uncertainties. Our assumption of dust obscuration is again self-consistent.

Clearly, the real situation due to dust obscuration is much more complicated (see Fall & Pei 1993 for a more detailed treatment). For example, it is not clear what the relationship between dust content and metallicity should be and how large the scatter is. In addition, dust properties are largely uncertain. Nevertheless, such an approximate treatment evidently does have a quite large effect. It appears that the computed model would be in fair agreement with observations with regard to several major diagnostics shown above, if dust obscuration has played an important role. It is, however, too early to celebrate the existence of dust obscuration. There are other observations that indicate dust obscuration effect may not be as large as considered here (Ellison et al. 2001; Prochaska & Wolfe 2002). It is worth noting that we find that in the redshift range $`z=23`$ the fraction of the number of excluded DLAs due to dust obscuration (those that lie above the dashed lines in Figure 3) is 9.9%. Therefore, in the CORALS radio survey of Ellison et al. (2001) one would have expected to see, on average, only $`1.8`$ DLAs, that are dust obscured, in their sample of total $`19`$ DLAs. Clearly, a sample size of several-fold larger than the current one would greatly firm up the statistical significance of their important findings and determine the importance of dust obscuration unambiguously. Although the main conclusion of this paper that the metallicity of DLAs is substantially sub-solar and its evolution is relatively mild does not require the adoption of dust obscuration as proposed, other properties of DLAs, including column density distribution, evolution of neutral gas content and metal content in DLAs (see Figure 17 below and related text), are all seen to become consistent with observations under that assumption. Alternatively, there are at least two other possibilities to potentially explain the apparent discrepancies. First, the adopted cosmological model may contain flaws and require revision with regard to its properties related to the formation of DLAs. However, the adopted model is known to reproduce many other observations on scales that overlap with that related to DLAs, including galaxy abundance and clustering, although there is independent evidence that the adopted model may possess excessive power on small scales (see below). Second, it is possible that astrophysical processes related to star formation have not been treated adequately. In particular, metal enrichment of the intergalactic medium treated here is clearly oversimplified. We think this is a possibility that requires substantially more work and deserves further investigation with future simulations that possess more refined treatments of galaxy formation.

Let us now return to the main issue of metallicity evolution of DLAs. To compute this, we adopt two averaging methods used in analyzing observational data sets. The first method weights each DLA equally, $`<[Zn/H]>_1\frac{_N[Zn/H]}{N}`$, and the second method weights each DLA by its neutral hydrogen column density, $`<[Zn/H]>_2\mathrm{log}_{10}\left(\frac{N_{Zn}}{{\scriptscriptstyle N_{HI}}}\right)`$. In Figure 7 we show $`<[Zn/H]>_1`$ (thin solid curve) and $`<[Zn/H]>_2`$ (thin dashed curve) as a function of redshift, averaged over all DLAs identified in the simulation. The thick solid and thick dashed dotted curves in Figure 7 are $`<[Zn/H]>_1`$ and $`<[Zn/H]>_2`$ averaged over only the “unobscured” DLAs. Also shown in Figure 7 as a thin dotted curve is the metallicity averaged over all the gas in the universe (using the second averaging method).

Several interesting points may be extracted from Figure 7. The evolution of $`<[Zn/H]>_1`$ of all DLAs (thin solid curve) follows rather closely the universal metallicity evolution of all gas (thin dotted curve), which is mainly a coincidence. The evolution of $`<[Zn/H]>_2`$ of all DLAs (thin dashed curve), while approximately parallel to $`<[Zn/H]>_1`$ and evolution of all gas, has an amplitude about 0.5 dex higher. This indicates that formation of dense neutral cores occurs preferentially in regions that have higher metallicity than the average regions of the universe at all epochs, with the trend that deviation from the mean increases with redshift. However, this systematic difference between $`<[Zn/H]>_1`$ and $`<[Zn/H]>_2`$ is substantially reduced to $`0.2`$ dex, if dust obscuration has a significant selection effect on DLAs, as evident from a comparison between the two thick curves, except at the high redshift end ($`z3`$).

There appears to be a reasonable agreement between the model and observations with regard to the evolution of metallicity of unobscured DLAs for both neutral hydrogen column density weighted and unweighted means. There is some mild indication that the unweighted mean of the model may be a little lower than observed or equivalently, there is an overabundance of low metallicity low column density DLAs in the simulation. A possible explanation may be that the current simulation underestimates star formation in small galaxies due to limited resolution; with an ideal resolution, metal enrichment from additional stars in small galaxies would have raised the metallicity in very low metallicity regions ($`[Zn/H]3`$) in the present simulation to a higher level. Given the cost of the simulation, we will have to investigate this issue in future simulations. In addition, our oversimplified treatment of chemical evolution might have contributed to the discrepancy.

The most striking property that is borne out from both the simulation and observations is that a typical DLA has a metallicity substantially below the solar value and does not show significant evolution with redshift. Both can be understood by considering selection bias and the source and sink terms for objects termed DLAs. If one were to follow, in a Lagrangian sense, a typical region of dense gas one would find it steadily merging with other such units and finding itself (on average) in systems of higher and higher observed neutral hydrogen column density. Also, since star formation, whether steady or episodic, leads to a monotonic increase in metallicity, the fluid element would show an ever higher metallicity if observed. At some point in time it would be in a system with neutral column density $`N_{HI}>10^{20}`$cm-2 and, when observed would be called a “DLA” system. The physical system in which it is embedded would evolve to higher $`[Zn/H]`$ over time and ultimately be identified with either an optical galaxy or, if it should find itself in a very high temperature high density region and the gas be stripped/evaporated, to become an E or S0 system. Typically by the time the metallicity had reached solar values either gas would have been depleted (turned into stars), expelled (by supernovae), or stripped due to interaction with a hot phase intergalactic medium and so the system would no longer have $`N_{HI}>10^{20}`$cm-2. All the processes are included in our code, but of course limited by our numerical resolution. Subsequent fluid elements passing along this same sequence would have roughly the same metallicity $`[Zn/H]`$ at the same values of column density $`N`$. Thus the $`[Zn/H](N_{HI})`$ relation would end up evolving far more slowly as a function of redshift than would the metallicity of a given fluid element.

To further elucidate this physical picture, we show, in Figure 8, a comparison between the metallicity of DLA gas and the metallicity of stars. In addition to the globally averaged metallicity of DLA gas we show the stellar metallicity of all stars and stars in DLA galaxies, and gas in all galaxies and in massive (with stellar mass $`10^{10.5}\mathrm{M}_{}`$) galaxies, respectively. A few interesting features may be noted from Figure 8. First, at all epochs the average stellar metallicity is always higher than the gas metallicity of DLAs or all galaxies, which in turn is higher than the average gas in the universe (see Figure 7). This finding is very interesting and perhaps seems counter-intuitive at first instance: since stars are formed out of gas in the galaxy, one naively would expect that the gas should always have a higher metallicity than stars due to enrichment by subsequently formed stars. This occurs in the “closed box” model, where the metallicity of gas in galaxies is typically twice the metallicity of the stars in that galaxy. But, in the hierarchical structure formation model, fresh, lower metallicity gas is continuously accreted onto DLAs (or DLAs-to-be) or DLAs form from mergers of smaller systems with gas of lower metallicity. In addition, metal rich gas is steadily ejected from galaxies through the action of supernovae. This process has been invoked to solve the so-called G-dwarf problem (Ostriker & Thuan 1975). It is also interesting to note that the gas in the massive galaxies at high redshift ($`z>2.5`$) does have higher metallicities than their resident stars, probably indicating that these massive galaxies constitute relatively closed boxes.

Second, the gas in the most massive galaxies has noticeably higher metallicity than both gas in DLAs and in all galaxies. This is consistent with the picture outlined above in that higher metallicity systems can not maintain the required high neutral hydrogen column density due to either hotter environment or hotter internal gas.

Third, stars in DLAs (open triangles) generally have higher metallicities than typical stars (solid triangles), especially at redshift $`z>1`$. This would be consistent with observations that less massive galaxies generally have lower metallicity than higher mass galaxies, if DLAs are more massive than a typical galaxy at high redshift. This indeed turns out to be the case \[see Figure 9 below, where we show that DLAs are more massive than the $`L_{Sch}^{}`$ (one of the two parameters in the Schechter function) galaxies at $`z>2`$.\] At lower redshift, the trend is somewhat reversed (also consistent with Figure 9.)

Fourth, the metallicity of the gas in DLAs (solid squares) displays an interesting crossover at $`z3`$: at $`z>3`$ the average gas in all galaxies (open squares) has a higher metallicity than the average gas in DLAs, while the reverse is true at $`z<3`$. This may indicate that at $`z<3`$ star formation in low mass systems is less efficient than in more massive systems so that gas is less metal enriched. This trend is further consistent with the result that the gas in the most massive galaxies (solid dots in Figure 6) has a metallicity higher than DLAs (solid squares).

Finally, we see by $`z=0`$ that the gas even in the most massive galaxies has a metallicity lower than average stars or stars in DLAs (stars in the same massive galaxies have comparable and very slightly higher metallicity than stars in DLAs at $`z=0`$; not plotted). If one associates the most massive galaxies with large spiral galaxies such as the Milky Way, one would expect to see gas metallicity typically higher than stellar metallicity. What is happening here is that these most massive galaxies are elliptical systems and live in clusters where galactic gas has been stripped off; thus the gas metallicity for these massive galaxies is actually the gas metallicity of the intra-cluster gas, whose value ($`1/3`$ solar) seems consistent with observed cluster gas metallicity (Arnaud et al. 1994; Mushotzky et al. 1996; Tamura et al. 1996; Mushotzky & Lowenstein 1997), indicating that the adopted yield based on the theory of stellar interior is approximately correct empirically. Note that we do not include metal enrichment due to type Ia supernovae in the present simulations. Just as a consistency check, we find that 115 galaxies having \[Zn/H\]=\[-0.1, 0.2\], which is quite consistent with the number of $`L^{}`$ galaxies that one expects from a simulation box of this size, indicating that Milky Way-like galaxies can properly be accounted for.

Needless to say galaxy formation and DLA formation are intimately related. It is important to know what type of galaxies constitute the host galaxies of DLAs. In Figure 9 we show the median stellar mass (with quartile errorbars) of DLAs as a function of redshift. Not surprisingly, and consistent with the depicted description above, DLAs correspond to progressively larger galaxies towards lower redshift. A typical DLA follows closely the $`M_{}`$ at all redshifts, with a slight tilt such that at low redshift ($`z<2`$) a typical DLA is sub $`M_{}`$ and the reverse is true at high redshift ($`z>2`$). While about 30% of DLAs at $`z=0`$ are due to $`L^{}`$ galaxies (consistent with the visual inspection of Figure 2), consistent with local observations (Rao & Turnshek 2000), DLAs at $`z=35`$ are typically in galaxies less luminous than $`0.1L^{}(z=0)`$ with some fractian in very faint galaxies ($`L<0.01L_{}`$). The median luminosity of a DLA, $`L_{DLA}(z)`$, in units of typical galaxy luminosity at that redshift, $`L^{}(z)`$, that is, $`(L_{DLA}/L^{})_z`$ decreases from $`1.1`$ to $`0.5`$ as redshift declines from $`z=3`$ to $`z=0`$, but the absolute luminosity of the median DLA system increases in the same interval by a factor of five from $`0.1L^{}(z=0)`$ to $`0.5L^{}(z=0)`$. We note that the stellar masses for a typical DLA galaxy at $`z=3`$ is $`M_{stellar}10^9\mathrm{M}_{}`$, about one order of magnitude lower than the stellar masses deduced for LBGs (Shapley et al. 2001). This implies that DLAs and LBGs are different populations or at least that a typical DLA galaxy does not correspond to an LBG (Colbert & Malkan 2001).

Figure 10 shows the cumulative number density of DLAs as a function of velocity width at $`z=0,1,2,3,4,5`$ (curves). The velocity width is defined as $`0.67v_{\mathrm{rot}}`$, where as $`v_{\mathrm{rot}}`$ is the rotation velocity of the DLA galaxy at radius $`100h^1`$kpc. The agreement between the simulation and observations is reasonably good at the high width end but the disagreement at the low width end is quite severe, where a factor of $`24`$ overabundance in the simulation is seen. This discrepancy was seen earlier in the column density distribution in Figure 5, which, as we suggested, may indicate excessive small-scale power in the adopted model. Put aside the difference for a moment, it is clear that DLAs are composed of galaxies with a broad range of velocity structures. Our simulation says that most of DLAs are due to small, relatively isolated galaxies. In the redshift range where comparisons can be made the results found here are in good agreement with extant observations (Steidel et al. 1995; Lanzetta et al. 1997; Le Brun et al. 1997; Miller, Knezek, & Bregman 1999; Bunker et al. 1999; Cohen 2000; Kulkarni et al. 2000). However, this finding is not at odds with local 21cm observations, which find most of the neutral hydrogen gas in the most massive galaxies. To clarify, panel (a) of Figure 11 shows the cumulative fraction of neutral hydrogen mass as a function of galactic stellar mass. Panel (b) of Figure 11 shows in a slightly different way the same information: the solid curve indicates the median galaxy stellar mass and the two dashed curves indicate the quartile points of galaxy stellar masses with regard to the fraction of total neutral hydrogen mass. Indeed, we find that most of the neutral hydrogen gas mass at $`z=0`$ is in large (spiral) galaxies; $`70\%`$ of neutral hydrogen gas mass at $`z=0`$ is in $`L_{}`$ or larger galaxies, consistent with observations (Zwaan et al. 2001).

What are the environments surrounding DLAs? Figure 12 shows the average “environmental temperature” in the vicinity of DLAs, defined as the average temperature over 124 cells surrounding the central cell of a DLA, as a function of redshift. Note that the actual (internal) temperature of the DLA gas is low (at $`<10^4`$K). We see that DLAs tend to reside in a higher temperature environment at lower redshift, close to $`L^{}`$ galaxies at $`z=0`$ but are in low density regions and filaments at high redshift. Clearly, it is unlikely that DLAs should be found in environments like groups or clusters of galaxies at any redshift, as indicated visually earlier in Figure 2 for $`z=0`$. We show in Figure 13 the probability distribution of the local mass overdensity (with a Gaussian smoothing window of radius $`1h^1`$Mpc) of all galaxies and DLA galaxies at three redshifts. Consistent with Figures (2,12), DLAs at $`z=0`$ avoid the highest density regions of clusters of galaxies (top panel of Figure 11) as well as very low density regions, while all galaxies extend to both high density regions of clusters of galaxies and very low density regions. At high redshift, DLAs still avoid very low density regions but not the high density regions, because the highest density regions do not correspond to rich clusters of galaxies. Figure 14 shows the proper peculiar velocity distribution of all galaxies and DLA galaxies at $`z=0,2,4`$. In full agreement with Figure 13, DLA galaxies tend to be found in relatively cold (low velocity) regions of filaments.

How old are DLA host galaxies? Figure 15 shows the absolute ages of DLAs (top panel) and the relative age (to the age of the Universe) (bottom panel) at three redshifts (z=0,2,4). We see that DLAs at $`z=0`$ are old with their ages close to the age of the universe and no DLAs are expected to be younger than $`5`$Gyr. Going to high redshifts DLAs become relatively younger and at $`z=4`$ the median age of DLAs is about half the age of the universe at that redshift. In addition, the spread in age distribution becomes larger at high redshifts with the full width spanning about one decade at $`z=24`$. We then subdivide the DLAs according to their metallicity and compute separately the age distributions of the subsets. Figure 16 shows the age distributions for three subsets at three different redshifts. We see that at $`z=0`$ a trend is clearly visible: high metallicity DLAs are older than their low high metallicity counterparts and no DLAs with $`[Zn/H]>1.0`$ are expected to be younger than $`10`$Gyr at $`z=0`$. However, at high redshift ($`z=24`$) there does not appear to exist any significant age segregation with respect to DLA metallicity.

The fact that a typical DLA galaxy, nearly independent of its metallicity, would be $`0.82`$Gyr old at $`z=24`$ is interesting in light of recent observations. Pettini et al. (2002) find that, based on a high resolution, high S/N sample of 15 DLAs, the nitrogen abundance in DLAs with $`[O/H]=2`$ to $`1`$ has not saturated by $`<z>=2.4`$. The indicated age of DLAs at the relevant redshift in our model, if taken to interpret the observational data, would imply a requirement that the release of nitrogen from intermediate mass stars be delayed by $`1`$Gyr, the typical age of DLAs derived here. If the delay is much shorter than the typical age of DLAs, then the nitrogen abundance in DLAs would have already saturated.

The ages of DLA galaxies appear to be comparable to those of LBGs (see Figure 10a of Shapley et al. 2001). It may suggest that LBGs have higher star formation rate than DLA galaxies, which, of course, is entirely consistent with the lower metallicity of DLAs.

Finally, Figure (17) shows the fractions of metal mass in gaseous phase as a function of redshift, divided into three components: DLAs (with temperature $`T<10^4`$K; solid curve), intracluster gas ($`T>10^7`$K; dashed curve), and the remaining IGM gas ($`10^4T10^7`$K; dotted curve). We see that at $`z1`$ roughly one half of the metals in the gas phase is in DLAs and the other half in the general intergalactic medium, while the fraction of the metals in very hot, rich clusters is small. But by $`z=0`$ only about 25% of metals are in DLAs, while the fraction in clusters has increased dramatically, reaching $`30\%`$, consistent with the fact that most clusters form at $`z<1`$ in this model. The decline of metal mass fraction in DLAs at low redshift is partly due to metal-enriched gas being collected in the intra-cluster gas and partly due to star formation out of metal rich gas at low redshift (as indicated by the long-dashed curve). Our calculations indicate that roughly $`40\%`$ of all metals are initially ejected into the IGM through supernovae and then re-formed into and locked up in stars by $`z=0`$; this fraction is lower at $`20\%`$ at $`z=35`$.

We note that most of the metals in DLAs are in the dust “obscured” high metallicity, high column density DLAs, according to our model (see Figures 5,6 and related discussion). Nevertheless, the fraction of total metals in the “unobscured” DLAs (thin solid curve in Figure 13) ranging from $`5\%`$ at $`z=0`$ to $`20\%`$ at $`z=5`$, in good agreement with observational estimates (Pettini 1999; Pagel 2001). If we identify the dust “obscured” DLA population with SCUBA galaxies at high redshift, the implication is that SCUBA galaxies may contain $`70\%`$ of metals residing in galaxies or $`40\%`$ of all metals at high redshift, which appears to be in agreement with recent observations (Dunne et al. 2002). In addition, since both the computed metallicity ($`\mathrm{Z}_{}/3`$) and abundance of X-ray clusters at low redshift are found to be in agreement with observations, the computed fraction of metals in the intracluster gas should be approximately correct.

Before concluding the paper, one word on numerical resolution is appropriate. Since the real DLAs have sizes comparable to our numerical resolution at $`z=0`$ (30kpc comoving; note that our numerical resolution is better at high redshift scaling as 30kpc/(1+z) kpc proper), it seems inevitable that our limited numerical resolution must have affected the results. In particular, the actual density and hence the column density of the computed DLAs may have been underestimated. The countervailing effect is that lower resolution may artificially enlarge the central dense regions thus produce more DLAs (if the central density is sufficiently high). Therefore, the net effect is not immediately obvious. In any case, to test the resolution effect we artificially lower the column density threshold that we use to identify DLAs from $`2\times 10^{20}`$cm-2 to $`1\times 10^{20}`$cm-2. We find that this change does not change the results presented in all the figures noticeably, except that in Figure 5 the unweighted metallicity (thin solid curve) of all DLAs is slightly lower by $`0.1`$ dex at $`z=3`$ and $`0.2`$ dex at $`z=4`$.

We should still stress that we do have enough resolution to resolve the internal structures of galaxies. All that we are trying to present in the paper are global properties averaged over all galaxies and intergalactic medium on scales comparable to or larger than our resolution, that do not sensitively depend on details on the sub-grid scales. To support our view that our treatment of star formation and metal production rates are reasonable and empirically in agreement with observations over a wide range of observables, we will refer to published papers based on the same simulation (Cen & Ostriker 1999; Nagamine et al. 2000; Nagamine et al. 2001a,b). To summarize, there are at least three lines of evidence that would indicate that the simulation is adequate for tracking global star formation, galaxy formation and metal production. First, we have computed the global star formation rate as a function of redshift (Nagamine etal 2000, 2001a;), i.e., the “Madau plot”, and found reasonable agreement with observations over the redshift range that can be compared, given the uncertainties concerning both observations and simulations. This indicates that the observed star formation history is reasonably well produced by the simulation of this currently best cosmological model. This would remove most of the uncertainties in the simulation concerning the amount of stars formed. Second, the computed luminosity function of galaxies are found to be in reasonable agreement with observations with a flat faint end slope of $`1.15`$ down to a galaxy mass of $`10^{8.5}\mathrm{M}_{}`$; Nagamine etal 2001b). Once again, this indicates that our simulation has sufficient resolution as to identify galaxies and compute the mass and luminosity of galaxies (but not to resolve internal structures of galaxies). Third, the resulting metallicity of the computed intra-cluster gas ($`1/3`$ solar value) is in excellent with observations, implying that our adopted metal yield value ($`0.02`$, based on stellar interior theory, see Arnett 2000) agrees with empirically determined observational value. Since the star formation history and galaxy luminosity function are independent of metal yield and most of the metals ending up in the intracluster medium are due to galaxies close to those hosting DLAs, the yield is essentially fixed, on which metallicity of DLAs mainly depend. In addition, our previous, lower resolution simulations (Cen & Ostriker 1999) have yielded results over a wide range of environments (from Lyman alpha forest to clusters of galaxies) that are in reasonable agreement with observations, further supporting the conclusion that our overall metal production rates as well as those regions previously resolved has been correctly computed.

## 4 Conclusions

We use a latest high mass resolution hydrodynamic simulation of a $`\mathrm{\Lambda }`$CDM model to compute the metallicity evolution of damped Lyman alpha systems. Contrary to a naive expectation, but consistent with observations, the computed true mean metallicity of all DLAs shows a very slow evolution with time. Observations appear to pick out damped Lyman alpha systems to be at an intermediate stage in galactic evolution. They tend to live in the filaments, where metal enriched gas from previous generations of small star forming subgalactic units has collected via the hierarchical merger process. As star formation and merging proceeds, they tend to turn into normal massive galaxies ultimately often residing in high density high temperature regions where the gas component is stripped off.

It is important to note that DLAs are not a simple population of galaxies but due to sightlines through a variety of different systems. The mixture of galaxy hosts of DLAs changes with redshift. The median luminosity of a DLA, $`L_{DLA}(z)`$, in units of typical galaxy luminosity at that redshift, $`L^{}(z)`$, that is, $`(L_{DLA}/L^{})_z`$ decreases from $`1.1`$ to $`0.5`$ as redshift declines from $`z=3`$ to $`z=0`$, but the absolute luminosity of the median DLA system increases in the same interval by a factor of $`5`$ from $`0.1L^{}(z=0)`$ to $`0.7L^{}(z=0)`$. We predict that DLAs should not be found in hot environments such as groups and clusters of galaxies, but tend to live in regions moderate overdensity with relatively cold flows.

The age of the DLA host galaxies approaches the age of the universe at $`z=0`$, while it is approximately half the age of the universe at $`z=34`$. At high redshift ($`z=24`$) there appears to exist no significant age segregation with respect to DLA metallicity.

We find that about 50% all metals in gaseous phase is in DLAs from $`z=5`$ to $`z=1`$, then decreasing rapidly to $`25\%`$ by $`z=0`$, as these metals are swept into hotter X-ray emitting cluster gas and recycled into newly formed stars. It is noted that about 50% of total metals are locked up in stars by $`z=0`$; this number is 20% at $`z3`$.

Observational selection effects (such as dust obscuration) may have significantly masked the true average metallicity of DLAs. Using only those DLAs that are located in the unobscured region of DLAs in the metallicity-column density plane (using a very simple dust obscuration model), we find that, while the trend of metallicity evolution remains little unchanged, the mean (column density-weighted) metallicity is reduced by 0.3 to 0.5 dex from $`z=0`$ to $`z=5`$, bringing the simulation results into good agreement with observations. The “unobscured” DLAs in the simulation also provide good matches to both the observed column density distribution, redshift evolution of the neutral gas content in DLAs and metal content in DLAs. Therefore, in order for the model to agree with observations for a range of variables either the assumed dust obscuration effect occurs in the real universe or the model has to be modified as to remove the high column density, high metallicity DLAs. The possibilities for the latter will be explored in the future. It should be emphasized, however, that current observations do not seem to require significant amount of dust in DLAs (Ellison et al. 2001; Prochaska & Wolfe 2002). We point out, however, that only $`1.8`$ DLAs are expected to be dust obscured (in optical) in the CORALS radio survey of Ellison et al. (2001) of total $`19`$ DLAs. Clearly, a sample several times larger is required to greatly firm up the statistical significance of their important findings and determines the importance of dust obscuration unambiguously. While our simulations strongly suggest the significance of dust obscuration with regard to DLA evolution, we will be forced to make a substantial revision of the computed model, if significant dust obscuration is firmly excluded observationally in the future. We stress, however, that some of the major conclusions of the present paper, including the mild metallicity evolution of DLAs, the large disparity of DLA host galaxies and the monotonic evolution of the increasing mass of DLA systems with time, are much less prone to these uncertainties and not expected to alter significantly.

Finally, we pointed out that that the adopted model predicts (approximately a factor of $`23`$) too many damped Lyman alpha systems, indicating either imperfections of the simulation, or the possibility that this $`\mathrm{\Lambda }`$CDM model has too much small scale power.

The work is supported in part by grants AST93-18185 and ASC97-40300. We thank K. Nagamine for kindly providing the catalog of simulated galaxies, Greg Bryan for kindly making a software program available to us which is used to produce Figure 1, and Max Pettini and Mike Fall for very useful comments.