text:
# Density Matrix Renormalization: A Review of the Method and its Applications

## 1 Introduction

The basics of the Density Matrix Renormalization Group were developed by S. White in 1992 and since then DMRG has proved to be a very powerful method for low dimensional interacting systems. Its remarkable accuracy can be seen for example in the spin-1 Heisenberg chain: for a system of hundreds of sites a precision of $`10^{10}`$ for the ground state energy can be achieved. Since then it has been applied to a great variety of systems and problems including, among others, spin chains and ladders, fermionic and bosonic systems, disordered models, impurities and molecules and 2D electrons in high magnetic fields. It has also been improved substantially in several directions like two (and three) dimensional (2D) classical systems, stochastic models, the presence of phonons, quantum chemistry, field theory, the inclusion of temperature and the calculation of dynamical and time-dependent properties. Some calculations have also been performed in 2D quantum systems. All these topics are treated in detail and in a pedagogical way in the book , where the reader can find an extensive review on DMRG. In this article we will attempt to cover the different areas where it has been applied without entering into details but in a few cases, where we have chosen some representative contributions. We suggest the interested reader to look for further information in the referenced work. Our aim here is to give the reader a general overview on the subject.

One of the most important limitations of numerical calculations in finite systems is the great amount of states that have to be considered and its exponential growth with system size. Several methods have been introduced in order to reduce the size of the Hilbert space to be able to reach larger systems, such as Monte Carlo, renormalization group (RG) and DMRG. Each method considers a particular criterion of keeping the relevant information.

The DMRG was originally developed to overcome the problems that arise in interacting systems in 1D when standard RG procedures were applied. Consider a block B (a block is a collection of sites) where the Hamiltonian $`H_B`$ and end-operators are defined. These traditional methods consist in putting together two or more blocks (e.g. B-B’, which we will call the superblock), connected using end-operators, in a basis that is a direct product of the basis of each block, forming $`H_{BB^{}}`$. This Hamiltonian is then diagonalized, the superblock is replaced by a new effective block $`B_{new}`$ formed by a certain number $`m`$ of lowest-lying eigenstates of $`H_{BB^{}}`$ and the iteration is continued (see Ref. ). Although it has been used successfully in certain cases, this procedure, or similar versions of it, has been applied to several interacting systems with poor performance. For example, it has been applied to the 1D Hubbard model keeping $`m1000`$ states. For 16 sites, an error of 5-10% was obtained . Other results were also discouraging. A better performance was obtained  by adding a single site at a time rather than doubling the block size. However, there is one case where a similar version of this method applies very well: the Kondo (and Anderson) model. Wilson mapped the one-impurity problem onto a one-dimensional lattice with exponentially descreasing hoppings. The difference with the method explained above is that in this case, one site (equivalent to an “onion shell”) is added at each step and, due to the exponential decrease of the hopping, very accurate results can be obtained.

Returning to the problem of putting several blocks together, the main source of error comes from the election of eigenstates of $`H_{BB^{}}`$ as representative states of a superblock. Since $`H_{BB^{}}`$ has no connection to the rest of the lattice, its eigenstates may have unwanted features (like nodes) at the ends of the block and this can’t be improved by increasing the number of states kept. Based on this consideration, Noack and White tried including different boundary conditions and boundary strengths. This turned out to work well for single particle and Anderson localization problems but, however, it did not improve the results significantly for interacting systems. These considerations led to the idea of taking a larger superblock that includes the blocks $`BB^{}`$, diagonalize the Hamiltonian in this large superblock and then somehow project the most favorable states onto $`BB^{}`$. Then $`BB^{}`$ is replaced by $`B_{new}`$. In this way, awkward features in the boundary would vanish and a better representation of the states in the infinite system would be achieved. White proposed the density matrix as the optimal way of projecting the best states onto part of the system and this will be discussed in the next section. Some considerations concerning the effect of boundary conditions in the nature of the states kept (and an analogy to the physics of black holes) is given in . The justification of using the density matrix is given in detail in Ref.. A very easy and pedagogical way of understanding the basic functioning of DMRG is applying it to the calculation of simple quantum problems like one particle in a tight binding chain .

In the following Section we will briefly describe the standard method; in Sect. 3 we will mention some of the most important applications; in Sect. 4 we review the most relevant extensions to the method and finally in Sect. 5 we concentrate on the way dynamical calculations can be performed within DMRG.

## 2 The Method

The DMRG allows for a systematic truncation of the Hilbert space by keeping the most probable states describing a wave function (e.g. the ground state) instead of the lowest energy states usually kept in previous real space renormalization techniques.

The basic idea consists in starting from a small system (e.g with $`N`$ sites) and then gradually increase its size (to $`N+2`$, $`N+4`$,…) until the desired length is reached. Let us call the collection of $`N`$ sites the universe and divide it into two parts: the system and the environment (see Fig. 1). The Hamiltonian is constructed in the universe and its ground state $`|\psi _0>`$ is obtained. This is considered as the state of the universe and called the target state. It has components on the system and the environment. We want to obtain the most relevant states of the system, i.e., the states of the system that have largest weight in $`|\psi _0`$. To obtain this, the environment is considered as a statistical bath and the density matrix is used to obtain the desired information on the system. So instead of keeping eigenstates of the Hamiltonian in the block (system), we keep eigenstates of the density matrix. We will be more explicit below.

Let’s define block \[B\] as a finite chain with $`l`$ sites, having an associated Hilbert space with, $`m`$ states where operators are defined (in particular the Hamiltonian in this finite chain, $`H_B`$ and the operators at the ends of the block, useful for linking it to other chains or added sites). Except for the first iteration, the basis in this block isn’t explicitly known due to previous basis rotations and reductions. The operators in this basis are matrices and the basis states are characterized by quantum numbers (like $`S^z`$, charge or number of particles, etc). We also define an added block or site as \[a\] having $`n`$ states. A general iteration of the method is described below:

i) Define the Hamiltonian $`H_{BB^{}}`$ for the superblock (the universe) formed by putting together two blocks \[B\] and \[B’\] and two added sites \[a\] and \[a’\] in this way: \[B a a’ B’ \] (the primes are only to indicate additional blocks, but the primed blocks have the same structure as the non-primed ones; this can vary, see the finite-size algorithm below). In general, blocks \[B\] and \[B’\] come from the previous iteration. The total Hilbert space of this superblock is the direct product of the individual spaces corresponding to each block and the added sites. In practice a quantum number of the superblock can be fixed (in a spin chain for example one can look at the total $`S^z=0`$ subspace), so the total number of states in the superblock is much smaller than $`(mn)^2`$. In some cases, as the quantum number of the superblock consists of the sum of the quantum numbers of the individual blocks, each block must contain several subspaces (several values of $`S^z`$ for example). Here periodic boundary conditions can be attached to the ends and a different block layout should be considered (e.g. \[B a B’ a’ \]) to avoid connecting blocks \[B\] and \[B’\] which takes longer to converge. The boundary conditions are between \[a’\] and \[B\]. For closed chains the performance is poorer than for open boundary conditions .

ii) Diagonalize the Hamiltonian $`H_{BB^{}}`$ to obtain the ground state $`|\psi _0`$ (target state) using Lanczos or Davidson algorithms. Other states could also be kept, such as the first excited ones: they are all called target states.

iii) Construct the density matrix:

$$\rho _{ii^{}}=\underset{j}{}\psi _{0,ij}\psi _{0,i^{}j}$$
(1)
on block \[B a\], where $`\psi _{0,ij}=ij|\psi _0`$, the states $`|i`$ belonging to the Hilbert space of the block \[B a\] and the states $`|j`$ to the block \[B’ a’\]. The density matrix considers the part \[B a\] as a system and \[B’ a’\], as a statistical bath. The eigenstates of $`\rho `$ with the highest eigenvalues correspond to the most probable states (or equivalently the states with highest weight) of block \[B a\] in the ground state of the whole superblock. These states are kept up to a certain cutoff, keeping a total of $`m`$ states per block. The density matrix eigenvalues sum up to unity and the truncation error, defined as the sum of the density matrix eigenvalues corresponding to discarded eigenvectors, gives a qualitative indication of the accuracy of the calculation.

iv) With these $`m`$ states a rectangular matrix $`O`$ is formed and it is used to change basis and reduce all operators defined in \[B a\]. This block \[B a\] is then renamed as block \[Bnew\] or simply \[B\] (for example, the Hamiltonian in block \[B a\], $`H_{Ba}`$, is transformed into $`H_B`$ as $`H_B=O^{}H_{Ba}O`$).

v) A new block \[a\] is added (one site in our case) and the new superblock \[B a a’ B’\] is formed as the direct product of the states of all the blocks.

vi) This iteration continues until the desired length is achieved. At each step the length is $`N=2l+2`$ (if \[a\] consists of one site).

When more than one target state is used, i.e more than one state is wished to be well described, the density matrix is defined as:

$$\rho _{ii^{}}=\underset{l}{}p_l\underset{j}{}\varphi _{l,ij}\varphi _{l,i^{}j}$$
(2)
where $`p_l`$ defines the probability of finding the system in the target state $`|\varphi _l`$ (not necessarily eigenstates of the Hamiltonian).

The method described above is usually called the infinite-system algorithm since the system size increases at each iteration. There is a way to increase precision at each length $`N`$ called the finite-system algorithm. It consists of fixing the lattice size and zipping a couple of times until convergence is reached. In this case and for the block configuration \[B a a’ B’ \], $`N=l+1+1+l^{}`$ where $`l`$ and $`l^{}`$ are the number of sites in $`B`$ and $`B^{}`$ respectively. In this step the density matrix is used to project onto the left $`l+1`$ sites. In order to keep $`N`$ fixed, in the next block configuration, the right block $`B^{}`$ should be defined in $`l1`$ sites such that $`N=(l+1)+1+1+(l1)^{}`$. The operators in this smaller block should be kept from previous iterations (in some cases from the iterations for the system size with $`N2`$).

The calculation of static properties like correlation functions is easily done by keeping the operators in question at each step and performing the corresponding basis change and reduction, in a similar manner as done with the Hamiltonian in each block. The energy and measurements are calculated in the superblock. In Ref. an interpretation of the correlation functions of systems at criticality is given in terms of wave function entanglement, conjecturing a modification of DMRG for these cases that preserves the entanglement.

A faster convergence of Lanczos or Davidson algorithm is achieved by choosing a good trial vector. An interesting analysis on DMRG accuracy is done in Ref. . Fixed points of the DMRG and their relation to matrix product wave functions were studied in  and an analytic formulation combining the block renormalization group with variational and Fokker-Planck methods in . The connection of the method with quantum groups and conformal field theory is treated in . There are also interesting connections between the density matrix spectra and integrable models via corner transfer matrices. These articles give a deep insight into the essence of the DMRG method.

## 3 Applications

Since its development, the number of papers using DMRG has grown enormously and other improvements to the method have been performed. We would like to mention some applications where this method has proved to be useful. Other applications related to further developments of the DMRG will be mentioned in Sect. 4.

A very impressive result with unprecedented accuracy was obtained by White and Huse  when calculating the spin gap in a $`S=1`$ Heisenberg chain obtaining $`\mathrm{\Delta }=0.41050J`$. They also calculated very accurate spin correlation functions and excitation energies for one and several magnon states and performed a very detailed analysis of the excitations for different momenta. They obtained a spin correlation length of 6.03 lattice spacings. Simultaneously Sørensen and Affleck also calculated the structure factor and spin gap for this system up to length 100 with very high accuracy, comparing their results with the nonlinear $`\sigma `$ model. In a subsequent paper they applied the DMRG to the anisotropic $`S=1`$ chain, obtaining the values for the Haldane gap. They also performed a detailed study of the $`S=1/2`$ end excitations in an open chain. Thermodynamical properties in open $`S=1`$ chains such as specific heat, electron paramagnetic resonance (EPR) and magnetic susceptibility calculated using DMRG gave an excellent fit to experimental data, confirming the existence of free spins 1/2 at the boundaries. A related problem, i.e. the effect of non-magnetic impurities in spin systems (dimerized, ladders and 2D) was studied in . In addition, the study of magnon interactions and magnetization of $`S=1`$ chains was done in , supersymmetric spin chains modelling plateau transitions in the integer quantum Hall effect in  and ESR studies in these systems was considered in . For larger integer spins there have also been some studies. Nishiyama and coworkers calculated the low energy spectrum and correlation functions of the $`S=2`$ antiferromagnetic Heisenberg open chain. They found $`S=1`$ end excitations (in agreement with the Valence Bond Theory). Edge excitations for other values of $`S`$ have been studied in Ref. . Almost at the same time Schollwöck and Jolicoeur calculated the spin gap in the same system, up to 350 sites, ($`\mathrm{\Delta }=0.085J`$), correlation functions that showed topological order and a spin correlation length of 49 lattice spacings. More recent accurate studies of $`S=2`$ chains are found in  and of $`S=1`$ chains in staggered magnetic fields  including a detalied comparison to the non-linear sigma model in . In Ref.  the dispersion of the single magnon band and other properties of the $`S=2`$ antiferromagnetic Heisenberg chains were calculated.

Concerning $`S=1/2`$ systems, DMRG has been crucial for obtaining the logarithmic corrections to the $`1/r`$ dependence of the spin-spin correlation functions in the isotropic Heisenberg model . For this, very accurate values for the energy and correlation functions were needed. For $`N=100`$ sites an error of $`10^5`$ was achieved keeping $`m=150`$ states per block, comparing with the exact finite-size Bethe Ansatz results. For this model it was found that the data for the correlation function has a very accurate scaling behaviour and advantage of this was taken to obtain the logarithmic corrections in the thermodynamic limit. Other calculations of the spin correlations have been performed for the isotropic  and anisotropic cases . Luttinger liquid behaviour with magnetic fields have been studied in , field-induced gaps in , anisotropic systems in  and the Heisenberg model with a weak link in . An analysis of quantum critical points and critical behaviour in spin chains by combining DMRG with finite-size scaling was done in .

Similar calculations have been performed for the $`S=3/2`$ Heisenberg chain . In this case a stronger logarithmic correction to the spin correlation function was found. For this model there was interest in obtaining the central charge $`c`$ to elucidate whether this model corresponds to the same universality class as the $`S=1/2`$ case, where the central charge can be obtained from the finite-size scaling of the energy. Although there have been previous attempts, these calculations presented difficulties since they involved also a term $`1/\mathrm{ln}^3N`$. With the DMRG the value $`c=1`$ was clearly obtained.

In Ref. , DMRG was applied to an effective spin Hamiltonian obtained from an SU(4) spin-orbit critical state in 1D. Other applications were done to enlarged symmetry cases with SU(4) symmetry in order to study coherence in arrays of quantum dots, to obtain the phase diagram for 1D spin orbital models and dynamical properties in a magnetic field.

Dimerization and frustration have been considered in Refs.  and alternating spin chains in .

The case of several coupled spin chains (ladder models) have been investigated in , spin ladders with cyclic four-spin exchanges in  and Kagome antiferromagnets in . Zigzag spin chains have been considered in  and spin chains of coupled triangles in . As the DMRG’s performance is optimal in open systems, an interesting analysis of the boundary effect on correlation functions is done in . Magnetization properties and plateaus for quantum spin ladder systems have also been studied. An interesting review on the applications to some exact and analytical techniques for quantum magnetism in low dimension, including DMRG, is presented in .

There has been a great amount of applications to fermionic systems such as 1D Hubbard and t-J models , Luttinger liquids with boundaries, the Falicov-Kimball model , the quasiperiodic Aubry-Andre chain and Fibonacci-Hubbard models. It has also been applied to field theory. The method has been very successful for several band Hubbard models, Hubbard ladders  and t-J ladders. Also several coupled chains at different dopings have been considered  as well as flux phases in these systems . Time reversal symmetry-broken fermionic ladders have been studied in  and power laws in spinless-fermion ladders in . Long-range Coulomb interactions in the 1D electron gas and the formation of a Wigner crystal was studied in . Several phases including the Wigner crystal, incompressible and compressible liquid states, stripe and pairing phases, have been found using DMRG for 2D electrons in high magnetic fields considering different Landau levels. Persistent currents in mesoscopic systems have been considered in .

Quite large quasi-2D systems can be reached, for example in  where a 4x20 lattice was considered to study ferromagnetism in the infinite-U Hubbard model; the ground state of a 4-leg t-J ladder in ; the one and two hole ground state in 9x9 and 10x7 t-J lattices in ; a doped 3-leg t-J ladder in ; the study of striped phases in ; domain walls in 19x8 t-J systems in ; the 2D t-J model in  and the magnetic polaron in a 9x9 t-J lattice in . Also big $`CaV_4O_9`$ spin-1/2 lattices reaching 24x11 sites have been studied. There have been some recent attempts to implement DMRG in two and higher dimensions  but the performance is still poorer than in 1D. A recent extension using a two-step DMRG algorithm for highly anisotropic spin systems has shown promising results.

Impurity problems have been studied for example in one-  and two-impurity  Kondo systems, in spin chains  and in Luttinger Liquids . There have also been applications to Kondo and Anderson lattices , Kondo lattices with localized $`f^2`$ configurations , the two-channel Kondo lattice on a ladder , a t-J chain coupled to localized Kondo spins and ferromagnetic Kondo models for manganites .

## 4 Other extensions to DMRG

There have been several extensions to DMRG like the inclusion of symmetries to the method such as spin and parity. Total spin conservation and continuous symmetries have been treated in  and in interaction-round a face Hamiltonians, a formulation that can be applied to rotational-invariant sytems like $`S=1`$ and 2 chains. A momentum representation of this technique  that allows for a diagonalization in a fixed momentum subspace has been developed as well as applications in dimension higher than one and Bethe lattices. The inclusion of symmetries is essential to the method since it allows to consider a smaller number of states, enhance precision and obtain eigenstates with definite quantum numbers. Other recent applications have been in nuclear shell model calculations where a two level pairing model has been considered and in the study of ultrasmall superconducting grains, in this case, using the particle (hole) states around the Fermi level as the system (environment) block.

A very interesting and successful application is a recent work in High Energy Physics. Here the DMRG is used in an asymptotically free model with bound states, a toy model for quantum chromodynamics, namely the two dimensional delta-function potential. For this case an algorithm similar to the momentum space DMRG was used where the block and environment consist of low and high energy states respectively. The results obtained here are much more accurate than with the similarity renormalization group and a generalization to field-theoretical models is proposed based on the discreet light-cone quantization in momentum space. Below we briefly mention other important extensions, leaving the calculation of dynamical properties for the next Section.

### 4.1 Classical systems

The DMRG has been very successfully extended to study classical systems. For a detailed description we refer the reader to Ref. . Since 1D quantum systems are related to 2D classical systems, it is natural to adapt DMRG to the classical 2D case. This method is based on the renormalization group transformation for the transfer matrix $`T`$ (TMRG). It is a variational method that maximizes the partition function using a limited number of degrees of freedom, where the variational state is written as a product of local matrices. For 2D classical systems, this algorithm is superior to the classical Monte Carlo method in accuracy, speed and in the possibility of treating much larger systems. A recent improvement of this method considering periodic boundary conditions is given in  and a detailed comparison between symmetric and asymmetric targetting is done in . TMRG has also been successfully used to renormalize stochastic transfer matrices in a study of cellular automatons. The calculation of thermodynamical properties of 3D classical statistical systems has been proposed where the eigenstate of the transfer matrix with maximum eigenvalue is represented by the product of local tensors optimized using DMRG.

A further improvement to this method is based on the corner transfer matrix, the CTMRG and can be generalized to any dimension.

It was first applied to the Ising model and also to the Potts model, where very accurate density profiles and critical indices were calculated. Further applications have included non-hermitian problems in equilibrium and non-equilibrium physics. In the first case, transfer matrices may be non-hermitian and several situations have been considered: a model for the Quantum Hall effect, the $`q`$-symmetric Heisenberg chain related to the conformal series of critical models and the anisotropic triangular nearest and next-nearest neighbour Ising models. In the second case, the adaptation of the DMRG to non-equilibrium physics like the asymmetric exclusion problem and reaction-diffusion problems  has shown to be very successful. It has also been applied to stochastic lattice models like in  and to the 2D XY model .

### 4.2 Finite-temperature DMRG

The adaptation of the DMRG method for classical systems paved the way for the study of 1D quantum systems at non zero temperature, by using the Trotter-Suzuki method . In this case the system is infinite and the finiteness is in the level of the Trotter approximation. Standard DMRG usually produces its best results for the ground state energy and less accurate results for higher excitations. A different situation occurs here: the lower the temperature, the less accurate the results.

Very nice results have been obtained for the dimerized, $`S=1/2`$, XY model, where the specific heat was calculated involving an extremely small basis set ($`m=16`$), the agreement with the exact solution being much better in the case where the system has a substantial gap. It has also been used to calculate thermodynamical properties of the anisotropic $`S=1/2`$ Heisenberg model, with relative errors for the spin susceptibility of less than $`10^3`$ down to temperatures of the order of $`0.01J`$ keeping $`m=80`$ states. A complete study of thermodynamical properties like magnetization, susceptibility, specific heat and temperature dependent correlation functions for the $`S=1/2`$ and 3/2 Heisenberg models was done in . Other applications have been the calculation of the temperature dependence of the charge and spin gap in the Kondo insulator, the calculation of thermodynamical properties of ferrimagnetic chains and spin ladders, the study of impurity properties in spin chains, frustrated quantum spin chains, t-J and spin ladders and dimerized frustrated Heisenberg chains.

An alternative way of incorporating temperature into the DMRG procedure was developed by Moukouri and Caron. They considered the standard DMRG taking into account several low-lying target states (see Eq. 2) to construct the density matrix, weighted with the Boltzmann factor ($`\beta `$ is the inverse temperature):

$$\rho _{ii^{}}=\underset{l}{}e^{\beta E_l}\underset{j}{}\varphi _{l,ij}\varphi _{l,i^{}j}$$
(3)
With this method they performed reliable calculations of the magnetic susceptibility of quantum spin chains with $`S=1/2`$ and $`3/2`$, showing excellent agreement with Bethe Ansatz exact results. They also calculated low temperature thermodynamical properties of the 1D Kondo Lattice Model and of organic conductors . Zhang et al. applied the same method in the study of a magnetic impurity embedded in a quantum spin chain.

### 4.3 Phonons, bosons and disorder

A significant limitation to the DMRG method is that it requires a finite basis and calculations in problems with infinite degrees of freedom per site require a large truncation of the basis states. However, Jeckelmann and White developed a way of including phonons in DMRG calculations by transforming each boson site into several artificial interacting two-state pseudo-sites and then applying DMRG to this interacting system (called the “pseudo-site system”). The idea is based on the fact that DMRG is much better able to handle several few-states sites than few many-state sites. The key idea is to substitute each boson site with $`2^N`$ states into $`N`$ pseudo-sites with 2 states. They applied this method to the Holstein model for several hundred sites (keeping more than a hundred states per phonon mode) obtaining negligible error. In addition, up to date, this method is the most accurate one to determine the ground state energy of the polaron problem (Holstein model with a single electron).

An alternative method (the “Optimal phonon basis”) is a procedure for generating a controlled truncation of a large Hilbert space, which allows the use of a very small optimal basis without significant loss of accuracy. The system here consists of only one site and the environment has several sites, both having electronic and phononic degrees of freedom. The density matrix is used to trace out the degrees of freedom of the environment and extract the most relevant states of the site in question. In following steps, more bare phonons are included to the optimal basis obtained in this way. This method was successfully applied to study the interactions induced by quantum fluctuations in quantum strings, as an application to cuprate stripes and the dissipative two-state system . A variant of this scheme is the “four block method”, as described in . They obtain very accurately the Luttinger liquid-CDW insulator transition in the 1D Holstein model for spinless fermions.

The method has also been applied to pure bosonic systems such as the disordered bosonic Hubbard model, where gaps, correlation functions and superfluid density are obtained. The phase diagram for the non-disordered Bose-Hubbard model, showing a reentrance of the superfluid phase into the insulating phase was calculated in Ref. . It has also been used to study a chain of oscillators with optical phonon spectrum  and optical phonons in the quarter-filled Hubbard model for organic conductors.

The DMRG has also been generalized to 1D random and disordered systems, and applied to the random antiferromagnetic and ferromagnetic Heisenberg chains, including quasiperiodic exchange modulation and a detailed study of the Haldane and Griffiths phase in these systems. Strongly disordered spin ladders have been considered in . It has also been used in disordered Fermi systems such as the spinless model. In particular, the transition from the Fermi glass to the Mott insulator and the strong enhancement of persistent currents in the transition was studied in correlated one-dimensional disordered rings. Disorder-induced crossover effects at quantum critical points were studied in .

### 4.4 Molecules and Quantum Chemistry

There have been several applications to molecules and polymers, such as the Pariser-Parr-Pople (PPP) Hamiltonian for a cyclic polyene (where long-range interactions are included), magnetic Keplerate molecules , molecular Iron rings and polyacenes considering long range interactions . It has also been applied to conjugated organic systems (polymers), adapting the DMRG to take into account the most important symmetries in order to obtain the desired excited states. Also conjugated one-dimensional semiconductors  have been studied, in which the standard approach can be extended to complex 1D oligomers where the fundamental repeat is not just one or two atoms, but a complex molecular building block. Relatively new fields of application are the calculation of dynamical properties in the Rubinstein-Duke model for reptons  and excitons in dendrimer molecules.

Recent attempts to apply DMRG to the ab initio calculation of electronic states in molecules have been successful. Here, DMRG is applied within the conventional quantum chemical framework of a finite basis set with non-orthogonal basis functions centered on each atom. After the standard Hartree-Fock (HF) calculation in which a Hamiltonian is produced within the orthogonal HF basis, DMRG is used to include correlations beyond HF, where each orbital is treated as a “site” in a 1D lattice. One important difference with standard DMRG is that, as the interactions are long-ranged, several operators must be kept, making the calculation somewhat cumbersome. However, very accurate results have been obtained in a check performed in a water molecule (keeping up to 25 orbitals and $`m200`$ states per block), obtaining an offset of 0.00024Hartrees with respect to the exact ground state energy, a better performance than any other approximate method.

In order to avoid the non-locality introduced in the treatment explained above, White introduced the concept of orthlets, local, orthogonal and compact wave functions that allow prior knowledge about singularities to be incorporated into the basis and an adequate resolution for the cores. The most relevant functions in this basis are chosen via the density matrix. An application based on the combination with the momentum version of DMRG is used in  to calculate the ground state of several molecules.

## 5 Dynamical correlation functions

The DMRG was originally developed to calculate static ground state properties and low-lying energies. However, it can also be useful to calculate dynamical response functions. These are of great interest in condensed matter physics in connection with experiments such as nuclear magnetic resonance (NMR), neutron scattering, optical absorption, photoemission, etc. We will describe three different methods in this Section.

A recent development for calculating response functions in single impurity systems in the presence of a magnetic field was done in  by using the DMRG within Wilson’s NRG to obtain the Green’s function.

An interesting extension of DMRG to tackle time-dependent quantum many-body systems out of equilibrium was considered by Cazalilla and Marston.

### 5.1 Lanczos and correction vector techniques

An effective way of extending the basic ideas of this method to the calculation of dynamical quantities is described in Ref.. It is important to notice here that due to the particular real-space construction, it is not possible to fix the momentum as a quantum number. However, we will show that by keeping the appropriate target states, a good value of momentum can be obtained.

We want to calculate the following dynamical correlation function at $`T=0`$:

$$C_A(tt^{})=\psi _0|A^{}(t)A(t^{})|\psi _0,$$
(4)
where $`A^{}`$ is the Hermitean conjugate of the operator $`A`$, $`A(t)`$ is the Heisenberg representation of $`A`$, and $`|\psi _0`$ is the ground state of the system. Its Fourier transform is:

$$C_A(\omega )=\underset{n}{}|\psi _n|A|\psi _0|^2\delta (\omega (E_nE_0)),$$
(5)
where the summation is taken over all the eigenstates $`|\psi _n`$ of the Hamiltonian $`H`$ with energy $`E_n`$, and $`E_0`$ is the ground state energy.

Defining the Green’s function

$$G_A(z)=\psi _0|A^{}(zH)^1A|\psi _0,$$
(6)
the correlation function $`C_A(\omega )`$ can be obtained as

$$C_A(\omega )=\frac{1}{\pi }\underset{\eta 0^+}{lim}\mathrm{Im}G_A(\omega +i\eta +E_0).$$
(7)
The function $`G_A`$ can be written in the form of a continued fraction:

$$G_A(z)=\frac{\psi _0|A^{}A|\psi _0}{za_0\frac{b_1^2}{za_1\frac{b_2^2}{z\mathrm{}}}}$$
(8)
The coefficients $`a_n`$ and $`b_n`$ can be obtained using the following recursion equations :

$$|f_{n+1}=H|f_na_n|f_nb_n^2|f_{n1}$$
(9)
where

$`|f_0`$ $`=`$ $`A|\psi _0`$
$`a_n`$ $`=`$ $`f_n|H|f_n/f_n|f_n,`$
$`b_n^2`$ $`=`$ $`f_n|f_n/f_{n1}|f_{n1};b_0=0`$ (10)
For finite systems the Green’s function $`G_A(z)`$ has a finite number of poles so only a certain number of coefficients $`a_n`$ and $`b_n`$ have to be calculated. The DMRG technique presents a good framework to calculate such quantities. With it, the ground state, Hamiltonian and the operator $`A`$ required for the evaluation of $`C_A(\omega )`$ are obtained. An important requirement is that the reduced Hilbert space should also describe with great precision the relevant excited states $`|\psi _n`$. This is achieved by choosing the appropriate target states. For most systems it is enough to consider as target states the ground state $`|\psi _0`$ and the first few $`|f_n`$ with $`n=0,1\mathrm{}`$ and $`|f_0=A|\psi _0`$ as described above. In doing so, states in the reduced Hilbert space relevant to the excited states connected to the ground state via the operator of interest $`A`$ are included. The fact that $`|f_0`$ is an excellent trial state, in particular, for the lowest triplet excitations of the two-dimensional antiferromagnet was shown in Ref. . Of course, if the number $`m`$ of states kept per block is fixed, the more target states considered, the less precisely each one of them is described. An optimal number of target states and $`m`$ have to be found for each case. Due to this reduction, the algorithm can be applied up to certain lengths, depending on the states involved. For longer chains, the higher energy excitations will become inaccurate. Proper sum rules have to be calculated to determine the errors in each case.

As an application of the method we calculate

$$S^{zz}(q,\omega )=\underset{n}{}|\psi _n|S_q^z|\psi _0|^2\delta (\omega (E_nE_0)),$$
(11)
for the 1D isotropic Heisenberg model with spin $`S=1/2`$.

The spin dynamics of this model has been extensively studied. The lowest excited states in the thermodynamic limit are the des Cloiseaux-Pearson triplets , having total spin $`S^T=1`$. The dispersion of this spin-wave branch is $`\omega _q^l=\frac{J\pi }{2}|\mathrm{sin}(q)|`$. Above this lower boundary there exists a two-parameter continuum of excited triplet states that have been calculated using the Bethe Ansatz approach  with an upper boundary given by $`\omega _q^u=J\pi |\mathrm{sin}(q/2)|`$. It has been shown , however, that there are excitations above this upper boundary due to higher order scattering processes, with a weight that is at least one order of magnitude lower than the spin-wave continuum.

In Fig. 2 we show the spectrum for $`q=\pi `$ and $`N=24`$ for different values of $`m`$, where exact results are available for comparison. The delta peaks of Eq. (11) are broadened by a Lorentzian for visualizing purposes. As expected, increasing $`m`$ gives more precise results for the higher excitations. This spectra has been obtained using the infinite-system method and more precise results are expected using the finite-system method, as described later.

In Fig. 3 we show the spectrum for two systems lengths and $`q=\pi `$ and $`q=\pi /2`$ keeping $`m=200`$ states and periodic boundary conditions. For this case it was enough to take 3 target states, i. e.  $`|\psi _0`$, $`|f_0=S_\pi ^z|\psi _0`$ and $`|f_1`$. Here we have used $`40`$ pairs of coefficients $`a_n`$ and $`b_n`$, but we noticed that if we considered only the first ($`10`$) coefficients $`a_n`$ and $`b_n`$, the spectrum at low energies remains essentially unchanged. Minor differences arise at $`\omega /J2`$. This is another indication that only the first $`|f_n`$ are relevant for the low energy dynamical properties for finite systems.

In the inset of Fig. 3 the spectrum for $`q=\pi /2`$ and $`N=28`$ is shown. For this case we considered 5 target states i. e.  $`|\psi _0`$, $`|f_0=S_{\pi /2}^z|\psi _0`$, $`|f_nn=1,3`$ and $`m=200`$. Here, and for all the cases considered, we have verified that the results are very weakly dependent on the weights $`p_l`$ of the target states (see Eq.(2)) as long as the appropriate target states are chosen. For lengths where this value of $`q`$ is not defined we took the nearest value.

Even though we are including states with a given momentum as target states, due to the particular real-space construction of the reduced Hilbert space, this translational symmetry is not fulfilled and the momentum is not fixed. To check how the reduction on the Hilbert space influences the momentum $`q`$ of the target state $`|f_0=S_q^z|\psi _0`$, we calculated the expectation values $`\psi _0|S_q^{}^zS_q^z|\psi _0`$ for all $`q^{}`$. If the momenta of the states were well defined, this value is proportional to $`\delta _{qq^{}}`$ if $`q0`$. For $`q=0`$, $`_rS_r^z=0`$.

The momentum distribution for $`q=\pi `$ is shown in Fig. 4 in a semilogarithmic scale where the $`y`$-axis has been shifted by .003 so as to have well-defined logarithms. We can see here that the momentum is better defined, even for much larger systems, but, as expected, more weight on other $`q^{}`$ values arises for larger $`N`$.

As a check of the approximation we calculated the sum rule

$$\frac{1}{4\pi ^2}_0^{\mathrm{}}𝑑\omega _{q=0}^{2\pi }S^{zz}(q,\omega )\psi _0|(S_{r=0}^z)^2|\psi _0=\frac{1}{4}$$
(12)
for $`N=28`$, 5 target states and $`m=200`$. We obtain a relative error of 0.86%.

Recently, important improvements to this method have been published : By considering the finite system method in open chains, Kühner and White obtained a higher precision in dynamical responses of spin chains. In order to define a momentum in an open chain and to avoid end effects, they introduce a filter function with weight centered in the middle of the chain and zero at the boundaries.

Recent applications of this method include the calculation of excitations in spin-orbital models (SU(4)) in a magnetic field, spin dynamics in models for cuprate spin ladders including cyclic spin exchange, optical conductivity of the ionic Hubbard model, excitations in the one-dimensional Bose-Hubbard model and the optical response in 1D Mott insulators.

In this section we have presented a method of calculating dynamical responses with DMRG. Although the basis truncation is big, this method keeps only the most relevant states and, for example, even by considering a $`0.1\%`$ of the total Hilbert space (for $`N=28`$ only $``$ 40000 states are kept) a reasonable description of the low energy excitations is obtained. We show that it is also possible to obtain states with well defined momenta if the appropriate target states are used.

#### 5.1.1 Correction vector technique

Introduced in Ref.  in the DMRG context and improved in Ref. , this method focuses on a particular energy or energy window, allowing a more precise description in that range and the possibility of calculating spectra for higher energies. Instead of using the tridiagonalization of the Hamiltonian, but in a similar spirit regarding the important target states to be kept, the spectrum can be calculated for a given $`z=w+i\eta `$ by using a correction vector (related to the operator $`A`$ that can depend on momentum $`q`$).

Following (6), the (complex) correction vector $`|x(z)`$ can be defined as:

$$|x(z)=\frac{1}{zH}A|\psi _0$$
(13)
so the Green’s function can be calculated as

$$G(z)=\psi _0|A^{}|x(z)$$
(14)
Separating the correction vector in real and imaginary parts $`|x(z)=|x^r(z)+i|x^i(z)`$ we obtain

$$((Hw)^2+\eta ^2)|x^i(z)=\eta A|\psi _0$$
(15)
and

$$|x^r(z)=\frac{1}{\eta }(wH)|x^i(z)$$
(16)
The former equation is solved using the conjugate gradient method. In order to keep the information of the excitations at this particular energy the following states are targeted in the DMRG iterations: The ground state $`|\psi _0`$, the first Lanczos vector $`A|\psi _0`$ and the correction vector $`|x(z)`$. Even though only a certain energy is focused on, DMRG gives the correct excitations for an energy range surrounding this particular point so that by running several times for nearby frequencies, an approximate spectrum can be obtained for a wider region .

A variational formulation of the correction vector technique which leads to more accurate excited energies and spectral weights has been developed in . It has been successfully applied to calculate the optical conductivity of Mott insulators.

### 5.2 Moment expansion

This method relies on a moment expansion of the dynamical correlations using sum rules that depend only on static correlation functions which can be calculated with DMRG. With these moments, the Green’s functions can be calculated using the maximum entropy method.

The first step is the calculation of sum rules. As an example, and following , the spin-spin correlation function $`S^z(q,w)`$ of the Heisenberg model is calculated where the operator $`A`$ of Eq. (4) is $`S^z(q)=N^{1/2}S^z(l)\mathrm{exp}(iql)`$ and the sum rules are:

$`m_1(q)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dw}{\pi }}{\displaystyle \frac{S^z(q,w)}{w}}={\displaystyle \frac{1}{2}}\chi (q,w=0)`$
$`m_2(q)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dw}{\pi }}w{\displaystyle \frac{S^z(q,w)}{w}}={\displaystyle \frac{1}{2}}S^z(q,t=0)`$
$`m_3(q)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dw}{\pi }}w^2{\displaystyle \frac{S^z(q,w)}{w}}={\displaystyle \frac{1}{2}}[[H,S^z(q)],S^z(q)]`$ (17)
$`=`$ $`2[1\mathrm{cos}(q)]{\displaystyle \underset{i}{}}S_i^+S_{i+1}^{}+S_i^{}S_{i+1}^+`$
where $`\chi (q,w=0)`$ is the static susceptibility. These sum rules can be easily generalized to higher moments:

$`m_l(q)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dw}{\pi }}w^{l1}{\displaystyle \frac{S^z(q,w)}{w}}`$ (18)
$`=`$ $`{\displaystyle \frac{1}{2}}[[H,\mathrm{},[H,S^z(q)]\mathrm{}],S^z(q)]`$
for $`l`$ odd. A similar expression is obtained for $`l`$ even, where the outer square bracket is replaced by an anticommutator and the total sign is changed. Here $`H`$ appears in the commutator $`l2`$ times.

Apart from the first moment which is given by the static susceptibility, all the other moments can be expressed as equal time correlations (using a symbolic manipulator). The static susceptibility $`\chi `$ is calculated by applying a small field $`h_q_in_i\mathrm{cos}(qi)`$ and calculating the density response $`n_q=1/N_in_q\mathrm{cos}(qi)`$ with DMRG. Then $`\chi =n_q/h_q`$ for $`h_q0`$. These moments are calculated for several chain lengths and extrapolated to the infinite system. Once the moments are calculated, the final spectra is constructed via the Maximum Entropy method (ME), which has become a standard way to extract maximum information from incomplete data (for details see Ref.  and references therein). Reasonable spectra are obtained for the XY and isotropic models, although information about the exact position of the gaps has to be included. Otherwise, the spectra are only qualitatively correct. This method requires the calculation of a large amount of moments in order to get good results: The more information given to the ME equations, the better the result.

### 5.3 Finite temperature dynamics

In order to include temperature in the calculation of dynamical quantities, the Transfer Matrix RG described above (TMRG) was extended to obtain imaginary time correlation functions. After Fourier transformation in the imaginary time axis, analytic continuation from imaginary to real frequencies is done using maximum entropy (ME). The combination of the TMRG and ME is free from statistical errors and the negative sign problem of Monte Carlo methods. Since we are dealing with the transfer matrix, the thermodynamic limit can be discussed directly without extrapolations. However, in the present scheme, only local quantities can be calculated.

A systematic investigation of local spectral functions is done in Ref.  for the anisotropic Heisenberg antiferromagnetic chain. The authors obtain good qualitative results especially for high temperatures but a quantitative description of peaks and gaps are beyond the method, due to the severe intrinsic limitation of the analytic continuation. This method was also applied with great success to the 1D Kondo insulator. The temperature dependence of the local density of states and local dynamic spin and charge correlation functions were calculated.

## 6 Conclusions

We have presented here a very brief description of the Density Matrix Renormalization Group technique, its applications and extensions. The aim of this article is to give the unexperienced reader an idea of the possibilities and scope of this powerful, though relatively simple method. The experienced reader can find here an extensive (however incomplete) list of references covering most applications to date using DMRG in a great variety of fields such as Condensed Matter, Statistical Mechanics and High Energy Physics.

## Acknowledgments

The author acknowledges hospitality at the Centre de Recherches Mathematiques, University of Montreal and at the Physics Department of the University of Buenos Aires, Argentina, where this work has been performed. We thank S. White for a critical reading of the manuscript and all those authors that have updated references and sent instructive comments. K. H. is a fellow of CONICET, Argentina. Grants: PICT 03-00121-02153 and PICT 03-00000-00651.


format: text

text:
# Implications of WMAP Observations On the Population III Star Formation Processes

## 1 Introduction

The cosmological reionization process is likely to be quite complex (e.g., Barkana & Loeb 2001; Madau 2002; Wyithe & Loeb 2003a; Mackey, Bromm, & Hernquist 2002; Cen 2003; Venkatesan, Tumlinson, & Shull 2003). It seems rather remarkable that we can now start observationally probing this process. A major milestone was laid down by the recent observations of high redshift quasars from the Sloan Digital Sky Survey (SDSS), indicating that the final reionization episode came to completion at $`z6`$ (e.g., Fan et al. 2001; Becker et al. 2001; Barkana 2002; Cen & McDonald 2002; Litz et al. 2002). Given the final reionization epoch and the fact that the density fluctuations at that epoch is sufficiently well constrained in the context of the standard cold dark matter model (Bahcall et al. 1999), it becomes possible to show (Cen 2003), under reasonable assumptions, that the universe has been reionized twice, first reionized at $`z15`$ by metal-free, massive Population III (Pop III) stars and second reionized at $`z=6`$, provided that Pop III IMF is top-heavy. In this picture a larger Thomson optical depth, $`\tau _e0.07`$, was expected.

Another major milestone was set very recently by the WMAP polarization observations, which probes the ionization state of the gas in the universe at high redshift during the cosmological reionization period. The WMAP observations measured $`\tau _e=0.17\pm 0.04`$ ($`68\%`$) (Kogut et al. 2003). One can immediately draw the conclusion that the universe did not just experience an impulsive reionization event at $`z6`$, which would have yielded $`\tau _e0.03`$, in the standard cosmological constant dominated cold dark matter model. Evidently, the universe became (and possibly had to stay) significantly ionized from a much higher redshift, probably from $`z15`$. In this Letter we compute the detailed reionization process and compare to WMAP observations, investigating what can be learned about Pop III star formation processes at high redshift.

## 2 Thomson Optical Depth $`\tau _e`$

A spatially flat cold dark matter cosmological model with $`\mathrm{\Omega }_M=0.27`$, $`\mathrm{\Omega }_b=0.047`$, $`\mathrm{\Lambda }=0.73`$, $`H_0=72`$km/s/Mpc, $`n_s=0.99`$ and $`\sigma _8=0.90`$, symbolized as LCDMAP model, is used. This model is the best fit model constrained by WMAP observations with a fixed power-law index (Spergel et al. 2003). To set a plausible upper bound we also compute a model with $`n_s=1.03`$ and $`\sigma _8=1.0`$, denoted as LCDMAP+ model, which is consistent with current WMAP data (Spergel et al. 2002) and maximizes the small-scale structures responsible for the reionization process. We employ an efficient method to compute the coupled evolution of gas and star formation in the high redshift universe, described in §4 of Cen (2003), which simultaneously computes the dynamic evolution of all gaseous phases including HII regions, HI regions and partially ionized regions. Such a careful treatment is necessary given the various short time scales involved at high redshift for important processes, including recombination and cooling. Throughout, Pop III stars are defined to be metal-free stars with a top-heavy initial mass function (IMF) and Pop III galaxies are defined to be galaxies where Pop III stars form. We define “minihalos” as those whose virial temperature is less than $`1\times 10^4`$K where only $`\mathrm{H}_2`$ cooling is possible in the absence of metals, and large halos as those with virial temperature above $`1\times 10^4`$K capable of cooling via atomic lines.

The most important factors determining the ionization process are $`\mathrm{\Psi }c_{}\times f_{\mathrm{esc}}`$ and $`C`$, where $`c_{}`$ is the star formation efficiency, defined as the fraction of gas in halos formed into stars; $`f_{\mathrm{esc}}`$ is the fraction of ionizing photons that escapes from galaxy halos and gets into the IGM; $`C`$ is the clumping factor of the intergalactic medium (IGM). Basically, $`\mathrm{\Psi }`$ is the source of ionizing photons from stars in galaxies and $`C`$ determines the sink of ionizing photons due to the IGM. The competition between the two primarily dictates the ionization state of the IGM. The initial phase of ionization front propagating through neutral medium does not play a significant role near the end of the reionization process, and cosmological effects are unimportant throughout the reionization period. But all these processes are included in the calculation. Theoretically, $`C`$ is reasonably well determined, if we adopt the standard cold dark matter cosmological model. In addition, $`C`$ may be constrained directly by observations of quasar absorption spectra at redshifts close to $`z=6`$ (e.g., Cen & McDonald 2002; Djorgovski et al. 2001; Cen & Haiman 2000). We set the IGM clumping factor at $`z=6`$ to be $`50.0`$, according to Gnedin & Ostriker (1997), which requires $`C_{\mathrm{halo}}=702`$ with the adopted cosmological model (see Cen 2003 for discussion). Results do not sensitively depend on $`C`$ at $`z=6`$; for example, a change in $`C`$ at $`z=6`$ by $`25\%`$ only causes $`5\%`$ change in $`\tau _e`$. Both theoretically and observationally, our direct knowledge of $`\mathrm{\Psi }`$ for high redshift galaxies is very little. However, as pointed out in Cen (2003), the fact that the universal reionization ends just at $`z6`$, as observations seem to indicate (e.g., Fan et al. 2001; Becker et al. 2001; Barkana 2002; Cen & McDonald 2002; Litz et al. 2002), tightly constrains $`\mathrm{\Psi }`$, knowing $`C`$. Such a normalization point is proven to be very powerful and unambiguously determines $`\mathrm{\Psi }`$ at $`z6`$ (see Figure 7 of Cen 2003) for Pop II stars, which are likely responsible for completing the long reionization process at $`z6`$. For all the models computed below, the normalization at $`z=6`$ requires $`\mathrm{\Psi }(II)c_{}(II)\times f_{esc}(II)=0.0150.017`$ for Pop II galaxies at $`z6`$. For models computed we assume $`\mathrm{\Psi }(II)`$ to be constant with redshift, except indicated otherwise. On the other hand, $`c_{}(III)`$ and $`f_{esc}(III)`$ are less constrained.

We compute seven models with varying parameters for Pop III star formation processes and cosmological parameters, listed in Table 1, to show the possible range for the Thomson optical depth due to the IGM at reionization. In Table 1, $`c_{}(H_2,III)`$ indicates the star formation efficiency in minihalos with virial temperature less than $`8,000`$K, where $`\mathrm{H}_2`$ cooling dominates for metal-free gas; $`c_{}(HI,III)`$ indicates the star formation efficiency in large halos with more efficient atomic cooling; $`f_{crit}`$ indicates the fractional amount of gas formed into Pop III stars at the transitional epoch from Pop III to Pop II (which is assumed to determine the metallicity of the IGM; see Cen 2003); $`f_{esc}(III)`$ is the mean fraction of ionizing photons produced by galaxies that escapes from the halos of Pop III galaxies and into the IGM; the last column $`\tau _e`$ is the integrated (from $`z=0`$ to $`z=1000`$) Thomson optical depth due to the IGM, where multiple entries correspond to respective multiple entries for $`f_{esc}(III)`$.

Models #1,2 represents the case where the IMF at high redshift is assumed to remain as the Salpeter (1955) function with the same low-mass cut-off ($`0.1\mathrm{M}_{}`$), as at $`z=6`$, which is the likely minimum model in terms of $`\tau _e`$, for LCDMAP and LCDMAP+ model, respectively. Models #3,4 have a Pop III star formation efficiency in minihalos $`c_{}(H_2,III)=0.002`$, suggested by numerical simulations of Abel, Bryan & Norman (2002) with a reasonable star formation efficiency for Pop III galaxies with large halos of $`c_{}(HI,III)=0.10`$, for LCDMAP and LCDMAP+ model, respectively. Models #5,6 assume a higher Pop III star formation efficiency in minihalos $`c_{}(H_2,III)=0.01`$, for LCDMAP and LCDMAP+ model, respectively; in this case, ionizing photons from Pop III minihalos dominate over those from Pop III large halos for the first reionization process, regardless of the value of $`c_{}(HI,III)`$. Model #7 is similar to Models #1, except that in the former we arbitrarily set the transition epoch from Pop III to Pop II at the time when a fraction $`f_{\mathrm{crit}}=5\times 10^4`$ of total gas has formed Pop III stars instead of $`f_{\mathrm{crit}}=1\times 10^4`$ in the latter set. We note that $`f_{\mathrm{crit}}=1\times 10^4`$ is derived based on the assumption that most of the Pop III metal-free massive stars experience supernovae and most of the ejected metals get out of galaxies to enrich the IGM (Oh et al. 2001; Cen 2003). If, for some reason, the fractional amount of metals enriching the IGM is significantly smaller than the amount of gas formed into massive stars, the efficiency of metal enrichment by Pop III stars would be lower; Model #7 serves to illustrate this possibility. Possible physical processes that may reduce the metal enrichment efficiency include a large fraction of Pop III stars being very massive stars ($`M300\mathrm{M}_{}`$) which may collapse to black holes instead of exploding as supernovae (Rakavy, Shaviv, & Zinamon 1967; Bond, Arnett, & Carr 1984; Glatzel, Fricke, & El Eid 1985; Woosley 1986) or not all ejected metals are able to be transported to IGM. We note that Model #1 is close to the model used in Cen (2003), having perhaps more conventional (thus thought to be more reasonable) values for the relevant parameters. We also test a variant of Model #1 by using $`\mathrm{\Psi }(II)(1+z)^{1/2}`$ instead of being constant with redshift and find $`\tau _e=0.08`$ for that case versus $`0.07`$ for a constant $`\mathrm{\Psi }(II)`$ (Model #1).

From Table 1 we can draw several conclusions. First, without Pop III stars (Models #1,2), it is likely that $`\tau _e<0.09`$ for both LCDMAP and LCDMAP+ models, inconsistent with WMAP results at $`2\sigma `$ level (Kogut et al. 2003). The conclusions reached here with regard to the requirement of a top-heavy IMF for Pop III stars in order to explain the observed high $`\tau _e`$ value are consistent with those by Wyithe & Loeb (2003b), Haiman & Holder (2003) and Sokasian et al. (2003). Second, for reasonable ranges for the star formation efficiency \[$`c_{}(HI,III)0.1`$, $`c_{}(H_2,III)0.002`$\] and ionizing photon escape fraction ($`f_{esc}(III)0.20`$) from Pop III galaxies, we expect that $`\tau _e0.12`$ for LCDMAP model (Model #3) and $`\tau _e0.15`$ for LCDMAP+ model (Model #4). Third, prolonging the Pop III era by making the metal enrichment of the IGM from Pop III galaxies less efficient increases the Thomson optical depth incrementally. Finally, while $`\tau _e=0.17`$ is the mean value determined by WMAP observations, it seems hard to reach in the LCDMAP model. There are two possible ways to achieve a Thomson optical depth as high as $`\tau _e0.17`$: either the ionizing photon escape fraction is large ($`f_{esc}(III)0.3`$) and the underlying cosmological model has a blue power spectrum with a positive tilt to $`n1.03`$ (Model #5) or Pop III star formation efficiency in minihalos is substantially larger than current simulations seem to indicate, requiring $`c_{}(H_2,III)>0.01`$.

Figures 1,2 show the reionization history and cumulative Thomson optical depth, respectively, for four representative models listed in Table 1. We see that without Pop III stars (Model #1) reionization event only occurs once at $`z6`$ and the reionization history is monotonic. With Pop III stars the first reionization occurs at $`z=1520`$ (Models #3,4). In the last case (Model #5) where Pop III star formation efficiency in minihalos is substantially larger than suggested by current simulations of Abel, Bryan & Norman (2002), $`0.01`$ versus $`0.002`$, the first reionization occurs at $`z2025`$. While both Model #4 and Model #5 have comparable total Thomson optical depth, $`\tau _e=0.150.16`$, their respective reionization histories at early times ($`z15`$) are somewhat different. We find that for Models #4, minihalos and large halos contribute to 62% and 38% of all ionizing photons to reionize the universe at the first time, by which 9 ionizing photons per baryon have been produced; the corresponding numbers for Model #5 are 61%, 39% and 11 photons per baryon. In contrast, for Model #3, minihalos and large halos contribute to 20% and 80% of all ionizing photons to the first reionization and 6 ionizing photons per baryon have been produced by then. We see, as expected, that an earlier first reionization requires a larger number of ionizing photons per baryons and a larger contribution from minihalos appears likely if the first reionization occurs at $`z20`$ in the context of the currently favored $`\mathrm{\Lambda }`$ cosmology. It is possible that future CMB experiments such as Planck survey might be able to distinguish between these two scenarios (e.g., Holder et al. 2003). We will defer a more detailed analysis of this aspect to a separate paper.

## 3 Discussion

Since it appears that a significant fraction ($`f_{esc}(III)0.1`$) of ionizing photons produced by Pop III galaxies is required to escape from galaxies in order to sufficiently ionize the universe to be consistent with WMAP observations, it is worthwhile to check what kind of escape fraction is possible. Here we give an order of magnitude assessment of the situation. If a fraction, $`c_{}(\mathrm{III})`$, of gas in a high redshift Pop III galaxy formed into very massive stars emitting at $`1.6\times 10^{48}`$ hydrogen ionizing photons per solar mass of stars per second (e.g., Bromm, Kudritzki, & Loeb 2001), then we can obtain the ratio of the total number of hydrogen ionizing photons emitted during the $`3\times 10^6`$ yrs of lifetime of the Pop III galaxy (assumed to be equal to the lifetime of the Pop III massive stars) to the number of baryons in the galaxy,

$$\frac{N_{\mathrm{ph}}}{N_\mathrm{H}}=1.3\times 10^5c_{}(\mathrm{III}).$$
(1)
It can also be shown that the ratio of the number of recombinations of hydrogen atoms (assuming a temperature of $`1.5\times 10^4`$K at $`z=20`$) to the number of baryons in the galaxy is

$$\frac{N_{\mathrm{rec}}}{N_\mathrm{H}}=0.061C_\mathrm{g},$$
(2)
where $`C_\mathrm{g}`$ is the effective clumping factor of the gas in the galaxy, taking into account all possible effects including shielding and clumpiness. Thus, the ratio of the total number of ionizing photons to the total number of recombinations of hydrogen atoms during the lifetime of the galaxy at redshift $`z=20`$ is

$$\frac{N_{\mathrm{ph}}}{N_{\mathrm{rec}}}=2.1\times 10^6C_\mathrm{g}^1c_{}(\mathrm{III}).$$
(3)
The time it takes for the ionization front from very massive Pop III stars to break out is roughly

$$t_{\mathrm{IF}}=\frac{4\pi r^3n}{L_{\mathrm{ph}}}=6.8(\frac{M_{}}{100\mathrm{M}_{}})^1(\frac{r}{1\mathrm{p}\mathrm{c}})^3(\frac{n}{10^4\mathrm{cm}^3})\mathrm{yrs},$$
(4)
which is clearly much shorter than the lifetime of the galaxy. In Equation (4) $`M_{}`$ is the amount of Pop III stars formed, $`r`$ the size of the galaxy disk and $`n`$ the interstellar medium density. Thus, ionization of the gas inside the galaxy may be treated as instantaneous and the ratio in Equation (3) is an instructive gauge.

If $`c_{}(HI,\mathrm{III})0.1`$, as appears to be required (as discussed in §2) in order to produce a large enough $`\tau _e`$ ($`0.1`$) to be consistent with WMAP observations, then Equation (3) implies that $`C_g`$ has to be less than $`2.1\times 10^5`$ in order for a significant fraction of ionizing photons to escape from galaxy halos. On one hand, if the clumping factor is dominated by gas in the halo with overdensity of order $`100`$, $`C_g`$ could be less than $`2.1\times 10^5`$. On the other hand, the gas on the disk or at radii much smaller than the virial radius would be much denser and their contribution to the clumping factor would critically depend on the density, geometry and self-shielding effects. A reliable estimate of the escape fraction for ionizing photons is unattainable without detailed modeling. Nevertheless, it appears that it is not entirely impossible to have a significant fraction of ionizing photons escape under special conditions. Possible favorable scenarios include a large fraction of star formation activities occur at a substantial height above the galactic disk and/or off-center within star formation sites embedded in very clumpy interstellar medium.

It is clearly seen from the above analysis that the smaller the value of $`c_{}(\mathrm{III})`$, the less likely that a significant fraction of ionizing photon can escape. Consequently, a high Pop III star formation efficiency may be the key to a high overall ionizing photon production rate in Pop III galaxies with respect to the receiving IGM, with the latter depending on the former in a probably nonlinear way. The question is then: how high can $`c_{}(\mathrm{III})`$ be?

The simulation by Abel, Bryan & Norman (2002; see also Bromm, Coppi, & Larson 2002) for a halo of mass $`10^6\mathrm{M}_{}`$ with $`6\times 10^4\mathrm{M}_{}`$ baryons forms a massive star of mass $`100\mathrm{M}_{}`$ at the center. This would give a star formation efficiency of $`0.0017`$, taken at face value, which is much smaller than required value of order $`0.01`$. For a discussion why Pop III star formation in minihalos may continue unimpeded throughout the Pop III era, see Cen (2003). The value of $`0.1`$ required for Pop III stars in large halos to ionize the universe at $`z=1520`$ is still higher by a factor of ten. It should be pointed out that that the central protostar in the Abel et al. (2002) simulation appears to be still accreting matter rapidly at the end of their simulation, so it is possible that the stellar mass may grow substantially. We also note, however, that their simulated galaxy has a virial temperature below $`1\times 10^4`$K hence does not possess the efficient atomic cooling. It is possible that, for galaxies with efficient atomic cooling, star formation efficiency may be much higher than what Abel et al. (2002) simulation indicates. Detailed simulations will be invaluable.

Finally, we note that, if star formation in cooling shells produced by exploding massive Pop III stars, the “Pop II.5” stars proposed by Mackey et al. (2003), is efficient, they may make a non-negligible contribution to the pool of ionizing photons. These stars might form at off-center locations and thus might possess a relatively larger ionizing photon escape fraction. The counteracting factor is that these stars might be substantially less massive than the first generation metal-free stars, since the metallicity of the gas in the cooling shells may be quite high, hence are less efficient ionizing photon emitters.

## 4 Conclusions

Adopting the best fit standard cold dark matter model with a fixed power-law index by WMAP observations with $`\mathrm{\Omega }_M=0.27`$, $`\mathrm{\Omega }_b=0.047`$, $`\mathrm{\Lambda }=0.73`$, $`H_0=72`$km/s/Mpc, $`n_s=0.99`$ and $`\sigma _8=0.90`$ (Spergel et al. 2003) and based on the observed Thomson optical depth due to intergalactic medium by WMAP (Kogut et al. 2003), we are able to draw several relatively secure conclusions with regard to Pop III star formation processes at very high redshift.

(1) The combination of the normal Salpeter IMF for Pop III metal-free stars and the absence of a dramatic upturn in the star formation efficiency and/or ionizing photon escape fraction at high redshift ($`z>6`$) would produce a Thomson optical depth due to IGM at reionization of $`\tau _e0.09`$, inconsistent with the observed $`\tau _e=0.17\pm 0.04`$ ($`68\%`$) (Kogut et al. 2003) at $`2\sigma `$ level.

(2) A top-heavy IMF for the Pop III metal-free stars and plausible star formation efficiency and ionizing photon escape, as gauged by the corresponding values for Pop II galaxies required in order to achieve the second reionization finale at $`z6`$ yield $`\tau _e0.12`$;

(3) In the event that the metal enrichment efficiency of the intergalactic medium by Pop III stars is very low thus Pop III era is prolonged, one may be able to obtain $`\tau _e=0.15`$.

(4) It seems quite improbable to reach $`\tau _e0.17`$ even with very massive Pop III metal-free stars, unless (i) the cosmological model power index $`n`$ is positively tilted to $`n1.03`$ and/or (ii) Pop III star formation in minihalos with molecular hydrogen cooling has an efficiency $`c_{}(H_2,III)>0.01`$ (still requiring ionizing photon escape fraction greater than $`0.3`$) or (iii) alternatively, there may be unknown, non-stellar ionizing sources at very high redshift in the context of the $`\mathrm{\Lambda }`$CDM cosmology.

It is expected that the errorbars on the detected Thomson optical depth should shrink significantly with time using more data. Thus, if the current observed value of Thomson optical depth ($`\tau _e0.09`$ \- currently $`2\sigma `$ lower bound) withstands future data, we will have strong observational evidence for the existence of massive metal-free Pop III stars. One possible outcome will be that $`0.09\tau _e0.12`$ ($`1.252.0\sigma `$ away from the current observed central value), in which case no drastic increase in the star formation efficiency in Pop III galaxies would be needed. On the other hand, if $`\tau _e0.17`$ is firmed up in the future, there will be a degeneracy between possible competing scenarios. Future CMB experiments such as the Planck surveyor may be able to distinguish between them. On the theoretical side, more accurate calculations of Pop III metal-free star formation in both minihalos and large halos will be extremely useful to shed light on these vital parameters.

I thank Zoltan Haiman, Lars Hernquist, Eve Ostriker, David Spergel and Jonathan Tan for useful discussion. This research is supported in part by NSF grant AST-0206299.


text:
# Rapidity Asymmetry in High-energy 𝑑+𝐴 Collisions

## I Introduction

Recent experiments  at the Relativistic Heavy-ion Collider (RHIC) have shown a significant suppression of high $`p_T`$ hadron spectra in central $`Au+Au`$ collisions that was predicted  as a consequence of parton energy loss or jet quenching in dense matter. In addition, the same mechanism is predicted to produce azimuthal anisotropy in high $`p_T`$ hadron spectra  that was also observed in experiments at RHIC . This is a dramatic departure from the heavy-ion collisions at the SPS energies where no significant suppression of high $`p_T`$ spectra is observed . Since theoretical studies  of parton propagation in a dense medium all show that the parton energy loss induced by multiple scattering is proportional to the gluon density, RHIC data thus indicate an initial gluon density in central $`Au+Au`$ collisions at the RHIC energies that is much higher than that in a large cold nucleus .

More precise extraction of the parton energy loss from the final hadron suppression in $`A+A`$ collisions, however, requires the understanding of normal nuclear effects in $`p+A`$ collisions. As pointed out by many early  and recent  studies, the high $`p_T`$ hadron spectra can also be modified by initial multiple scatterings in $`p+A`$ and $`A+A`$ collisions giving rise to the observed Cronin effect. Within a multiple scattering model, the high $`p_T`$ hadron spectra are normally enhanced relative to $`p+p`$ collisions, except in the kinematic region where the EMC effect  is important (The EMC effect is the depletion of parton distributions in $`x0.20.8`$ in nuclei caused by the nuclear binding effect). Such a normal nuclear enhancement will not affect the interpretation of the hadron suppression in central $`A+A`$ collisions as a consequence of jet quenching, though inclusion of it is important for more precise extraction of parton energy loss and the initial gluon density.

A very different mechanism was recently proposed for the observed high $`p_T`$ hadron suppression based on the parton saturation model  which also predicts a similar suppression in $`p+A`$ collisions at RHIC, in contrast to the predicted enhancement by the multiple parton scattering model. While such a model is not yet checked against the existing $`p+A`$ collisions for energies up to $`\sqrt{s}=40`$ GeV where enhancement of hadron spectra for $`p_T>2`$ GeV/$`c`$ has been successfully explained by the multiple scattering model, the up-coming data of $`d+A`$ collisions at RHIC will attest the relevance of parton saturation at the RHIC energies.

In this letter, we will point out an additional feature in the rapidity dependence of the Cronin enhancement due to multiple parton scattering in $`d+A`$ collisions. Such a rapidity dependence was also studied recently by Vitev . However, we predict here a unique rapidity asymmetry of high $`p_T`$ hadron spectra due to stronger Cronin enhancement in the forward (projectile) region as a result of the transverse momentum broadening of the initial partons inside the projectile. The shape of the rapidity asymmetry will also depend on the nuclear modification of the parton distributions inside a nucleus, in particular at small $`x`$. As one decreases $`p_T`$, the parton shadowing will reduce the hadron spectra in the forward region, thus changing the rapidity asymmetry. When soft and coherent interactions become dominant at very low $`p_T<1`$ GeV/$`c`$, the shape of the rapidity asymmetry will be reversed because of the strong suppression of hadron production in the projectile region relative to a superposition of binary $`p+p`$ collisions. We will calculate the rapidity asymmetry within a perturbative QCD (pQCD) parton model and study the effect of nuclear modification of parton distributions.

In a Glauber multiple parton scattering model, large $`p_T`$ spectra are generally enhanced relative to the binary model of hard scattering. It can be shown  that a combination of absorptive corrections and the power-law behavior of perturbative parton cross section leads to a nuclear enhancement at high $`p_T`$ that decreases as $`1/p_T^2`$. The same absorptive processes suppress the spectra relative to the binary model at low $`p_T`$, where soft processes dominate and the $`p_T`$ spectra deviate from a power-law behavior. Since the pQCD model cannot be applied to soft processes, one has to resort to phenomenological models like the string model in which coherent particle production is modeled by a string excitation for each participant nucleon. Shown in Fig. 1 is the rapidity distribution (integrated over transverse momentum) of charged hadrons in minimum-biased $`d+Au`$ collisions at $`\sqrt{s}=200`$ GeV from the HIJING Monte Carlo model , which employs the string model for soft and coherent interactions. Because most of the soft particles are produced through string fragmentation, their number should then be proportional to the number of participants. The rapidity distributions of hadrons from the string should also follow their parent nucleons. Since there are more target nucleon participants than projectile nucleons in $`d+A`$ collisions, the hadron rapidity distribution is thus asymmetrical with respect to $`\eta =0`$. For high $`p_T`$ hadrons, the underlying processes are hard parton scatterings. The multiplicity should be approximately proportional to the number of binary scatterings and the rapidity distributions should be approximately symmetric with respect to $`\eta =0`$. This is roughly the case for high $`p_T`$ hadrons in HIJING model as shown in Fig. 1. One, however, can notice some asymmetric effect in the rapidity distribution of large $`p_T`$ hadrons. This is partially due to the coherence between transverse jets and the beam remnants, which is responsible for the asymmetric pedestal effect underlying a jet event in $`p+A`$ collisions. In addition, energy and quark number conservation will also cause some asymmetric effects. These effects are most important in the large rapidity region. It is expected that they are small in the central rapidity region.

Though the HIJING model has incorporated hard processes, it has not included transverse momentum broadening of initial partons. For this purpose, we use the pQCD model as employed in Ref. . We will use a lowest order (LO) pQCD-inspired parton model in which the inclusive particle production cross section in $`pp`$ collisions is given by

$`{\displaystyle \frac{d\sigma _{pp}^h}{dyd^2p_T}}`$ $`=`$ $`K{\displaystyle \underset{abcd}{}}{\displaystyle 𝑑x_a𝑑x_bd^2k_{aT}d^2k_{bT}g_p(k_{aT},Q^2)g_p(k_{bT},Q^2)}`$ (2)
$`f_{a/p}(x_a,Q^2)f_{b/p}(x_b,Q^2){\displaystyle \frac{D_{h/c}^0(z_c,Q^2)}{\pi z_c}}{\displaystyle \frac{d\sigma }{d\widehat{t}}}(abcd),`$
where $`D_{h/c}^0(z_c,Q^2)`$ is the fragmentation function of parton $`c`$ into hadron $`h`$ as parameterized in Ref.  from $`e^+e^{}`$ data, and $`z_c`$ is the momentum fraction of a parton jet carried by a produced hadron. The $`K1.5`$ (at $`\sqrt{s}=200`$ GeV) factor is used to account for higher order QCD corrections to the jet production cross section. The parton distributions $`f_{a/N}(x,Q^2)`$ in a nucleon are given by the MRS D$`^{}`$ parameterization . The initial transverse momentum distribution $`g_N(k_T,Q^2)`$ is assumed to have a Gaussian form,

$$g_N(k_T,Q^2)=\frac{1}{\pi k_T^2_N}e^{k_T^2/k_T^2_N}.$$
(3)
Following , we choose a $`Q`$-dependent average initial transverse momentum,

$$k_T^2_N(Q^2)=1.2(\mathrm{GeV}^2)+0.2\alpha _s(Q^2)Q^2,$$
(4)
which should include both the intrinsic and pQCD radiation-generated transverse momentum in this LO calculation. The form of the $`Q`$-dependence and the parameters are chosen to reproduce the experimental data , especially at low energies. Following the same approach as in Refs. , we choose $`Q^2`$ to be $`Q^2=2\widehat{s}\widehat{t}\widehat{u}/(\widehat{s}^2+\widehat{t}^2+\widehat{u}^2)`$.

To take account of multiple initial-state scattering, we assume that the inclusive differential cross section for large $`p_T`$ particle production is still given by a single hard parton-parton scattering. However, due to multiple parton scattering prior to the hard processes, we consider the initial transverse momentum $`k_T`$ of the beam partons to be broadened. Assuming that each scattering provide a $`k_T`$ kick which also has a Gaussian distribution, we can effectively change the width of the initial $`k_T`$ distribution. Then the single inclusive particle cross section in minimum-biased $`p+A`$ collisions is,

$`{\displaystyle \frac{d\sigma _{pA}^h}{dyd^2p_T}}`$ $`=`$ $`K{\displaystyle \underset{abcd}{}}{\displaystyle d^2bt_A(b)𝑑x_a𝑑x_bd^2k_{aT}d^2k_{bT}g_A(k_{aT},Q^2,b)g_p(k_{bT},Q^2)}`$ (6)
$`f_{a/p}(x_a,Q^2)f_{b/A}(x_b,Q^2,b){\displaystyle \frac{D_{h/c}^0(z_c,Q^2)}{\pi z_c}}{\displaystyle \frac{d\sigma }{d\widehat{t}}}(abcd),`$
where $`t_A(b)`$ is the nuclear thickness function normalized to $`d^2bt_A(b)=A`$. We will use the Woods-Saxon form of nuclear distribution for $`t_A(b)`$ throughout this paper unless specified otherwise. The parton distribution per nucleon inside the nucleus (with atomic mass number $`A`$ and charge number $`Z`$) at an impact parameter $`b`$,

$$f_{a/A}(x,Q^2,b)=S_{a/A}(x,b)\left[\frac{Z}{A}f_{a/p}(x,Q^2)+(1\frac{Z}{A})f_{a/n}(x,Q^2)\right],$$
(7)
is assumed to be factorizable into the parton distribution in a nucleon $`f_{a/N}(x,Q^2)`$ and the nuclear modification factor $`S_{a/A}(x,b)`$, for which we take both the new HIJING  and EKS  parameterizations. The initial parton transverse momentum distribution inside a projectile nucleon going through the target nucleus at an impact parameter $`b`$ is still a Gaussian with a broadened width

$$k_T^2_A(Q^2)=k_T^2_N(Q^2)+\delta ^2(Q^2)(\nu _A(b)1).$$
(8)
The broadening is assumed to be proportional to the number of scattering $`\nu _A(b)`$ the projectile suffers inside the nucleus, which is assumed to given by

$$\nu _A(b)=\sigma _{NN}t_A(b)=\sigma _{NN}\frac{3A}{2\pi R_A^2}\sqrt{1b^2/R_A^2}$$
(9)
in a hard sphere nuclear distribution, where $`R_A=1.12A^{1/3}`$ fm and $`\sigma _{NN}`$ is the inelastic nucleon-nucleon cross section. We also assume that $`k_T`$ broadening during each nucleon-nucleon collision $`\delta ^2`$ also depends on the hard momentum scale $`Q=P_T^{\mathrm{jet}}`$ in the parameterized form,

$$\delta ^2(Q^2)=0.225\frac{\mathrm{ln}^2(Q/\mathrm{GeV})}{1+\mathrm{ln}(Q/\mathrm{GeV})}\mathrm{GeV}^2/c^2,$$
(10)
which is chosen to best fit the existing experimental data in $`p+A`$ collisions  up to $`\sqrt{s}=40`$ GeV. The predictive power of this model lies in the energy and flavor dependence of the hadron spectra. It is straightforward to also calculate hadron spectra in $`d+A`$ collisions, incorporating the initial $`k_T`$ broadening and parton shadowing. We have tried different forms of nuclear distribution for deuteron and find little difference in the final results. So we will still use the Woods-Saxon distribution for $`t_d(b)`$. In heavy nuclear $`A+A`$ collisions, one can incorporate parton energy loss induced by the dense medium through modified fragmentation functions .

The pQCD model described above has been compared to experimental data for $`pp`$, $`p\overline{p}`$ and $`pA`$ collisions at various energies . One found that both the intrinsic $`k_T`$ and the nuclear broadening are very important to describe the existing data in $`pp`$ and $`pA`$ collisions, especially at around SPS energies. One can find some more detailed description of this pQCD-inspired model in Ref. . We will restrict ourselves to the study of the rapidity dependence of the nuclear modification in $`d+A`$ collisions at RHIC and LHC energies in this paper.

Neglecting initial transverse momentum, the initial momentum fractions are related to the transverse momentum and rapidities of the final jets by $`x_1=x_T(e^{y_1}+e^{y_2})/2`$, $`x_2=xT(e^{y_1}+e^{y_2})/2`$, $`x_T=2E_T/\sqrt{s}`$, where $`E_T`$, $`y_1`$ and $`y_2`$ are the transverse momentum and rapidities of the produced jets, respectively. Large positive rapidities, therefore, correspond to large parton fractional momentum $`x_1`$ from the projectile and small momentum fraction $`x_2`$ from the target. Conversely, negative rapidities correspond to small $`x_1`$ and large $`x_2`$. In our parton model calculation, we assume the final hadron rapidity to be the same as that of the fragmenting jets. Shown in Fig. 2 are the rapidity distributions of the charged hadron spectra for four different values of $`p_T`$. Without nuclear modification of the parton distributions (shown as dotted lines), the $`k_T`$ broadening of the projectile partons enhances the particle spectra. The enhancement is the strongest in the forward (projectile) region, thus giving rise to a rapidity asymmetry. The EKS  parameterization has a strong anti-shadowing at $`x_20.10.2`$ for partons from the nuclear target. Such strong anti-shadowing further enhances the rapidity asymmetry (shown as solid lines). However, the HIJING  parameterization has mostly shadowing in this region and thus reduces the rapidity asymmetry in the spectra (dot-dashed line) which is still visible for $`p_T>3`$ GeV/$`c`$. At $`p_T>10`$ GeV/$`c`$, the target parton distribution at $`x_20.20.8`$ is suppressed due to the nuclear binding, which is known as the EMC effect . This reduces the hadron spectra in the backward (target) rapidity region and thus further enhances the rapidity asymmetry.

As one decreases the transverse momentum so that parton shadowing in the target becomes significant, the hadron spectra in the forward region are strongly suppressed. It can even overcome the enhancement caused by transverse momentum broadening. In this case the rapidity asymmetry is reversed, as shown in Fig. 2 for $`p_T=1.5`$ GeV/$`c`$. This is much like the asymmetric rapidity distribution of soft particles produced via soft and coherent interactions as shown in Fig. 1. The validity of the parton model at $`p_T=1.5`$ GeV/$`c`$ may be questionable. However, it clearly shows the trend of the rapidity asymmetry as one changes the value of $`p_T`$.

We also show in Fig. 3 similar rapidity distributions for $`d+Pb`$ collisions at $`\sqrt{s}=5.5`$ TeV where one has access to much larger values of transverse momentum. At such high energies and transverse momenta, the parton scattering cross is much flatter in $`p_T`$ than at lower energies. Therefore, the final hadron spectra are less sensitive to the transverse momentum broadening due to initial multiple scatterings. The rapidity asymmetry caused by the $`k_T`$ broadening is thus very weak, as shown by the dotted lines. The nuclear modification of the parton distributions in the target is causing most of the rapidity asymmetry in the spectra at these energies. Here we used only EKS  parameterization of nuclear modification of the parton distributions. The HIJING  parameterization does not have any scale dependence, which is very important at LHC energies since the shadowing still has significant effects on jet production with large transverse momentum.

To demonstrate the $`p_T`$ dependence of the effect of parton shadowing and $`k_T`$ broadening, we show in Fig. 4 and 5 the nuclear modification factors defined as the ratio of charged hadron spectra in $`d+A`$ over that in $`p+p`$ normalized to the averaged number of binary nucleon collisions,

$$R_{dA}(p_T,y)\frac{d\sigma _{AB}^h/dyd^2p_T}{2Ad\sigma _{pp}^h/dyd^2p_T}.$$
(11)
At the RHIC energy, the transverse momentum broadening enhances the hadron spectra for $`p_T=28`$ GeV/$`c`$ in both forward and backward rapidity regions. The exact enhancement in this $`p_T`$ region depends on the parton shadowing and anti-shadowing. The enhancement in the forward rapidity is larger than in the backward region as demonstrated by the rapidity asymmetry in Fig. 2. For $`p_T>8`$ GeV/$`c`$, the EMC effect of the nuclear modification of parton distribution starts to suppress the hadron spectra. Such a suppression is stronger in the target rapidity region than in the projectile region. We also find that the suppression due to the EMC effect is stronger for kaons than pions. For $`p_T<2`$ GeV/$`c`$, parton shadowing suppresses the hadron spectra in the case of HIJING parameterization. However, the pQCD model might not be valid anymore in this small $`p_T`$ region for a quantitative calculation. The nuclear modification factor for $`d+Pb`$ collisions at the LHC energy shown in Fig. 5 has much smaller variation and the Cronin enhancement is also smaller.

In summary, we have studied the rapidity distribution of hadron spectra in high-energy $`d+A`$ collisions, in particular the rapidity asymmetry caused by multiple parton scatterings. The effects of initial multiple parton scatterings are incorporated via an impact-parameter dependent nuclear shadowing of parton distributions and the broadening of the initial $`k_T`$ carried by partons before they collide and produce high $`p_T`$ hadrons. At low $`p_T`$ parton shadowing suppresses the hadron spectra in the projectile rapidity region, giving rise to a rapidity asymmetry much like the soft particle production via soft and coherent interactions in a string model. However, as one increases $`p_T`$, hadron production is dominated by hard parton scatterings and the rapidity distribution is becoming more symmetric. Within a pQCD model, transverse momentum broadening via initial multiple scattering enhances the hadron spectra in the projectile region at moderate $`p_T=38`$ GeV/$`c`$, causing a rapidity asymmetry opposite to that of soft hadrons. This is in sharp contrast to the prediction of the parton saturation model . The coming data of $`d+Au`$ collisions from RHIC can easily distinguish those two models and will be important to verify whether jet quenching due to parton energy loss is truly the underlying mechanism for the observed hadron suppression in central $`A+A`$ collisions.

We should also caution that our pQCD model cannot take into account the non-perturbative effects that can cause some rapidity asymmetry at very large rapidities. Such non-perturbative effect is responsible for the asymmetric pedestal effect underlying a jet event in $`p+A`$ collisions. Our conclusions are more robust within the central rapidity region where such non-perturbative effect and constraints by total energy conservation are not yet important.

This work was supported by the the Director, Office of Energy Research, Office of High Energy and Nuclear Physics, Division of Nuclear Physics, and by the Office of Basic Energy Science, Division of Nuclear Science, of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.