text:
# Tsunamis in Galaxy Clusters: Heating of Cool Cores by Acoustic Waves

## 1 Introduction

For many years, it was thought that radiative losses via X-ray emission in clusters of galaxies leads to a substantial gas inflow, which was called a “cooling flow” (Fabian, 1994, and references therein). However, X-ray spectra taken with ASCA and XMM-Newton fail to show line emission from ions having intermediate or low temperatures, implying that the cooling rate is at least 5 or 10 times less than that previously assumed (e.g. Ikebe et al., 1997; Makishima et al., 2001; Peterson et al., 2001; Tamura et al., 2001; Kaastra et al., 2001; Matsushita et al., 2002). Chandra observations have also confirmed the small cooling rates (e.g. McNamara et al., 2000; Johnstone et al., 2002; Ettori et al., 2002; Blanton, Sarazin, & McNamara, 2003).

These observations suggest that a gas inflow is prevented by some heat sources that balance the radiative losses. There are two popular ideas about the heating sources. One is energy injection from a central AGN of a cluster (Tucker & Rosner, 1983; Böhringer & Morfill, 1988; Rephaeli, 1987; Binney & Tabor, 1995; Soker et al., 2001; Ciotti & Ostriker, 2001; Böhringer et al., 2002; Churazov et al., 2002; Soker, Blanton, & Sarazin, 2002; Reynolds, Heinz, & Begelman, 2002; Kaiser & Binney, 2003). Recent Chandra observations show that AGNs at cluster centers actually disturb the intracluster medium (ICM) around them (Fabian et al., 2000; McNamara et al., 2000; Blanton et al., 2001; McNamara et al., 2001; Mazzotta et al., 2002; Fujita et al., 2002; Johnstone et al., 2002; Kempner, Sarazin, & Ricker, 2002), although some of them were already discovered by ROSAT (Böhringer et al., 1993; Huang & Sarazin, 1998). Numerical simulations suggest that buoyant bubbles created by the AGNs mix and heat the ambient ICM to some extent (Churazov et al., 2001; Quilis, Bower, & Balogh, 2001; Saxton, Sutherland, & Bicknell, 2001; Brüggen & Kaiser, 2002; Basson & Alexander, 2003). The other possible heat source is thermal conduction from the hot outer layers of clusters (Takahara & Takahara, 1979, 1981; Tucker & Rosner, 1983; Friaca, 1986; Gaetz, 1989; Böhringer & Fabian, 1989; Sparks, 1992; Saito & Shigeyama, 1999; Narayan & Medvedev, 2001).

However, it has already been known that the ICM heating by AGNs or thermal conduction has problems. For the AGN heating, the efficiency of the heating must be quite high (Fabian, Voigt, & Morris, 2002). Moreover, the intermittent activity of an AGN makes the temperature profile of the host cluster irregular, which is inconsistent with observations (Brighenti & Mathews, 2003). For the thermal conduction, stability is the most serious problem; either the observed temperature gradient disappears or the conduction has a negligible effect relative to radiative cooling (Bregman & David, 1988; Brighenti & Mathews, 2003; Soker, 2003). Moreover, thermal conduction alone cannot sufficiently heat the central regions of some clusters (Voigt et al., 2002; Zakamska & Narayan, 2003). Although a “double heating model” that incorporates the effects of simultaneous heating by both the central AGN and thermal conduction may alleviate the stability problem (Ruszkowski & Begelman, 2002), Brighenti & Mathews (2003) indicates that the conductivity must still be about $`0.35\pm 0.10`$ times the Spitzer value.

In this paper, we consider another natural heating source. In the ICM, fluid turbulence is generated by substructure motion or cluster mergers. From numerical simulations, Nagai, Kravtsov, & Kosowsky (2003) showed that the turbulent velocities in the ICM is about 20%–30% of the sound speed even when a cluster is relatively relaxed. Such turbulence generates acoustic waves in the ICM. Compressive characters of the acoustic waves with a relatively large amplitude inevitably lead to the the steepening of the wave fronts to form shocks. As a result, the waves can heat the surrounding ICM through the shock dissipation. A similar heating mechanism has also been proposed in the solar corona; the waves are excited by granule motions of surface convection (Osterbrock, 1961; Ulmschneider, 1971; McWhirter, Thonemann, & Wilson, 1975). The idea of wave heating in the ICM was proposed by Pringle (1989), but the study was limited to order-of-magnitude estimates. In this paper, we study the wave heating by an analytical model and numerical simulations. We use cosmological parameters of $`\mathrm{\Omega }_0=0.3`$, $`\lambda =0.7`$, and $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ unless otherwise mentioned.

## 2 Analytical Approach

### 2.1 Models

In the ICM the magnetic pressure is generally negligible against the gas pressure (Sarazin, 1986). Therefore, acoustic waves (strictly speaking, fast mode waves in high-$`\beta `$ plasma) could carry much larger amount of energy than other modes of magnetohydrodynamical waves. We expect that turbulence in the ICM excites acoustic waves that propagate in various directions. In this paper, we focus on the acoustic waves traveling inward, which play an important role in the heating of the cluster center. These waves, having a relatively large but finite amplitude, eventually form shocks to shape sawtooth waves (N-waves) and directly heat the surrounding ICM by dissipation of their wave energy. We adopt the heating model for the solar corona based on the weak shock theory (Suzuki, 2002; Stein & Schwartz, 1972). In this section, we assume that a cluster is spherically symmetric and stationary. The equation of continuity is

$$\dot{M}=4\pi r^2\rho v,$$
(1)
where $`\dot{M}`$ is the mass accretion rate, $`r`$ is the distance from the cluster center, $`\rho `$ is the gas density, and $`v`$ is the gas velocity. The equation of momentum conservation is

$$v\frac{dv}{dr}=\frac{GM(r)}{r^2}\frac{1}{\rho }\frac{dp}{dr}\frac{1}{\rho c_s\{1+[(\gamma +1)/2]\alpha _w\}}\frac{1}{r^2}\frac{d}{dr}(r^2F_w)$$
(2)
where $`G`$ is the gravitational constant, $`M(r)`$ is the mass within radius $`r`$, $`p`$ is the gas pressure, $`c_s`$ is the sound velocity, $`\gamma (=5/3)`$ is the adiabatic constant, and $`\alpha _w`$ is the wave velocity amplitude normalized by the ambient sound velocity ($`\alpha _w=\delta v_w/c_s`$). For the actual calculations, we ignore the term $`vdv/dr`$ because the velocity is much smaller than the sound velocity except for the very central region of a cluster where the weak shock approximation is not valid ($`\alpha _w1`$; see §2.2). The wave energy flux, $`F_w`$, is given by

$$F_w=\frac{1}{3}\rho c_s^3\alpha _w^2\left(1+\frac{\gamma +1}{2}\alpha _w\right).$$
(3)
Note that the sign of equation (3) is the opposite of equation (7) of Suzuki (2002), because we consider waves propagating inwards contrary to those in Suzuki (2002). The energy equation is written as

$$\rho v\frac{d}{dr}\left(\frac{1}{2}v^2+\frac{\gamma }{\gamma 1}\frac{k_BT}{\mu m_H}\right)+\rho v\frac{GM(r)}{r^2}+\frac{1}{r^2}\frac{d}{dr}[r^2(F_w+F_c)]+n_e^2\mathrm{\Lambda }(T,Z)=0,$$
(4)
where $`k_B`$ is the Boltzmann constant, $`T`$ is the gas temperature, $`\mu (=0.61)`$ is the mean molecular weight, $`m_H`$ is the hydrogen mass, $`n_e`$ is the electron number density, and $`\mathrm{\Lambda }`$ is the cooling function. The term $`𝑭_w`$ indicates the heating by the dissipation of the waves. We adopt the classical form of the conductive flux for ionized gas,

$$F_c=f_c\kappa _0T^{5/2}\frac{dT}{dr}$$
(5)
with $`\kappa _0=5\times 10^7`$ in cgs units. The factor $`f_c`$ is the ratio of actual thermal conductivity to the classical Spitzer conductivity. The cooling function is a function of temperature $`T`$ and metal abundance $`Z`$, and is given by

$`\mathrm{\Lambda }(T,Z)`$ $`=`$ $`2.1\times 10^{27}\left(1+0.1{\displaystyle \frac{Z}{Z_{\mathrm{}}}}\right)\left({\displaystyle \frac{T}{\mathrm{K}}}\right)^{0.5}`$ (7)
$`+\mathrm{\hspace{0.33em}8.0}\times 10^{17}\left(0.04+{\displaystyle \frac{Z}{Z_{\mathrm{}}}}\right)\left({\displaystyle \frac{T}{\mathrm{K}}}\right)^{1.0}`$
in units of $`\mathrm{ergs}\mathrm{cm}^3\mathrm{s}^1`$. This is an empirical formula derived by fitting to the cooling curves calculated by Böhringer & Hensler (1989). We assume that wave injection takes place at radii far distant from the cluster center, and thus there is no source term of waves in equation (4).

The equation for the evolution of shock wave amplitude is given by

$$\frac{d\alpha _w}{dr}=\frac{\alpha _w}{2}\left[\frac{1}{p}\frac{dp}{dr}+\frac{2(\gamma +1)\alpha _w}{c_s\tau }\frac{2}{r}\frac{1}{c_s}\frac{dc_s}{dr}\right],$$
(8)
where $`\tau `$ is the period of waves, which we assume to be constant (Suzuki, 2002). We give the period by $`\tau =\lambda _0/c_{s0}`$, where $`c_{s0}`$ is the sound velocity at the average temperature of a cluster ($`T_{\mathrm{av}}`$), and $`\lambda _0`$ is the wave length given as a parameter. The second term of the right side of equation (8) denotes dissipation at each shock front of the N-waves. We note that the sign of the term is the opposite of equation (6) of Suzuki (2002), because we consider waves propagating inwards contrary to those in Suzuki (2002).

For the mass distribution of a cluster, we adopt the NFW profile (Navarro, Frenk, & White, 1997). The mass profile is written as

$$M(r)\left[\mathrm{ln}\left(1+\frac{r}{r_s}\right)\frac{r}{r_s(1+r/r_s)}\right],$$
(9)
where $`r_s`$ is the characteristic radius of the cluster. The normalization can be given by $`M(r_{\mathrm{vir}})=M_{\mathrm{vir}}`$, where $`r_{\mathrm{vir}}`$ and $`M_{\mathrm{vir}}`$ are the virial radius and mass of a cluster, respectively. We ignore the self-gravity of the ICM.

### 2.2 Results

We show that our model can reproduce observed ICM density and temperature profiles of clusters. We choose A1795 and Ser 159–03 clusters to be compared with our model predictions. Zakamska & Narayan (2003) showed that thermal conduction alone can explain the density and temperature profiles for A1795; $`f_c=0.2`$ is enough and other heat sources are not required. On the other hand, the profiles for Ser $`15903`$ cannot be reproduced by thermal conduction alone (Zakamska & Narayan, 2003). The parameters of the mass profiles for the clusters are the same as those adopted by Zakamska & Narayan (2003) and are shown in Table 1. The concentration parameter of a cluster, $`c=r_{\mathrm{vir}}/r_s`$ is given by

$$c=\frac{1}{r_s}\left[\frac{3M_{\mathrm{vir}}}{4\pi \mathrm{\hspace{0.25em}200}\rho _{\mathrm{crit}}}\right]^{1/3},$$
(10)
where $`\rho _{\mathrm{crit}}`$ is the critical density of the universe. We fix the metal abundance profiles. For A1795, we assume $`Z(r)=0.8\mathrm{exp}(r/170\mathrm{kpc})Z_{\mathrm{}}`$ (Ettori et al., 2002), and for Ser 159–03, we assume $`Z(r)=0.51\mathrm{exp}(r/171\mathrm{kpc})Z_{\mathrm{}}`$ (Kaastra et al., 2001).

We carry out the modeling of ICM heating as follows. First, we select values of $`f_c`$, $`\dot{M}`$, and $`\lambda _0`$. Then, we set the boundary conditions of the equations (1), (2), (4), and (8) at $`r_i=1`$ kpc, that is, well inside the central cD galaxy. From $`r=r_i`$, we integrate the equations outward and compare the model profiles of $`n_e(r_i)`$ and $`T(r_i)`$ with the data. While we fix the value of $`\alpha _w(r_i)`$, we adjust $`n_e(r_i)`$ and $`T(r_i)`$ to be consistent with the observed profiles. We restrict ourselves to a comparison by eye, since neither the data nor the models are reliable enough for a detailed $`\chi ^2`$ fit. If we do not have satisfactory fits, we change the values of $`f_c`$, $`\dot{M}`$, and $`\lambda _0`$ and repeat the process. We show the values of $`f_c`$, $`\dot{M}`$, and $`\lambda _0`$ that we finally adopted in Table 2, and briefly summarize how the results depend on the choice of them as follows.

For A1795, we choose $`f_c=2\times 10^3`$, because Zakamska & Narayan (2003) have already shown that ICM heating only by thermal conduction with $`f_c0.2`$ is consistent with the observations. In this study, we will show that even when $`f_c`$ is much smaller than 0.2, the observed profiles can be reproduced if wave heating is included. However, we found that if $`f_c`$ is too small, the obtained temperature profile is too steep to be consistent with the observation. For Ser 159–03, we adopt $`f_c=0.2`$, which is suggested by Narayan & Medvedev (2001) in a turbulent MHD medium. If we take $`f_c`$ much smaller than this, the model cannot reproduce the relatively flat temperature distribution observed in this cluster.

For mass accretion rates $`\dot{M}`$, we take about 1/10 times the value claimed before the Chandra and XMM-Newton era. For A1795, $`\dot{M}500M_{\mathrm{}}\mathrm{yr}^1`$ ($`h=0.5`$) was reported (Edge, Stewart, & Fabian, 1992; Peres et al., 1998). Thus, we adopt $`\dot{M}=50M_{\mathrm{}}\mathrm{yr}^1`$, which is consistent with a recent XMM-Newton observation ($`<150M_{\mathrm{}}\mathrm{yr}^1`$; Tamura et al., 2001). For Ser 159–03, $`\dot{M}300M_{\mathrm{}}\mathrm{yr}^1`$ ($`h=0.5`$) was reported (White, Jones, & Forman, 1997; Allen & Fabian, 1997). Thus, we adopt $`\dot{M}=30M_{\mathrm{}}\mathrm{yr}^1`$. We note that if we assume that wave heating is effective and that $`\dot{M}`$ is much smaller than the above values, we cannot reproduce both density and temperature profiles obtained by X-ray observations; we get too high temperature and too low density.

Typical wave length, $`\lambda _0`$, should be comparable to the typical eddy size of turbulence in ICM. From numerical simulations, Roettiger, Stone, & Burns (1999) showed that the typical eddy size is the core scale of a cluster. Thus, we take $`\lambda _0=100`$ kpc for A1795. For Ser 159–03, we use a smaller value of $`\lambda _0=70`$ kpc because of its small mass (Table 1). Smaller $`\lambda _0`$ means a smaller distance that waves propagate before dissipation.

Among three of the parameters for the boundary conditions at $`r=r_i`$ ($`\alpha _w`$, $`n_e`$, and $`T`$), we fix $`\alpha _w=3`$ to reduce the number of fitting parameters. If we assume much smaller $`\alpha _w`$, wave heating becomes negligible. On the other hand, if we assume much larger $`\alpha _w`$, the region where the weak shock approximation is invalid ($`\alpha _w1`$) extends.

Figure 1 shows the model fits for the two clusters. The boundary conditions are presented in Table 2. The temperature $`T(r_i)`$ is especially required to be fine-tuned for the fit. We use the Chandra data of A1795 obtained by Ettori et al. (2002) and the XMM-Newton data of Ser 159–03 obtained by Kaastra et al. (2001). The XMM-Newton data of A1795 are also obtained by Tamura et al. (2001) and they are similar to those obtained by Ettori et al. (2002), although the former is not deprojected contrary to the latter. The good agreement between the model and the data suggests that wave heating is a promising candidate of the mechanism that solves the cooling flow problem. In Figure 1, densities go to infinity and temperatures go down to zero at $`r200`$ kpc. This suggests that waves injected outside of this radius cannot reach the cluster center.

In Figure 2a, we present the wave velocity amplitude normalized by the sound velocity ($`\alpha _w`$). As the N-waves propagate into the central regions of the cluster, $`\alpha _w`$ increases rapidly. This is mainly because of the geometrical convergence to the cluster center, whereas the total wave luminosity (= energy flux times $`r^2`$) mostly dissipates through the inward propagation in itself. Since $`\alpha _w>1`$ at $`r4`$ kpc for both A1795 and Ser 159–03, the results may not be quantitatively correct there. We note that the gas velocity for the region of $`r4`$ kpc is very small compared with the sound velocity and thus our ignorance of $`vdv/dr`$ in equation (2) is justified.

In Figure 2b, the ratio $`F_wf_c/F_c`$ is presented. The gaps at $`10`$ kpc reflect $`F_c<0`$. Assume that an observer made observations of the model clusters and the temperature distributions were exactly measured. If the observer assumed the classical conductivity, the heat flux measured by the observer should be $`F_c/f_c`$ ($`f_c<1`$) because of the definition of $`F_c`$ (equation ). Figure 2b shows that the observer would measure an X-ray emission much larger than that predicted by the classical thermal conduction ($`F_wf_c/F_c>1`$) if the energy swallowed by the black hole at the cluster center is small. Such large X-ray emissions have actually been estimated in some clusters (Voigt et al., 2002). Wave heating model can account for the observations without the help of heating by AGNs.

## 3 One-Dimensional Numerical Simulations

### 3.1 Models

In order to be compared with the results in the previous section, we performed one-dimensional numerical simulations. We solve the following equations:

$$\frac{\rho }{t}+\frac{1}{r^2}\frac{}{r}(r^2\rho v)=0,$$
(11)
$$\frac{(\rho v)}{t}+\frac{1}{r^2}\frac{}{r}(r^2\rho v^2)=\rho \frac{GM(r)}{r^2}\frac{p}{r},$$
(12)
$$\frac{e}{t}+\frac{1}{r^2}[r^2v(p+e)]=\frac{1}{r^2}\frac{}{r}\left[r^2\kappa (T)\frac{T}{r}\right]n_e^2\mathrm{\Lambda }(T)\rho u\frac{GM(r)}{r^2},$$
(13)
where the total energy is defined as $`e=p/(\gamma 1)+\rho v^2/2`$, and $`\kappa (T)=\kappa _0T^{5/2}`$. We ignore the self-gravity of ICM. For numerical simulations, we adopt the cooling function based on the detailed calculations by Sutherland & Dopita (1993),

$$n_e^2\mathrm{\Lambda }=[C_1(k_B)^\alpha +C_2(k_BT)^\beta +C_3]n_in_e,$$
(14)
where $`n_i`$ is the ion number density and the units for $`k_BT`$ are keV. For an average metallicity $`Z=0.3Z_{\mathrm{}}`$ the constants in equation (14) are $`\alpha =1.7`$, $`\beta =0.5`$, $`C_1=8.6\times 10^3`$, $`C_2=5.8\times 10^2`$, and $`C_3=6.4\times 10^2`$, and we can approximate $`n_in_e=0.704(\rho /m_H)^2`$. The units of $`\mathrm{\Lambda }`$ are $`10^{22}\mathrm{ergs}\mathrm{cm}^3`$ (Ruszkowski & Begelman, 2002). Note that equation (13) does not include the energy dissipation term contrary to equation (4). This is because the energy dissipation at shocks is included automatically in numerical simulations if shocks are resolved.

The hydrodynamic part of the equations is solved by a second-order advection upstream splitting method (AUSM) based on Liou & Steffen (1993, see also ). We use 500 unequally spaced meshes in the radial coordinate to cover a region with a radius of 300 kpc. The inner boundary is set at $`r=1`$ kpc. The innermost mesh has a width of $`15`$ pc, and the width of the outermost mesh is $`3`$ kpc. The following boundary conditions are adopted:

1. Variables except velocity have zero gradients at the center.
2. The inner edge is assumed to be a perfectly reflecting point.
3. The density and pressure at the outermost mesh are equal to specified values.

Waves are injected at the outermost mesh as

$$v(t)=\alpha _{w0}c_{s0}\mathrm{sin}\left(\frac{2\pi c_{s0}t}{\lambda _0}\right),$$
(15)
where $`\alpha _{w0}`$ is the parameter.

### 3.2 Numerical Results

The mass distribution we assumed is the same as that of A1795 in Table 1 except for $`T_{\mathrm{av}}`$. We assume that the ICM is isothermal and in pressure equilibrium at $`t=0`$. For the NFW profile (equation ), the gas initial density profile is written as

$$\rho (r)=\rho _0\mathrm{exp}[Bf(r/r_s)],$$
(16)
where

$$B=\frac{1}{m(1)}\frac{G\mu m_HM(r_s)}{r_sk_BT_{\mathrm{av}}},$$
(17)
$$m(x)=\mathrm{ln}(1+x)\frac{x}{1+x},$$
(18)
$$f(x)=1\frac{1}{x}\mathrm{ln}(1+x)$$
(19)
(Suto, Sasaki, & Makino, 1998). The initial ICM temperature is $`T_{\mathrm{av}}`$. The density and pressure of the outer boundary are fixed at the initial values. We finish the calculations when the temperature of some of the meshes becomes zero because we do not treat mass dropout from the hot ICM. We define the time as $`t_{\mathrm{cool}}`$. If the temperature does not go to zero, we finish the calculations at $`t=7`$ Gyr. In this section, we set $`n_e(0)=0.017\mathrm{cm}^3`$, $`T_{\mathrm{av}}=7`$ keV, and $`\lambda _0=100`$ kpc. Other model parameters are shown in Table 3.

Figures 3 and 4 show the temperature and density profiles for Models CF and CO, respectively. For reference, the data of A1795 are shown (Ettori et al., 2002), although the detailed comparison is premature. Contrary to Figure 1, we use a linear scale for the distance from the cluster center for both temperature and density profiles to see shock structures. Model CF is a genuine cooling flow model; the temperature goes to zero at the cluster center at $`t_{\mathrm{cool}}=2.7`$ Gyr. In Model CO, conduction dominates cooling and the solution is stable for a long time. The suppression factor of the conductivity, $`f_c=0.2`$, in Model CO is the same as that of Zakamska & Narayan (2003). The temperature gradient at the cluster center is less steep than that obtained by Zakamska & Narayan (2003). This may be because of the differences of the boundary conditions or the initial conditions.

Figures 5 and 6 show the temperature and density profiles for Models W1 and W2, respectively. As we predicted in Figure 2, waves are amplified in the central region. A shock is seen at $`r=13`$ kpc in Model W1. Both temperature and density diverge at the cluster center because these are one-dimensional simulations and waves focus on the cluster center. In the outer region, the shape of waves is complicated. This is because we adopted the reflection condition at the inner boundary and the inward and outward waves interact. For Model W1, $`t_{\mathrm{cool}}=4.7`$ Gyr, and for Model W2, $`t_{\mathrm{cool}}>7`$ Gyr. They are much larger than $`t_{\mathrm{cool}}`$ for Model CF. This means the wave heating is effective.

In Models M1 and M2, we include both thermal conduction and waves. The results of Model M1 and Model M2 are almost the same as those of Model W1 and CO, respectively. The cooling time-scale $`t_{\mathrm{cool}}`$ of Model M1 is shorter than that of Model W1, because shocks are weakened by the thermal conduction. In model M2, the conduction dominates the wave heating.

## 4 Discussion

Our analytical model and numerical simulations show that acoustic waves are amplified at the cluster center and can heat the cluster core. One should note that our assumption, that is, the spherical symmetry, could affect the amplification quantitatively. For more realistic modeling, we should consider that real clusters are not exactly spherically symmetric. However, Pringle (1989) indicated that even if a cluster is not spherically symmetric, the lower temperature and smaller sound velocity at the cluster center should have waves focus on the center. In order to study this effect, we need to perform high-resolution multi-dimensional numerical simulations (Wada, Fujita, &, Suzuki 2003, in preparation). Since the focusing effect depends on the temperature gradient, it may solve the fine-tuning problem of the heating of cluster cores. If the cooling dominates heating, the temperature at the cluster center decreases and the temperature gradient in the core increases. This strengthens the focusing effect and wave heating becomes more efficient. Multi-dimensional simulations have another benefit; it is free from the inner boundary conditions that we set in our one-dimensional simulations.

If waves are actually responsible for the heating in cluster cores, weak shocks should be observed there in some clusters. In our simulations, waves are amplified at $`r1050`$ kpc (Figures 5 and 6). Thus, if the wave length is large ($`100`$ kpc), the shocks are not necessarily observed in all clusters. However, as pointed out by Churazov et al. (2003), the passing of waves may cause gas-sloshing around cluster cores. Observations of fine structures in many cluster cores will be useful to understand the heating mechanism there.

In the central regions of clusters, turbulence generated by interaction of amplified waves from various directions could be observed as broadened metal lines by future X-ray observatories. Observational studies of optical emission lines revealed that there is warm gas at some cluster centers and that the velocity widths range from 100 to $`1000\mathrm{km}\mathrm{s}^1`$ (Hu, Cowie, & Wang, 1985; Johnstone, Fabian, & Nulsen, 1987; Heckman et al., 1989). These warm gas may be embodied in and move with the hot turbulent ICM (see Loewenstein & Fabian, 1990).

## 5 Conclusions

Through analytical and numerical approaches, we have shown that acoustic waves generated by turbulence in the ICM in the outer region of a cluster can effectively heat the central part of the cluster. The heat flux by the waves may exceed that by thermal conduction. The process presented here is phenomenologically analogous to the collapse of a ”tsunami” (a seismic sea wave) at seashore owing to the change of depth of the sea. As in tsunamis, even if the waves have small amplitude at their origin, they could bring huge damage at a distant point, namely the cluster center. Of course, one should note that tsunamis are gravity-driven waves and not acousitc waves. In the analytical studies, we have obtained time-independent solutions and compared the predicted density and temperature profiles with the observed ones; they are consistent with each other. In general, we have confirmed the results obtained by the analytical studies by one-dimensional numerical simulations. Since we assumed that a cluster is spherically symmetric and the assumption leads to artificial focusing of waves, one should take the quantitative results with care. However, it has been indicated that even if a cluster is not spherically symmetric, waves are focused by the temperature gradient at the cluster center. Thus, it is worthwhile to study the wave heating by multi-dimensional analyses.

Y. F. and K. W were supported in part by a Grant-in-Aid from the Ministry of Education, Culture, Sports, Science, and Technology of Japan (Y. F.:14740175; K. W.:15684003). T. K. S. is financially supported by the JSPS Research Fellowship for Young Scientists, grant 4607.


format: text

text:
# Control of Dynamical Localization

## I Introduction

The quantum kicked rotor (KR) and its classical analog, the standard map, have long served as a paradigm for quantum and classical chaos casatibook . Due to its atom optics realization raizenetc , the KR has recently attracted renewed interest. The KR is also of considerable interest in a variety of other fields such as condensed matter physics fishmanprl ; casatiprl , molecular physics fishman ; averbukh , and quantum information science facchi ; georgeot .

One well-known quantum effect in KR is “dynamical localization” (DL) casatibook . That is, although a classical kicked rotor displays unrestricted diffusive energy increase due to classical chaos, only a finite number of rotational states are excited in the quantum dynamics, with the quantum excitation probability versus the rotational quantum number typically assuming a characteristic exponential lineshape. DL is a pure quantum coherence effect and is therefore very sensitive to decoherence. For example, previous studies have demonstrated that noise noiseKR ; cohen ; graham , nonperiodicity in the kicking sequences casati90 ; Abal , and quantum measurements facchi  can destroy DL.

As a coherence effect, DL is also indicative of the possibility of controlling the KR dynamics via quantum effects ricebook ; brumerreview ; brumerbook . Indeed, we recently showed that the quantum phases describing the initial rotor quantum state can be manipulated to effectively control quantum fluctuations in quantum chaos and thus enhance or suppress quantum chaotic diffusion gongcontrol ; gongjcp  in KR. However, manipulating quantum phases in initial states cannot change the unitary evolution operator of the system. Hence, neither the average dynamical localization length nor the characteristic lineshape of dynamical localization can be altered in this way.

Motivated by interest in controlled classically chaotic quantum dynamics gongcontrol ; gongjcp ; gong02 , and to gain more insights into the nature of DL, we consider controlling DL and the associated energy absorption via a modified kicked rotor (MKR) system, in which the phase of the kicking field, or the timing of the kicking sequences (hence the evolution operator), is actively manipulated. In particular, we consider an MKR system in which the sign of the kicking potential is periodically changed, or alternatively time delayed, after a certain number of kicks. As shown below, such a slight modification of KR has profound effects on the dynamics: whereas periodically changing the sign of the kicking potential does not destroy DL, it dramatically changes the quantum diffusion dynamics of KR as well as the nature of DL, offering entirely new opportunities for controlling dynamics, as well as understanding quantum fluctuations and correlations in quantum chaos in periodically kicked systems. For example, we demonstrate that, with the sign of the kicking potential of KR periodically changed: (i) the dynamical localization length is significantly increased so that the energy absorption is strongly enhanced; (ii) classical anomalous diffusion (which can be slower than quantum anomalous diffusion under certain conditions gongpre03 ) can enhance control even when the effective Planck constant is about an order of magnitude larger than the relevant classical phase space structures, and (iii) DL may display strong deviations from purely exponential lineshapes.

This paper is organized as follows: In Sec. II we introduce the modified kicked-rotor model. We then present results in Sec. III on the enhancement of dynamical localization length, with qualitative explanations based upon a known result from the Band Random Matrix theory casatiBRM ; gongBRM ; izrailev . In Sec. IV we show control of DL in a different regime, where the dynamics can be tied to a new mechanism, i.e., the creation of new structures in classical phase space. Section V contains the results on coherent manipulations of the lineshapes for DL. We conclude the paper with a brief discussion in Sec. VI.

## II A Modified Kicked-rotor Model

The KR Hamiltonian is given by

$`H^{KR}(\widehat{L},\theta ,t)=\widehat{L}^2/2I+\lambda \mathrm{cos}(\theta ){\displaystyle \underset{n}{}}\delta (t/Tn),`$ (1)
where $`\widehat{L}`$ is the angular momentum operator, $`\theta `$ is the conjugate angle, $`I`$ is the moment of inertia, $`\lambda `$ is the strength of the kicking field, and $`T`$ is the time interval between kicks. The basis states of their Hilbert spaces are given by $`|m`$, with $`\widehat{L}|m=m\mathrm{}|m`$. The quantum KR map operator for propagating from time $`(N0^+)T`$ to time $`(N+10^+)T`$ is given by

$`\widehat{F}_{KR}=\mathrm{exp}\left[i{\displaystyle \frac{\tau }{2}}{\displaystyle \frac{^2}{\theta ^2}}\right]\mathrm{exp}[ik\mathrm{cos}(\theta )],`$ (2)
with dimensionless parameters $`k=\lambda T/\mathrm{}`$ and the effective Planck constant $`\tau =\mathrm{}T/I`$. For later use we also define the dimensionless scaled rotational energy as $`\stackrel{~}{E}\widehat{L}^2\tau ^2/2\mathrm{}^2`$, where $``$ represents the average over the quantum ensemble. In the $`|m`$ representation, $`\widehat{F}_{KR}`$ takes the following form izrailev :

$`m_1|\widehat{F}_{KR}|m_2=\mathrm{exp}\left(i{\displaystyle \frac{\tau }{2}}m_1^2\right)i^{m_1m_2}J_{m_1m_2}(k),`$ (3)
where $`J_{m_1m_2}(k)`$ is the Bessel function of the first kind of order $`(m_1m_2)`$.

The classical limit of the KR quantum map, i.e., the standard map, depends only on one parameter $`\kappa k\tau `$ and is given by:

$`\stackrel{~}{L}_N`$ $`=`$ $`\stackrel{~}{L}_{N1}+\kappa \mathrm{sin}(\theta _{N1});`$
$`\theta _N`$ $`=`$ $`\theta _{N1}+\stackrel{~}{L}_N,`$ (4)
where $`\stackrel{~}{L}L\tau /\mathrm{}`$ is the scaled c-number angular momentum and $`(\stackrel{~}{L}_N,\theta _N)`$ represents the phase space location of a classical trajectory at $`(N+10^+)T`$. For later discussion we note that for particular values of $`\kappa `$, the classical map Eq. (4) can generate accelerating trajectories whose momentum increases (or decreases) linearly with time (at least on the average). These trajectories are called transporting trajectories schanz . To see this consider the initial conditions: $`(\stackrel{~}{L}=2\pi l_1,\theta =\pm \pi /2)`$ for $`\kappa =2\pi l_2`$, where $`l_1`$ and $`l_2`$ are integers. Clearly, these phase space points are shifted by a constant value ($`\pm 2\pi l_2`$) in $`\stackrel{~}{L}`$ after each iteration, resulting in a quadratic increase of rotational energy. These transporting trajectories are rather stable insofar as they may persist for values of $`\kappa `$ close to $`2\pi l_2`$ (with their average momentum shift after each iteration oscillating around the constant value $`\pm 2\pi l_2`$), thus giving rise to transporting regular islands schanz , i.e., the accelerator modes in the KR case. If classical trajectories are launched from the accelerator modes, they simply jump to other similar islands located in adjacent phase space cells. For trajectories initially outside the accelerator modes, the “stickiness” of the boundary between the accelerator modes and the chaotic sea induces anomalous diffusion over the energy space, i.e., energy increases in a nonlinear fashion, but not quadratically. This is intrinsically different from the case of normal chaotic diffusion in which energy increases linearly with the number of kicks.

We introduce here a slightly modified kicked rotor (MKR) system whose Hamiltonian is given by:

$`H^{MKR}(\widehat{L},\theta ,t)=\widehat{L}^2/2I+\lambda \mathrm{cos}(\theta ){\displaystyle \underset{n}{}}f_M(n)\delta (t/Tn),`$ (5)
where $`f_M(n)`$ is real, $`|f_M(n)|=1`$, and $`f_M(n)`$ changes sign after every $`M`$ kicks. That is, the only difference between KR and MKR is that in the MKR the sign of the kicking potential is changed after every $`M`$ kicks. The effect of changing the sign of the kicking potential can be further understood in terms of the time-evolving wave function, which can be expanded as a superposition of $`|m`$ states: $`_mC_m\theta |m`$, with the expansion coefficients $`C_m`$. Since $`\mathrm{cos}(\theta +\pi )=\mathrm{cos}(\theta )`$, and $`C_m\theta +\pi |m`$=$`_m(1)^mC_m\theta |m`$, changing the sign of the kicking potential is seen to be equivalent to adding a $`\pi `$ phase difference between all neighboring basis sates. Compare this now to the effect of a time delay $`t_d=2\pi T/\tau `$ between two neighboring kicks. Due to the free evolution of the rotor any two angular momentum eigenstates $`|m`$ and $`|m+1`$ will acquire in time $`t_d`$ an additional relative quantum phase given by $`\mathrm{exp}\{i\tau t_d[(m+1)^2m^2]/(2T)\}=\mathrm{exp}[i(2m+1)\pi ]=\mathrm{exp}(i\pi )`$. As such, the MKR can be also realized by introducing the time delay $`t_d`$ after every $`M`$ kicks.

For times $`(N0^+)T`$ to $`(N+M0^+)T`$, the MKR quantum propagator can be written as

$`\widehat{F}_{MKR}=\mathrm{exp}\left(i\pi {\displaystyle \frac{^2}{\theta ^2}}\right)\widehat{F}_{KR}^M\widehat{D}\widehat{F}_{KR}^M.`$ (6)
where this equation defines $`\widehat{D}`$ as the free evolution operator over time $`t_d`$. Here $`\widehat{F}_{KR}^M`$ denotes $`M`$ applications of $`\widehat{F}_{KR}`$. From Eq. (6), one sees that the only difference in time propagation between KR and MKR for every $`M`$ kicks is the $`\widehat{D}`$ operator, whose matrix elements $`m_1|\widehat{D}|m_2`$ are given by

$`m_1|\widehat{D}|m_2=(1)^{m_1}\delta _{m_1m_2}.`$ (7)
Thus we have

$`m_1|\widehat{F}_{MKR}|m_2=(1)^{m_1}m_1|\widehat{F}_{KR}^M|m_2.`$ (8)
Note also that the classical limit of the MKR quantum map \[Eq. (6)\] is given by

$`\stackrel{~}{L}_{N+1}`$ $`=`$ $`\stackrel{~}{L}_N+\kappa f_M(N)\mathrm{sin}(\theta _N),`$
$`\theta _N`$ $`=`$ $`\theta _N+\stackrel{~}{L}_{N+1}.`$ (9)
With respect to model systems in the literature, the MKR here can be regarded as a specific realization of the so-called generalized kicked-rotor model, which was first introduced in Ref. danapre , in the context of quantum anti-resonance. However, it is very different from the amplitude-modulated kicked rotor systems shepelyansky2 ; casati90  previously studied because (i) the kicking field strength here remains constant, and therefore any interesting results arise from pure phase modulations; and (ii) as shown below, the dynamical localization is significantly altered, but not destroyed.

## III Enhanced Dynamical Localization Length

As an example, consider the case of $`k=4.0`$, $`\tau =2.0`$, and $`M=50`$, whose classical limit for both the KR and MKR is fully chaotic and displays normal chaotic diffusion. We demonstrate below that the dynamical localization (DL) and the related energy absorption associated with $`\widehat{F}_{MKR}`$ are dramatically enhanced over that associated with $`\widehat{F}_{KR}^M`$.

Figure 1a displays the angular momentum distribution $`P(m)`$ after $`4\times 10^5`$ kicks, starting with the initial state $`|0`$, for both KR and the MKR with $`M=50`$. The exponential line shape of $`P(m)`$ shown in Fig. 1a indicates that DL occurs in both cases. By fitting $`P(m)`$ with exponentials $`P(m)\mathrm{exp}[|m|/l_{KR}]`$ and $`P(m)\mathrm{exp}[|m|/l_{MKR}]`$ for KR and MKR, respectively, one obtains that the dynamical localization length $`l_{MKR}140.0`$ is significantly larger than $`l_{KR}7.0`$. This clear difference in dynamical localization length is also reflected in the energy absorption shown in Fig. 1b. In particular, while the energy absorption of KR saturates after a few kicks, the MKR system continues to absorb energy in a more or less linear manner for as long as $`10^4`$ kicks. Enhancement of dynamical localization length and energy absorption is also observed for other values of $`M`$, and for a wide range of parameters $`k`$ and $`\tau `$.

As is well known, the DL of KR can be traced back to the localization properties of the eigenstates of the quantum map operator \[Eq. (2)\]. Note first that, due to the rapid decay of $`J_{m_1m_2}(k)`$ with increasing $`|m_1m_2|`$ and the pseudo-random nature of the function $`\mathrm{exp}(i\tau m_1^2/2)`$ in $`m_1`$, the quantum map operator $`\widehat{F}_{KR}`$ in general assumes a band structure and behaves in a pseudo-random manner in the $`|m`$ representation. Hence, below we qualitatively consider the MKR results in terms of a well-known feature from Band-Random-Matrix theory casatiBRM ; gongBRM ; izrailev , namely, that the larger the band width of the quantum map operator the larger the dynamical localization length.

Note first that the matrix $`m|\widehat{F}_{KR}|m^{}`$ is pseudo-random. However, we do not expect $`m|\widehat{F}_{KR}^M|m^{}`$ to be a pseudo-random matrix of the same type since multiplying $`\widehat{F}_{KR}`$ by itself $`M`$ times is expected to establish correlations between the matrix elements $`m_1|\widehat{F}_{KR}^M|m_2`$. Nevertheless, we do assume that the matrix $`m|\widehat{F}_{MKR}|m^{}`$ is banded and pseudo-random since the eigenstates of $`\widehat{F}_{MKR}`$ have no simple connection with those of $`\widehat{F}_{KR}`$. Consider now an arbitrary matrix element $`m_1|\widehat{F}_{KR}^M|m_2`$. Due to the quantum diffusive dynamics within the $`M`$ kicks, the $`m_1|\widehat{F}_{KR}^M|m_2`$ with $`|m_1m_2|>>1`$ should be much greater than $`m_1|\widehat{F}_{KR}|m_2`$ (nevertheless, both of them can be very small). According to Eq. (8), this implies that the matrix element $`m_1|\widehat{F}_{MKR}|m_2`$ is also much greater than $`m_1|\widehat{F}_{KR}|m_2`$. In this sense, we expect that the $`m|\widehat{F}_{MKR}|m^{}`$ matrix should display a wider band than does $`m|\widehat{F}_{KR}|m^{}`$.

The band width of the matrix $`m|\widehat{F}_{KR}|m^{}`$ can be defined by choosing a cut-off value for its matrix elements. One traditional choice is $`|m_1m_2|k`$. In this case, $`m_1|\widehat{F}_{KR}|m_2`$ at the boundary of the band is on the order of $`J_k(k)`$, a number which is sufficiently small. However, this cut-off value, if applied to the matrix $`m|\widehat{F}_{MKR}|m^{}`$, would yield almost the same band width as $`\widehat{F}_{KR}`$. Hence, quantitatively characterizing the band structure of the matrix $`m|\widehat{F}_{MKR}|m^{}`$ is subtle, since the very small matrix elements $`m_1|\widehat{F}_{MKR}|m_2`$ must play an important role in enhancing the dynamical localization length of MKR.

To gain more insight we computationally examined the $`M`$-dependence of the dynamical localization length of the MKR. Specifically, we numerically diagonalized $`\widehat{F}_{MKR}`$, where each matrix is generated using $`2^{14}`$ basis states and is then truncated at dimension $`d`$ chosen below. We characterize the average dynamical localization length by the Shannon entropy $`S_{MKR}`$ sentropy  (note that the Shannon entropy is simply proportional to the dynamical localization length sentropy ) averaged over all approximate eigenstates $`|\varphi _j`$ of $`\widehat{F}_{MKR}`$, i.e.,

$`S_{MKR}={\displaystyle \frac{2}{\alpha d}}{\displaystyle \underset{j=1}{\overset{d}{}}}\mathrm{exp}\left[{\displaystyle \underset{m=d/2}{\overset{d/2}{}}}|m|\varphi _j|^2\mathrm{ln}|m|\varphi _j|^2\right],`$ (10)
where the constant $`\alpha `$ equals 0.96, and $`d`$ is chosen to be 2700. The results for the $`M`$-dependence of $`S_{MKR}`$, for $`M=2`$ to $`M=400`$ is shown in a $`\mathrm{log}\mathrm{log}`$ plot in Fig. 2a. $`S_{MKR}`$ is seen to behave initially as a smooth increasing function of $`M`$, and then to saturate at $`M50`$. To explain the results from the perspective of Band-Random-Matrix theory, we choose a cut-off value for matrix elements $`m_1|\widehat{F}_{MKR}|m_2`$ so as to define the band width $`b_{MKR}`$. Interestingly, we find that this cut-off value must be extremely small (roughly speaking, at least as small as $`10^{10}`$) in order that the $`M`$-dependence of $`b_{MKR}`$ resembles that of $`S_{MKR}`$. For example, Fig. 2b displays $`b_{MKR}`$ as a function of $`M`$ for a cut-off value of $`10^{20}`$. The evident similarities between Fig. 2a and Fig. 2b suggest (i) that we can indeed qualitatively explain the enhanced dynamical localization length in terms of the Band-Random-Matrix theory, and (ii) that even extremely small quantum fluctuations in the values of the matrix elements of the MKR quantum map operator affect its dynamical localization length.

In particular, comparing Fig. 2a with Fig. 2b, one sees that the saturation behavior of $`S_{MKR}`$ for large $`M`$ is consistent with the saturation behavior of $`b_{MKR}`$. The latter reflects the saturation of the matrix elements $`m_1|\widehat{F}_{MKR}|m_2`$ and therefore the matrix elements $`m_1|\widehat{F}_{KR}^M|m_2`$. As such, we infer that the saturation behavior of $`S_{MKR}`$ is simply a result of DL in KR.

Detailed studies on numerous other cases with varying $`k`$ and $`\tau `$ show that the above result, i.e., $`S_{MKR}`$ first increases with $`M`$ and then saturates, is quite general, as long as the system has a completely chaotic classical limit and quantum correlations are insignificant. On the other hand, if $`S_{MKR}`$ behaves differently, then there are two possible origins: either the system is in the deep quantum regime or there are non-negligible regular islands in the classical phase space. For example, in the next section we show cases in which the energy absorption in the MKR with $`M=2`$ is appreciably larger than that in the MKR with $`M=3`$. In these cases one can obtain even more significant changes in DL.

## IV Classical Anomalous Diffusion Assisted Control

In the previous section we studied cases where both KR and MKR essentially have a fully chaotic classical limit. However, as shown below, the MKR can also display new non-chaotic classical phase space structures that are absent in KR. In particular, regular islands with very interesting transporting properties can be induced in MKR. In such cases one can achieve even more dramatic alteration of the DL than that shown above.

Consider first the classical MKR map Eq. (9) for $`M=2`$. Interestingly, in this case there exists transporting trajectories that are different from those in KR. In particular, we have previously observed gongpre03  that for $`\kappa =(2l_2+1)\pi `$, trajectories emanating from $`\stackrel{~}{L}=(2l_1+1)\pi ,\theta =\pm \pi /2`$ will be shifted by a constant value \[$`\pm (2l_2+1)\pi `$\] in $`\stackrel{~}{L}`$ after each kick. This observation suggests that new transporting regular islands that differ from the accelerator modes in KR can be created by changing the sign of the kicking potential after every two kicks. This is confirmed in our extensive numerical studies, both here and in Ref. gongpre03 .

Note that in our previous work gongpre03 , we were most interested in the quantum-classical comparison in anomalous diffusion and considered relatively small effective Planck constants. In that case we found that quantum anomalous diffusion induced by the new transporting regular islands can be much faster than the underlying classical anomalous diffusion. Here, to make a closer connection to atom optics experiments we consider larger $`\tau 1.0`$. In these cases the effective Planck constant is about an order of magnitude larger than the area of the phase space structures associated with classical anomalous diffusion. Intuitively, one would anticipate that such transporting regular islands are too small to be relevant to the quantum dynamics. Surprisingly, this intuition is incorrect, as shown below.

To be more specific, consider first the case of $`\kappa =5.0`$. Figure 3 displays the classical phase space structures of both KR and the MKR with $`M=2`$. While the regular islands seen in Fig. 3a (the KR case) are not transporting (that is, the momentum of the trajectories launched from these islands is bounded and oscillates periodically), a simple computation reveals that the small islands seen in Fig. 3b in the MKR case are transporting. That is, classical trajectories launched from the right (left) transporting regular island shown in Fig. 3b have their momentum shifted by $`\pi `$ ($`\pi `$) on the average after each kick, indicating that these islands originate from the marginally stable point $`\stackrel{~}{L}=(2l_1+1)\pi ,\theta =\pm \pi /2`$ with $`\kappa =\pi `$. Hence, in this case phase manipulation in going from KR to the MKR has both destroyed the nontransporting regular islands of KR and induced new transporting regular islands. Consider a second case with $`\kappa =10.0`$. The corresponding classical phase space structures are shown in Fig. 4a (KR) and Fig. 4b (MKR) (Note that, to clearly display the transporting regular islands, only a part of one phase space cell is shown here). While there are hardly any regular islands seen in Fig. 4a, two small transporting regular islands are seen in Fig. 4b. The average momentum shift for each kick associated with these two islands is found to be $`\pm 3\pi `$, consistent with the fact that $`\kappa =10.0`$ is close to $`3\pi `$.

Consider now the quantum dynamics of these systems. There have been only a few studies on the quantum dynamics of delta-kicked systems where classical chaos coexists with transporting regular islands. Of particular relevance is the previous result that the accelerator modes of KR enhance deviations from the normal DL behavior in KR sundaram ; zaslavsky02 , even for systems far from the semiclassical limit. Since MKR displays new transporting islands, we therefore anticipate that DL may be strongly affected by modifying the Hamiltonian from the KR to MKR system. This is indeed seen below. The results are, however, counter-intuitive, since the classical transporting regular islands created by phase manipulation of the kicking field are found to have an area that is much smaller than the effective Planck constant.

For example, for each of KR and MKR, Figs. 5 and 6 display energy absorption for the cases of $`\tau =1.0`$, $`k=5.0`$ and $`\tau =1.0`$, $`k=10.0`$ (corresponding to $`\kappa =5`$ and $`\kappa =10`$), respectively. Also shown is the MKR case with $`M=3`$, discussed below. It is seen that the energy absorption associated with the MKR with $`M=2`$ (upper dashed line) is much larger than that of KR (solid line). Consider, for example, $`\stackrel{~}{E}`$ at a specific time $`t=1500T`$ when the energy absorption of both KR and MKR has clearly shown signs of saturation (e.g., the average rotational energy may decrease with time due to statistical fluctuations). In the first case (Fig. 5), $`\stackrel{~}{E}=67.7`$ for KR and $`\stackrel{~}{E}=1269.7`$ for MKR. In the second case (Fig. 6), $`\stackrel{~}{E}=798.0`$ for KR and $`\stackrel{~}{E}=10049.0`$ for MKR. In both cases a control factor larger than an order of magnitude has been achieved in going from the KR Hamiltonian to the MKR case with $`M=2`$.

This is not the case for $`M=3`$, $`\tau =1.0`$, $`k=5.0`$ shown in Fig. 5. Here the energy absorption in the MKR with $`M=3`$ is only slightly larger than in the KR and far less than in the MKR with $`M=2`$. Similarly, for the case of $`\tau =1.0`$, $`k=10.0`$ shown in Fig. 6, although energy absorption in the MKR with $`M=3`$ is much enhanced (compared with KR), it is still not as significant as in the MKR with $`M=2`$. This indicates, as confirmed by directly examining $`P(m)`$ (not shown) after saturation, that for both cases the dynamical localization length of the MKR with $`M=3`$ is no larger than that of the MKR with $`M=2`$, contrary to what is observed in the previous section. This is because the underlying classical dynamics of MKR here is not completely chaotic, i.e., the classical dynamics displays characteristics of anomalous diffusion due to transporting regular islands in phase space. Thus one can expect statistical deviations from the Band Random Matrix theory, used previously to relate $`M`$ to the extent of control.

These results emphasize that the control mechanism here is uniquely based upon the transporting regular islands created by our control scenario. This is further supported by the lineshape for DL, which can be strongly nonexponential, as discussed below.

## V Nonexponential Lineshapes for Dynamical Localization

It was pointed out more than two decades ago fishman  that the DL of KR can be mapped onto the problem of Anderson localization in disordered systems. In particular, an exactly soluble case of disorder in tight-binding models, i.e., the Lloyd model Lloyd , suggests that dynamical localization should assume an exponential line shape, at least on the average. This has been confirmed by numerous computational studies on KR. For example, Fig. 1a clearly demonstrates, for both KR and MKR, that the distribution function $`P(m)`$ can be fit beautifully with an exponential function with a characteristic localization length.

However, the Hamiltonian nature of the KR and MKR implies that there always exist some subtle quantum phase correlations in the quantum dynamics. Hence, in addition to some universal properties of DL, the DL lineshape can display rich non-universal properties, e.g., the system may display nonexponential dynamical localization. Nonexponential dynamical localization has been previously observed in KR but its origins are still poorly understood satija ; sundaram .

Here, we demonstrate that the MKR with $`M=2`$ can display strongly nonexponential line shapes for DL, rarely seen in KR. We focus on the MKR with $`M=2`$ since the classical MKR with $`M=2`$ has transporting regular islands that are absent in KR. Second, transporting regular islands may induce large fluctuations in DL sundaram . However, we examine below cases with connections to anomalous diffusion as well as those without clear connections to anomalous diffusion. We have also studied other versions of MKR with $`M2`$, and have found that nonexponential lineshapes for DL in the latter case are much less common than in the $`M=2`$ case.

Figure 7 compares nonexponential lineshapes for DL in the MKR (upper dashed lines) to the analogous exponential lineshapes for DL in KR (solid lines), for four different values of $`k`$ and $`\tau `$. The lineshapes are obtained by propagating the quantum dynamics for 8000 kicks from the initial state $`|0`$, and will remain essentially the same for longer propagation times note . The huge difference between the KR and MKR lineshapes is striking. As shown in Fig. 7, in the MKR case, $`P(m)`$ plotted in the $`\mathrm{log}`$ scale displays structures that are far from a purely exponential lineshape. For example, one sees that the initial exponential decay rate of $`P(m)`$ with $`|m|`$ is considerably smaller than its large-$`|m|`$ exponential decay rate, suggesting that multiple characteristic lengths are needed to describe the DL of the MKR. Also seen is that the difference in $`P(m)`$ between KR and MKR can be as large as ten orders of magnitude or more. Further, Figs. 7b and 7c indicate that the decay rate of $`P(m)`$ for large $`|m|`$ in the MKR case is very similar to the KR case, although this is not the case in Figs. 7a and 7d.

To the best of our knowledge, this is the first demonstration that intriguing differences in the lineshapes for DL can be created by simply changing the sign of the kicking potential periodically.

The four cases of nonexponential lineshapes shown in Fig. 7 can be divided into two categories, based upon the properties of their underlying classical dynamics. The classical dynamics associated with the cases in Figs. 7a and 7b, shown in Figs. 3b and 4b, displays transporting regular islands. The presence of these classical structures implies an inhomogeneous classical phase space, and, as demonstrated in the previous section, may have a significant impact on the quantum dynamics even when their size is much smaller than the effective Planck constant. In this regard, the nonexponential lineshapes shown in Figs. 7a and 7b may not be totally surprising. However, for the other two MKR cases shown in Figs. 7c and 7d, we did not find any regular islands in their classical phase space even when examined on a very fine scale, suggesting that their classical dynamics is essentially fully chaotic. Thus, understanding the nonexponential lineshapes shown in Fig. 7c and Fig. 7d will be even more challenging.

The nonexponential lineshapes for DL arise from extremely small quantum fluctuations and residual quantum correlations in quantum chaos. To be able to resolve the nonexponential lineshape for DL, $`P(m)`$ has to be known with high precision. For example, the two shoulders shown in Fig. 7b involve occupation probabilities $`P(m)`$ as small as $`10^{14}`$. It is therefore not surprising that, while it is common to have strongly nonexponential lineshapes for DL in the MKR case, each individual lineshape is highly sensitive to the exact value of the effective Planck constant. For example, for the case shown in Fig. 7b, increasing the value of $`\tau `$ from $`2.0`$ to $`2.0+10^5`$ can completely destroy the nonexponential lineshape! This drastic change in the lineshape for DL even occurs without causing an obvious difference in energy absorption behavior. Evidently then, both experimental observations and theoretical predictions of nonexponential dynamical localization are far from trivial and are in need of further study.

## VI Discussion and Conclusions

This paper has dealt with control of dynamical localization in kicked rotor systems. In all cases we manipulate the external kicking field to alter the properties of the rotor system, i.e. the distribution of population amongst rotor energy levels after saturation as well as the energy absorption. In particular, we have examined the effect of introducing a reversal of the kicking field after $`M`$ kicks which, within the framework of quantum mechanics, corresponds to introducing a phase shift amongst rotor energy levels.

Two parameter regimes have been examined, one which shows enhanced DL lengths with increasing $`M`$, and the other which need not. This behavioral difference can be understood in terms of the character of the underlying classical phase space: the former systems are completely chaotic whereas the latter show a mixed phase space that includes transporting regular islands. Indeed, we have found that even if the transporting islands are tiny compared to the effective Planck constant, they still have a profound effect on the control of the DL. Further, a comparison of the traditional kicked rotor system to the modified kicked rotor systems shows that the latter is much more capable of displaying nonexponential dynamical localization. Thus, by modifying the kicking potential are able to control the dynamical localization in the kicked rotor.

The results of this paper are relevant to two fields of study: quantum control and kicked rotor dynamics. From the control perspective, modifying the kicking field changes the dynamics. However, this system does not obviously permit a picture in terms of interfering quantum pathways (the standard view of weak field coherent control brumerreview ; brumerbook ) since (a) the kicking field is always on, and (b) there are a multitude of interfering transitions responsible for the observed behavior. Indeed, it is even difficult to isolate the interfering pathways that are responsible for dynamical localization in the simple kicked rotor, whose dynamics is easier than that of the MKR.

From the viewpoint of kicked rotor studies in the field of quantum chaos, this paper provides insights into the quantum dynamics in the case displaying classical anomalous diffusion. As one of the results of this study, we find that classical transporting regular islands can dramatically affect the quantum dynamics even when their size is much smaller than the effective Planck constant. Further, the MKR system proposed in this paper provides a new model for the study of quantum dynamics where the underlying classical chaos coexists with transporting regular islands. By choosing proper system parameters, we can create transporting regular islands whose size varies from being much larger to being much smaller than that of the accelerator modes of KR. Encouraged by this, we plan in the near future to further use the MKR to study quantum tunneling between the transporting regular island and the chaotic sea gongpre03 ; zaslavsky02  and between transporting regular islands, and to study the phase space structure of quantum eigenstates schanz .

In the fully chaotic case, our approach suggests the need for additional studies of dynamical localization from the Band-Random-Matrix theory perspective. For example, we have qualitatively explained the results in Sec. III in terms of well known features of Band-Random-Matrix theory. However, we found that the band width of the quantum map operator had to be defined using an extremely small cut-off value for the matrix elements, suggesting that a quantitative understanding of the MKR results such as Fig. 2 may require new models of Band-Random-Matrix ensembles.

The strongly nonexponential lineshapes for DL found in the MKR with $`M=2`$ further demonstrate the need for more theoretical work on properties of DL. In particular, our results should motivate greater interest in characterizing and understanding nonexponential dynamical localization, with efforts directed towards explaining why nonexponential dynamical localization occurs for some system parameters and not for others. This is of importance in understanding the high sensitivity of nonexponential lineshapes for DL to the exact value of the effective Planck constant.

We have chosen the system parameters to be within the reach of current atom optics experiments on KR raizenetc . Although experimental studies of nonexponential lineshapes for DL are difficult, we believe that it is straightforward to experimentally observe the results of Sec. III and Sec. IV. Apart from the atom optics realization of KR and MKR, it is also interesting to consider a molecular version of KR and MKR, i.e., diatomics periodically kicked in strong microwave fields fishman ; gongjcp . Preliminary computational studies Hans  confirm that directly observing quantum control of dynamical localization in molecular rotational motion is possible, e.g., in the case of quantum anomalous diffusion. Another promising experimental realization of KR and MKR requires kicked particles in a square-well potential laksh01 . Along this direction an interesting model (which is very different from ours) has recently been proposed for the study of classical and quantum anomalous diffusion laksh02 .

In summary, consideration of control in classically chaotic quantum systems is of general interest and importance to both fields of quantum chaos and quantum control. In this paper we have demonstrated, via a modified kicked rotor model, that dynamical localization, perhaps the best known phenomenon in quantum chaos, can be modified over a wide range. The results are of both experimental and theoretical interest.

This work was supported by the U.S. Office of Naval Research and the Natural Sciences and Engineering Research Council of Canada. HJW was partially supported by the Studien-stiftung des Deutschen Volkes and the Barth Fonds of ETH Zürich.


text:
# The Rest–Frame UV Luminosity Density of Star–Forming Galaxies at Redshifts 𝑧>3.5Based on observations obtained with the NASA/ESA Hubble Space Telescope obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. (AURA) under NASA contract NAS 5-26555.

## 1 INTRODUCTION

The amount of star formation that took place early in the cosmic history, e.g. at $`z>2`$, and its evolution with redshift are key pieces of information to constrain theories of galaxy evolution. Evidence of the possible signature of the trailing edge of the cosmic reionization epoch at redshifts $`z\stackrel{>}{}6`$ (e.g. Fan et al. 2002) has also recently renewed interest in measuring the star formation activity at these epochs, since starburst galaxies can contribute a significant fraction of the ionizing photons (Madau et al. 1999; Haiman et al. 2001; Steidel et al. 2001). Initial estimates based on the evolution of the cosmic neutral hydrogen mass density estimated from damped Ly$`\alpha `$ absorption systems (DLA, Fall et al. 1996) suggested that star formation activity peaked in the past at around $`z1`$ and then decreased at higher redshifts. The first direct measures based on the observed UV radiation of relatively unobscured star–forming galaxies in the local and intermediate–redshift universe and of Lyman–break galaxies (LBGs) at $`z>2`$ (Madau et al. 1996) seemed to confirm this general scenario. However, subsequent measures based on the systematic spectroscopic identification of hundreds of Lyman–break galaxies at $`z3`$ and $`z4`$ showed that the star formation density traced by these sources, after reaching a maximum value somewhere between $`z1`$ and $`2`$, remains nearly constant up to the highest redshift observed with some confidence, i.e. $`z4`$ (Steidel et al. 1999, S99 hereafter). More recent measure at $`z5`$ from deep, wide–area samples obtained with the Subaru telescope (Iwata et al. 2003) and with HST (Bouwens et al. 2003) show very mild evolution from $`z3`$, although possible evidence to the contrary has also been reported (Stanway et al. 2003a,b). Measures based on photometric redshifts also seem to be in general agreement with a constant or mildly decreasing activity of star formation as traced by starburst galaxies up to $`z6`$ (Thompson et al. 2001; Thompson 2003; Fontana et al. 2003; Kashikawa et al. 2003), or even suggest that the star–formation activity might actually increase with redshift (Lanzetta et al. 2002).

The current measures of UV specific luminosity density of LBGs come from the Hubble Deep Field survey (HDF, Williams et al. 1996, 2000) and other deep, but relatively small survey made with HST (Bouwens et al. 2003) or from ground–based surveys. The HST surveys are limited by its very small area, and thus cosmic volume, which makes it prone to systematic bias induced by cosmic variance (for example, S99 argue that the original result by Madau et al. (1996) was, at least in part, the effect of cosmic fluctuations). The ground–based surveys cover a much larger area, but at high redshifts the brightness of the night sky limits the sensitivity to only the bright end of the luminosity function. Unfortunately, LBGs are rather faint, with a characteristic magnitude $`z_{850}^{}=26.0`$ at $`z=6`$ (based on Adelberger & Steidel 2000’s UV luminosity function at $`z=3`$) if no luminosity evolution occurs between $`z3`$ and $`z6`$. Thus, one needs to reach rather faint flux limits, even for 8–m telescopes, to account for most of the light from LBGs once the redshift exceeds $`z4`$.

The GOODS/ACS observations offer a good compromise, covering an area $`33`$ times larger than the combined HDFs, with much deeper sensitivity out to wavelengths of nearly 1 $`\mu `$m than that of most ground–based samples. This makes it possible to detect galaxies with luminosity fainter than $`L_3^{}`$ up to $`z6.5`$. In this letter we present measures of the comoving specific luminosity density over the range $`3.5\stackrel{<}{}z\stackrel{<}{}6.5`$ based on the stack of three epochs of observations in both GOODS fields, the HDF–N and CDF–S. A companion letter by Dickinson et al. (2003, D03 hereafter) in this volume analyzes in greater detail the selection of galaxies at $`z6`$, including a discussion of the first spectroscopic identifications. Throughout, magnitudes are in the AB scale (Oke 1977), and the world model, when needed, is that of a flat universe with density parameters $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and Hubble constant $`H_0=70`$ km s-1 Mpc-1.

## 2 DATA AND SAMPLE SELECTION

The data used in this letter are comprised of the first three epochs of observations, out of a total of five, of the GOODS/ACS survey (Giavalisco et al. 2003). They consist of eight mosaics, four in each of the two survey fields, and one in each of the $`B_{435}`$, $`V_{606}`$, $`i_{775}`$ and $`z_{850}`$ bands. Each mosaic covers an area approximately 10’x16’ in size with a scale of 0.05 arcsec/pixel. Table 1 of Giavalisco et al. (2003) lists the relevant parameters of the data, including the sensitivity. Source catalogs have been extracted using the SExtractor software, and the procedures are also described in Giavalisco et al. (2003). The catalogs used to select high–redshift galaxies have been created performing the source detection in the $`z_{850}`$ band and then using the isophotes defined during this process as fixed apertures for photometry in the other bands. We derive colors from isophotal magnitudes, but use the AUTO aperture–corrected magnitudes (in the SExtractor parlance) to measure the $`z_{850}`$–band apparent magnitudes.

Samples of star–forming galaxies at high redshifts have been extracted from the matched catalogs using the Lyman–break technique (Madau et al. 1996; Dickinson 1998; S99, Steidel et al. 2003; see also Giavalisco 2002). Specifically, we have defined LBGs at $`z4`$ ($`B_{435}`$–band dropouts) using the colors equations:

$$(B_{450}V_{606})1.2+1.4\times (V_{606}z_{850})(B_{450}V_{606})1.2(V_{606}z_{850})1.2$$
$`(1.1),`$
and LBGs at $`z4`$ ($`V_{606}`$–band dropouts) by:

$$[(V_{606}i_{775})>1.5+0.9\times (i_{775}z_{850})]$$

$$[(V_{606}i_{775})>2.0](V_{606}i_{775})1.2(i_{775}z_{850})1.3$$
$`(1.2),`$
where $``$ and $``$ are the logical OR and AND operators, respectively. For bands which are entirely shortward of the Lyman discontinuity (e.g. the $`B_{435}`$ for the $`V_{606}`$– and $`i_{775}`$–band dropouts and the $`V_{606}`$ band for the $`i_{775}`$–band dropouts) we also request a non detection with $`S/N<2`$.

These selection criteria are largely based on our previous experience with Lyman–break galaxies (e.g. Dickinson 1998; S99; Steidel et al. 2003). We have visually fine–tuned them based on the observed colors of stars and galaxies in the ACS images, as well as on the ACS synthetic photometry from galaxy spectral templates to reject most interlopers from lower redshifts, while efficiently detecting typical UV–bright, star–forming galaxies at redshifts of interest.

Without a third band in the near–infrared, it is not possible to use a 2–color selection to define samples of LBGs at $`z6`$ in analogy to what done above. Instead, we have used the single–color threshold

$$(iz)1.3,$$
$`(1.3)`$
along with a non–detection requirement $`S/N<2`$ in the $`B_{435}`$ and $`V_{606}`$ bands to define our sample of $`z6`$ $`i_{775}`$–band dropouts. D03 discuss the selection and robustness of these dropouts in detail, including the use of available near–IR photometry as “a posteriori test” of the selection criteria. A number of galaxies meeting these criteria have already been confirmed spectroscopically to be at $`z6`$ (Bunker et al. 2003; D03; Stanway et al. 2003b). Later on we discuss the implications of this different selection of LBGs on the measure of the specific luminosity density of the $`i_{775}`$–band dropouts.

We have included in the samples only galaxies with $`S/N5`$ in the $`z_{850}`$ band, and we have visually inspected each of them and removed sources that were deemed as not real, either artifacts or spurious detections (estimated using counts of negative sources detected in the same data set). These amount to a negligible number for the $`B_{435}`$ and $`V_{606}`$ samples, and $`12`$% for the $`i_{775}`$ one. We have also eliminated from all samples sources with stellar morphology down to apparent magnitude $`z_{850}`$$`26`$, i.e. where such a morphological classification is reliable. These amount to the $`3.1`$%, $`8.3`$% and $`4.6`$% for the $`B_{435}`$, $`V_{606}`$ and $`i_{775}`$ samples, respectively. While the procedure bias our samples against LBGs that are unresolved by the ACS, it prevents contamination by galactic stars. In practice, it results in negligible changes to the measured specific luminosity density. Furthermore, while the vast majority of galaxies at $`z6`$ have $`i_{775}`$-$`i_{775}`$$`<1.3`$, photometric errors scatter some of these galaxies into our selection window. Thus, we have also purged the $`i_{775}`$–band sample from photometric–scatter interlopers using a statistical procedure. We have used a bright sub–sample of field galaxies, for which the photometry is very accurate, as a template of the color distribution of the sample as a whole and estimate the contamination using photometric errors as function of magnitude derived from the Monte Carlo simulations in each band (see D03).

Down to $`z_{850}`$$`26.5`$, roughly the $`50`$%–completeness flux limit for unresolved sources (see Figure 4 of Giavalisco et al. 2003), the culled samples include 1115, 275 and 122 $`B_{435}`$–band, $`V_{606}`$–band and $`i_{775}`$–band dropouts, respectively. With a survey area of 316 arcmin2, this corresponds to surface density $`\mathrm{\Sigma }=3.50\pm 0.10`$, $`0.87\pm 0.05`$ and $`0.39\pm 0.03`$ galaxies per arcmin2, respectively, where the error bars have been derived assuming Poisson fluctuations.

## 3 THE UV COMOVING SPECIFIC LUMINOSITY DENSITY

Because no systematic spectroscopic observations of the GOODS are available yet, we have used Monte Carlo simulations to estimate the redshift distribution function of the three samples. The technique consists of generating artificial LBGs distributed over a large redshift range (we used $`2.5z8`$) with assumed distribution functions of UV luminosity (we used a flat distribution, see later), SED, morphology and size. We have adjusted the input distribution functions of SED and size by requiring that the distribution functions recovered from the simulations using identical procedures as for the real Lyman–break galaxies match the observation at $`z4`$. In this way, both simulations and observations are subject to similar incompleteness, photometric errors in flux and color, blending, and other measurement errors. We have obtained the SED using the unobscured synthetic spectrum of a continuously star–forming galaxy with age $`10^8`$ yr, Salpeter IMF and solar metallicity (Bruzual & Charlot 2000), and reddened it with the starburst extinction law (Calzetti 2000) and $`E(BV)`$ randomly extracted from a gaussian distribution with $`\mu _{E(BV)}=0.15`$ and $`\sigma _{E(BV)}=0.15`$. For the cosmic opacity, we have used the Madau (1995) prescription extrapolated to $`z=8`$. We have used an equal number of $`r^{1/4}`$ and exponential profiles with random orientation, and size extracted from a log–normal distribution function, as described by Ferguson et al. 2003, this volume, Figure 2).

For a given dropout sample, the main output from the simulations is the probability function $`p(M,z,m)`$ that a LBG with absolute magnitude $`M`$ at redshift $`z`$ is observed as having apparent magnitude $`m`$. A commonly used technique to derive the specific luminosity density of LBGs in this case is that of the “effective volume” (S99), where the spatial volume occupied by galaxies with apparent magnitude $`m`$ in the sample is

$$V_{eff}(m)=p(M,z,m)𝑑M\frac{dV(z)}{dz}𝑑z.$$
$`(2)`$
The comoving specific luminosity density contributed by such galaxies is then estimated as

$$d(m)=\frac{n(m)L(m,\overline{z})}{V_{eff}(m)}dm,$$
$`(3)`$
where $`\overline{z}`$ is the average redshift of the simulated galaxies that have been selected into the sample, $`n(m)`$ are the number counts of real Lyman–break galaxies observed with magnitude $`m`$, and $`L(m,\overline{z})`$ is the specific luminosity of one such galaxy if placed at redshift $`\overline{z}`$. The specific luminosity density is then $`=𝑑(m)`$. The method provides a relatively accurate estimate in the case of the color–color selected $`B_{435}`$–band and $`V_{606}`$–band dropout samples, regardless of the assumptions about the UV luminosity distribution function used in the simulations, because these galaxies are selected based only on their observed SED, i.e. redshift, with little dependence of their intrinsic luminosity. In other words, these LBGs have a relatively tight correlation of $`m`$ and $`M`$ up to the the detection limit.

For the one–color selected $`i_{775}`$–band dropouts the $`V_{eff}`$ method underestimates the value of $``$, because in this case there is no tight correlation between absolute and apparent magnitude in the sample. Galaxies with a given apparent magnitude $`m`$ now have absolute magnitude $`M`$ distributed in a much larger interval, because whether or not they enter the sample depends on both their redshift and on their absolute magnitude. Specifically, the lower redshift bound is set by the color threshold, while the upper bound depends strongly on galaxy luminosity, as IGM opacity suppresses the $`z_{850}`$–band flux at higher redshifts. Hence, the effective volume $`V_{eff}(m)`$ is overestimated (unless in the simulation one uses the same intrinsic luminosity function fo the real galaxies), because in the simulations galaxies with widely different absolute magnitude are equally represented. In this case we have measured $``$ with a different procedure. We have done a $`\chi ^2`$ minimization to find the intrinsic luminosity function $`\varphi (M)`$ of the simulated galaxies such that once inserted into the real images and retrieved their number count

$$n_s(m)=p(M,z,m)\varphi (M)𝑑M𝑑z$$
$`(4)`$
is best–fit to the number count of the real galaxies. We use a Schechter function with slope fixed to the value found at $`z3`$ (Adelberger & Steidel 2000, AS00 hereafter), namely $`\alpha =1.6`$, and derive the parameter $`M^{}`$ and $`\varphi ^{}`$ from the fit. In practice, since our data only reach down to $`L0.2L_3^{}`$ at most, the assumption of a fixed value of $`\alpha `$ is a relatively minor one (we verified that values in the range $`1.4\alpha 1.8`$ change our results by at most $`12`$%). The specific luminosity density is then simply $`=L\varphi (M)𝑑M`$, and for each redshift value this is computed using the corresponding sample’s best–fit luminosity function (a detailed description of the fits and of the analysis of the result is the subject of a forthcoming paper).

We found the average redshift and standard deviation of the redshift distribution to be $`z_B=3.78`$ and $`𝒮_B=0.34`$ for the $`B_{435}`$–dropout sample, $`z_V=4.92`$ and $`𝒮_V=0.33`$ for the $`V_{606}`$ sample, and $`z_i=5.74`$ and $`𝒮_i=0.36`$ for the $`i_{775}`$ one. For each sample we have estimated $``$ at $`1500`$ Å ($`_{1500}`$) using K–corrections derived from a template LBG with the same UV color as the average one in the $`B_{435}`$–band dropout sample. We have derived $`_{1500}`$ integrating the luminosity function down to the “observed” limit, i.e. the absolute magnitude that for each sample corresponds to an apparent magnitude of $`z_{850}`$=26.5, and we have also estimated the (relatively small) corrections so that all the measures reach the same absolute magnitude, which we have chosen to be that of the B–band dropout sample ($`0.2L_3^{}`$). Finally, when using the $`V_{eff}`$ method, we have also estimated the magnitude–dependent corrections for light losses due to finite aperture photometry that we have derived from the Monte Carlo simulations.

The top panel of Figure 1 shows $`_{1500}`$ for the GOODS three samples derived using the two methods above, both the “observed” values and those corrected to the same absolute magnitude, compared to the measures at $`z3`$ and $`4`$ by S99. It can be seen that the two methods provide similar results for the $`B_{435}`$ and $`V_{606}`$ dropout samples, but differ for the $`i_{775}`$–band dropouts. Note that the luminosity function correction is relatively small for the $`V_{606}`$–band dropouts, but increases for the $`i_{775}`$–band ones, as the $`z_{850}`$=26.5 limit of the sample corresponds to a brighter intrinsic luminosity at higher redshift. The bottom panel of the figure shows the star–formation density as a function of redshift. For the GOODS points we have used the conversion factor by Madau et al. (1998), namely SFR$`=1.4\times 10^{28}`$ in units of M yr-1 Mpc-3; the other points are from S99, after conversion to our world model and correcting for the smaller range of absolute luminosity, i.e. $`>0.2L_3^{}`$ instead of $`>0.1L_3^{}`$. Both the observed values and those corrected for dust obscuration as suggested by AS00 are plotted.

## 4 DISCUSSION AND CONCLUSIONS

The specific luminosity density of LBGs appears to depend rather weakly on redshift over the range $`2.5\stackrel{<}{}z\stackrel{<}{}6.5`$. Integrating down to $`L0.2L_3^{}`$, S99 report $`_{1500}(z=3)=1.50\pm 0.10`$, while from the GOODS data we find $`_{1500}(z=4)=1.63\pm 0.05`$, $`_{1500}(z=5)=1.04\pm 0.08`$, and $`_{1500}(z=6)=1.15_{0.19}^{+0.24}`$ in units of $`\times 10^{26}`$ erg s-1 Hz-1 Mpc-3 (1–$`\sigma `$ error bars). However, while the points at $`z5`$ are relatively robust, the constraint at $`z6`$ is still somewhat weak, because it critically depends on the measure of the number counts near the sensitivity limit of the survey (see D03). For example, if we restrict the sample to $`z_{850}`$–band photometry with $`S/N6.5`$ (roughly $`z_{850}`$$`<26`$) then we found that the best-fit specific luminosity density drops down to $`_{1500}(z=6)=0.42_{0.30}^{+0.33}\times 10^{26}`$ erg s-1 Hz-1 Mpc-3 ($`0.42_{0.40}^{+0.60}`$ 2–$`\sigma `$). At face value this is still a relatively mild drop from the $`z3`$ value, a factor $`3.5`$, smaller than the factor $`7`$ proposed by Stanway et al. (2003), but the error is too large for the constraint to be meaningful. We note, though, that the $`z5`$ point, which is more robust than the $`z6`$ one and does not need as big an extrapolation down the luminosity function, is reasonably consistent with the mild evolution scenario, lending support to our best estimate at $`z6`$.

Thus, the result at $`z6`$ depends somewhat critically on the apparently substantial amount of light contributed by the faint galaxies. To the best of our knowledge the galaxies with $`5S/N6.5`$ are real $`z6`$ candidates. Based on our color selection, the observed dispersion of colors of faint galaxies, our simulations, and the expected colors of galaxies of various spectral types, we believe that our statistical corrections to the contamination is adequate and is not introducing a major systematic error in our measures. The very good agreement with the measure by S99 at $`z4`$, which is supported by systematic redshift identification of the galaxies, adds credence to this assertion. The current spectroscopic identifications at $`z6`$ (Bunker et al. 2003; D03; Stanway et al. 2003b) agree well with our predicted redshift distribution and also support this conclusion. The samples are still too sparse, however, to attempt a measure of the efficiency of the selection criteria and, thus, the contamination. Clearly, the measure of $`_{1500}(z=6)`$ clearly needs to be revisited with deeper data.

In any case, our measure at $`z6`$ is in overall good agreement with other similar measures from HST and ground–based data (Bouwens et al. 2003; Lehnert & Bremer 2003). Stanway et al. (2003a,b) report a factor $`7`$ decrease of the specific luminosity density, although this is also very likely in agreement with our resultof because they limit their measure to bright galaxies, i.e. $`L>L_3^{}`$. Actually, a direct quantitative comparison is difficult because of the difference in the sample selection and the pronounced dependence of $`_{1500}`$ on the faint counts. These authors also use the GOODS data, but they base their source detection on single–epoch images (we use a stack of three epochs), therefore necessarily reaching a shallower flux level and larger incompleteness111They also use a slightly different color selection criterion, namely $`i_{775}`$-$`z_{850}`$$`>1.5`$. As D03 observe, this misses one of the galaxies spectroscopically confirmed at $`z=5.8`$.. Down to $`z_{850}`$$`<25.6`$ they find a total of 14 galaxies over 350 arcmin2 (0.040 arcmin-2), while to the same flux level and using their same criteria we find 35 galaxies over 316 arcmin2 (0.11 arcmin-2). Also, note that they use the $`V_{eff}`$ method to derive $``$, which, as we have detailed above, underestimates the specific luminosity density for the $`i_{775}`$–band dropouts.

Our measures are in good quantitative agreement with another deep, large–area LBG survey at $`z5`$ (Iwata et al. 2003) and also with surveys based on photometric redshifts (Thompson et al. 2001; Thompson 2003; Fontana et al. 2003; Kashikawa et al. 2003). They differ from some of the photometric redshift estimates (from the HDF) by Lanzetta et al. (2002), who, after correcting for bias in the photometry due to the $`(1+z)^4`$ cosmological dimming, conclude that the specific star–formation density either remains constant up to $`z6`$, and then increases afterward, or increases monotonically. While we do explicitly account in our measure corrections for light losses due to finite aperture photometry (estimated from the Monte Carlo simulations), we do not find as large corrections as those suggested by Lanzetta et al (2002). Very likely, much of the difference is that our corrections are derived by requiring that the galaxy size distribution used as input in our Monte Carlo simulations is such that the output distribution matches the observations at $`z4`$ (this roughly corresponds to assuming that the sizes evolve as $`H(z)^1`$; see Ferguson et al., this volume), while Lanzetta et al.’s corrections are based on the light distribution of galaxies at $`1\stackrel{<}{}z\stackrel{<}{}1.5`$, i.e. at considerably lower redshifts.

It is not known if the UV specific luminosity density of Lyman–break galaxies at $`z>2.5`$ is a fair tracer of the cosmic star formation activity (e.g. AS00). Regardless, if the dust obscuration properties of these sources are similar to local starburst galaxies (Meurer et al. 1999; Calzetti 2000), and do not significantly change over the redshift range $`3<z<6.5`$, then our result can be rephrased by saying that the star–formation activity of LBGs decreases very mildly with increasing redshift, at $`z6`$ being $`25`$% lower than it was around its maximum at $`1\stackrel{<}{}z\stackrel{<}{}3`$, several Gyr earlier. Since LBGs are certainly sites of formation of a considerable amount of the present–day cosmic stellar mass density (AS00; Giavalisco 2002) it appears that the onset of substantial cosmic star formation takes place earlier than $`z>6`$, namely at less than $`7`$% of the cosmic age. Finally, Figure 1 also shows that there is an overall good agreement between the dust–corrected data points and the semi–analytical models by Somerville et al. (2001, blue curve), and between the uncorrected data points and the predictions based on the observed evolution of the neutral HI mass density as traced by DLA (which include the effect of dust obscuration) at least up $`z4`$, where these systems are known with goods statistics (Pei et al. 1999, red curves).

What is the number of hydrogen–ionizing photons expected from our $`i_{775}`$band dropouts at $`z6`$? Since the rest-frame UV continuum at 1500 Å is dominated by the same short-lived, massive stars that are responsible for the emission of photons shortward of the Lyman edge, the needed conversion factor, about one ionizing photon above 1 ryd for every 5 photons at 1500 Å, is fairly insensitive to the assumed initial mass function and is independent of the galaxy history for ages $`10^7`$ yr (Madau, Haardt, & Rees 1999; Haiman, Abel, & Madau 2001). Here we normalize the number of ionizing photons to the observed 1500 Å flux (rest-frame), i.e. we bypass the need for any correction due to dust extinction. A comoving luminosity density of $`_{1500}=10^{26}`$erg s-1 Hz-1 Mpc-3 implies then a comoving production rate of H-ionizing photons $`\dot{n}_{\mathrm{ion}}3\times 10^{51}`$ s-1 Mpc-3, or about 5$`f_{\mathrm{esc}}`$ photons per hydrogen atom per $`5\times 10^8`$ yr escaping into the intergalactic medium. Here $`f_{\mathrm{esc}}`$ is the average escape fraction of ionizing radiation from the galaxy HIlayers relative to the escape fraction at 1500 Å. Photoionization of intergalactic hydrogen requires more than one photon above 1 ryd per hydrogen atom, as extra photons are needed to keep the gas in overdense regions and filaments ionized against radiative recombinations. If $`f_{\mathrm{esc}}`$ is greater than a few tens of a percent, then our population of $`i_{775}`$ dropouts may contribute significantly to the UV metagalactic flux, and help to reionize the universe by redshift 6.

Support for the GOODS HST Treasury program was provided by NASA through grants HST-GO09425.01-A and HST-GO-09583.01 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, under NASA contract NAS5-26555. PM acknowledges support by NASA through grant NAG5-11513. We thank the referee, Dr. Adriano Fontana, for a very careful and thoughtful report.