text:
# Constraints to the evolution of Ly-𝛼 bright galaxies between z=3 and z=6Based on observations obtained at the ESO VLT, Paranal, Chile; ESO programs 67.A-0175 and 68.B-0088

## 1 Introduction

Important observational information about the epoch of galaxy formation can only be gained from the existence, number counts and properties of galaxies at very high redshift. In the past couple of years substantial progress has been made in discovering galaxies at redshifts $`z>4`$. Several dozens of Lyman break galaxies have been detected at $`z4`$ (e.g., Steidel et al. steidel99 (1999)). In addition, a substantial number of galaxies at $`z>5`$ have been detected and verified spectroscopically (e.g., Dey et al. dey (1998); Weyman et al. weyman (1998); van Breughel et al. vanBr (1998); Hu et al. hu99 (1999); Hu et al. hu02 (2002); Ellis et al. ellis (2001); Dawson et al. dawson (2002); Ajiki et al. ajiki (2002); and Rhoads et al. rhoads02 (2002)); one of these (van Breughel et al. vanBr (1998)) is a radio galaxy at $`z=5.2`$. Quasars have been identified out to $`z=6.3`$ (e.g., Pentericci et al. pent (2002); Fan et al. fan01 (2001)). However, there is a lack of *systematic* surveys for high redshift primeval galaxies. Many of the galaxies at $`z>5`$ have been detected serendipitously.

Moreover, as it has been pointed out, e.g., by Steidel et al. (steidel96 (1996)), the strong metal absorption lines in the spectra of Lyman break galaxies indicate that they have already formed at least one generation of stars which has chemically enriched the systems and lead to dust formation, making Ly-$`\alpha `$ emission weak or undetectable. Therefore, the Lyman break galaxies do not show the properties expected for primeval galaxies. In contrast to the Lyman break galaxies, all but one of the galaxies found at $`z>5`$ exhibit a very strong Ly-$`\alpha `$ emission line, in agreement with population synthesis models for primeval galaxies of Charlot & Fall (charlot (1993)), which predict that a young, dust free, star forming galaxy should show strong Ly-$`\alpha `$ emission with (intrinsic) equivalent widths in the range of 5-25 nm.

A direct probe of galaxy formation is to determine the number counts and redshift evolution of those nascent galaxies. Several systematic surveys aimed at discovering a significant number of $`z>4.5`$ galaxies have now been proposed and are being carried out, e.g., Cowie & Hu (cowie (1998)), Rhoads et al. (rhoads (2000)) and our survey, the Calar Alto Deep Imaging Survey (CADIS, Meisenheimer et al. meise98 (1998), meise03 (2003)). At $`z=4.5`$, a universe with standard parameters $`H_0=70`$ km s-1 Mpc-1, $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ is only 1.2 Gyr old, corresponding to a look-back time of $`>90`$% of the age of the universe. Any object observed at this or an earlier epoch should mark the early stages of galaxy formation.

Unlike the Lyman break technique (e.g., Steidel & Hamilton steidel92 (1992), steidel93 (1993)), which selects galaxies by the presence of young stars and intervening HI absorption, CADIS selects galaxies by the extreme equivalent widths of their Lyman-$`\alpha `$-emission ($`W_{\lambda _{rest}}>10\mathrm{nm}`$). It has been demonstrated by Kudritzki et al. (kudr (2000)) that the high Ly-$`\alpha `$ equivalent width of these galaxies can only arise in regions nearly free of dust, and therefore the Ly-$`\alpha `$ emission line is a good tracer of star formation from primordial material. Unlike similar programmes at 8$``$10 m telescopes, CADIS aims to find *the most luminous* starbursts, with SFR $`10\mathrm{M}_{}`$/year, which most likely led to the formation of the galactic bulge. Even in the age of 10-m telescopes, only these luminous starbursts are accessible for detailed spectroscopic follow-up.

This paper is structured in the following way: In section 2 we describe how we select Ly-$`\alpha `$ galaxy candidates at $`z>4.7`$ with CADIS. In section 3 we present the spectroscopic follow-up of some Ly-$`\alpha `$ candidates. Finally, in section 4 we discuss the comparison between observed and theoretical number counts of Ly-$`\alpha `$-emitting galaxies up to $`z6`$.

## 2 Lyman-$`\alpha `$ candidates from the CADIS survey

CADIS combines a moderately deep multi-band survey (10 $`\sigma `$ limit, R$`{}_{\mathrm{lim}}{}^{}=24`$) with a deep emission line survey employing an imaging Fabry-Perot-Interferometer ($`\mathrm{F}_{\mathrm{lim}}=310^{20}\mathrm{Wm}^2`$), which covers three waveband windows essentially free from atmospheric OH emission lines: A ($`\lambda =702\pm 6`$ nm), B ($`\lambda =819\pm 5`$ nm), and C ($`\lambda =920\pm 8`$ nm). The classification and analysis of emission line galaxies from CADIS are outlined in Hippelein et al. (hippe (2003)), Meisenheimer et al. (meise03 (2003)). Therefore, in this section we will present only those aspects which are relevant to the search of Ly-$`\alpha `$ galaxies using the CADIS (emission line) survey.

The optimal detection of a weak emission line galaxy superimposed on the bright night sky is reached when the instrumental resolution $`\delta \lambda `$ is adapted to the expected line width. Assuming $`\mathrm{\Delta }v300\mathrm{km}\mathrm{s}^1`$, a typical virial velocity of a Milky Way-like galaxy, yields $`\delta \lambda (1+z)(\mathrm{\Delta }v/c)\lambda _00.9`$ nm for a Lyman-$`\alpha `$ emission line at $`z6`$. On the other hand, the detection probability increases with the observed volume $`\mathrm{\Delta }V\mathrm{\Delta }z\mathrm{\Delta }\mathrm{\Omega }\mathrm{\Delta }\lambda \mathrm{\Delta }\mathrm{\Omega }`$, where $`\mathrm{\Delta }\mathrm{\Omega }`$ denotes the solid angle surveyed. Using a Fabry-Perot (FP) Etalon with $`\mathrm{\Delta }\lambda 2`$ nm, which is tuned in 9 steps and covers an interval of $`1215`$ nm, CADIS provides an optimal compromise of both sensitivity and search volume. The three wavelength intervals (A, B, and C, see Figure 1) which are scanned by CADIS correspond to Lyman-$`\alpha `$ redshifts of $`z4.8`$, 5.7, and 6.6, respectively. The field size is $`\mathrm{\Delta }\mathrm{\Omega }100`$ arcmin2. Model predictions (Thommes & Meisenheimer 1995) of the galaxy abundance down to a line flux limit F$`{}_{\mathrm{line}}{}^{}3.010^{20}\mathrm{Wm}^2`$ (which can be reached with the Calar Alto 2.2- and 3.5-m telescopes within a few hours of observing time) made clear, before the project had even started, that at most a handful of Lyman-$`\alpha `$ galaxies would be expected per field and wavelength interval in the case of an open universe, which is commonly assumed today ($`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$). Accordingly we decided to search every wavelength interval in five fields. After six years of observations on Calar Alto, 90% of the planned Fabry-Perot data have been obtained. The Fabry-Perot wavelength intervals which have been almost completely analysed are given in Table 1. The analysis of interval C has been postponed for the moment because of technical problems (strong “fringing” in the CCD images), and low expectations of success on the basis of the present results at $`z<6`$, which will be presented in this paper.

The following results are based on the data in the FP-windows in four CADIS fields presented in Table 1. It should be noted that the analysis of 69% of window B, and *only* 27% of window A has been finished (see Table 1), in terms of the final $`\mathrm{\Delta }\mathrm{\Omega }\times \mathrm{\Delta }z`$: $`\mathrm{\Delta }\mathrm{\Omega }_{final}`$ will be $`500\mathrm{arcmin}^2`$, corresponding to five CADIS fields, $`\mathrm{\Delta }z_{final}`$ for each Fabry-Perot window will be $`\mathrm{\Delta }z_{final}0.1`$.

The current classification of emission line galaxies shows that more than 97% of the emission line galaxies are foreground galaxies ($`z<1.2`$), in which we mainly detect the lines $`\mathrm{H}\alpha `$, $`\mathrm{H}\beta `$, $`[\mathrm{O}\mathrm{III}]\lambda \mathrm{\hspace{0.17em}5007}`$, and $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$. Therefore, the main challenge is to separate those foreground objects from the Ly-$`\alpha `$ candidates. We identify Lyman-$`\alpha `$-galaxy candidates among the emission line galaxies by the following exclusion criteria:

(1) There should be no flux below the Lyman limit of a Lyman-$`\alpha `$ galaxy, because of an intrinsic drop in the spectra of hot O and B stars at the Ly-limit (Charlot & Fall charlot (1993)), and absorption by the neutral interstellar medium either in the galaxy itself (Leitherer et al. 1995a , 1995b ) or in Lyman-limit foreground systems close to the primeval galaxy. For $`z>4.7`$ the Lyman limit lies at $`\lambda >520`$ nm. Therefore, no flux should be detected in the CADIS B filter ($`\lambda _{cent}=461`$ nm, FWHM$`=113`$ nm). Accordingly, we require that $`F_B<2\widehat{\sigma }_B`$ (“blue dropouts”), where $`F_B`$ is the flux measured in the B filter, and $`\widehat{\sigma }_B`$ its true overall error which takes into account both calibration errors and the scatter of the flux between the individual images. From 614 emission line galaxies selected in four fields in the FP-windows given in Table 1 this criterion yields $`\mathrm{N}_{\mathrm{Noblue}}=85`$ galaxies, i.e. 13.8% of the 614 emission line galaxies.

(2) Every single remaining Ly-$`\alpha `$ candidate has to be carefully checked on the FP, pre-filter and R images because of the possible contamination by nearby bright objects. Candidates which are closer than about 3′′ to a bright object are most likely spurious, since the changing atmospheric conditions change the tails of bright objects affecting the photometry of close-by galaxies. Therefore, such spurious objects are removed from the list of candidates. After this step, 5.7% of the emission line galaxies remain.

(3) Foreground galaxies, for which we detect one of the prominent emission lines $`\mathrm{H}\alpha `$, $`\mathrm{H}\beta `$, $`[\mathrm{O}\mathrm{III}]\lambda \mathrm{\hspace{0.17em}5007}`$, or $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ in the Fabry-Perot, can be excluded as soon a secondary line is detected in one of the veto filters, leaving 4.7% of the galaxies.

(4) The distinction between a (rare) galaxy with bright $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ emission but weak restframe UV continuum, i.e., undetectable in the CADIS B filter, and a Lyman-$`\alpha `$-galaxy is ambiguous based on the CADIS veto filters, because no line bluewards of $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ can be detected. Instead, we have to use the entire spectrum for this decision. Lyman-$`\alpha `$ galaxies at high redshift should distinguish themselves by a continuum step across the Lyman-$`\alpha `$ line, due to absorption from neutral hydrogen in the Ly-$`\alpha `$ clouds and Ly-limit systems along the line of sight (Madau madau (1995)). Therefore, the candidates which remain after step (3) are classified in two categories: likely Ly-$`\alpha `$, if they show almost no continuum on the blue side of the Ly-$`\alpha `$ emission line, e.g., no significant flux in the R filter, and likely $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ objects otherwise. However, we need spectroscopic follow up observations of candidates belonging to *both* categories, since it cannot be excluded for sure, that some Ly-$`\alpha `$ galaxies fall between the likely $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ objects. After this step, 16 Lyman-$`\alpha `$-candidates (classified as most likely Ly-$`\alpha `$) remain (see Table 2), i.e., 2.6% of the emission line galaxies we found in four CADIS fields. In addition we are left with 13 likely $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ candidates.

## 3 Spectroscopic observations

Due to the low abundance of Ly-$`\alpha `$ galaxies (see section 4), spectroscopic follow-up at large telescopes constitutes an important part of the survey: First, we have to search close to the detection limit, which makes contamination by noise nonnegilible, and second, even very rare and unlikely contaminants, like distant supernovae or other transient objects, reach surface densities comparable to that of the Ly-$`\alpha `$ candidates, and can be mistaken as emission line galaxies in our observations of the FPI scans which are spread over years in some fields.

The first goal of the spectroscopic follow-up of likely Ly-$`\alpha `$ galaxies with 8-m class telescopes is therefore to verify the emission line. Furthermore, the line shape and the continuum blue- and red-wards of the emission line could allow us to decide between Ly-$`\alpha `$ and $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$: If the line has been verified, but the resolution is not good enough to identify the emission line, a second step is the clear confirmation of the line using higher resolution grisms, in order to see if the line shows the asymetric profile expected for Ly-$`\alpha `$, or the double \[O II\]$`\lambda \lambda `$3726,3729 line with separation $`\mathrm{\Delta }\lambda (1+z)`$, where $`\mathrm{\Delta }\lambda =2.75`$Å . If galaxies with bright Ly-$`\alpha `$ line at $`z>4.7`$ are confirmed, we plan a detailed study of the corresponding galaxy as third step.

Thus, confirming Ly-$`\alpha `$ candidates by spectroscopy is a slow process. Nevertheless, the verification or non-verification of a line for every single Ly-$`\alpha `$ candidate is very important in order to set more robust upper limits to the number counts of galaxies at high redshift, which can set stringent constraints on theoretical models (see section 4).

Spectroscopic follow-up observations of Ly-$`\alpha `$ candidates in the 01h- and 23h-fields were obtained in the summer and autumn of 2001, using FORS 2 at the VLT. The slitmasks contained 1.0′′ wide slits of length between 10′′ and 20′′. The position angle of each mask was chosen such that the number of Ly-$`\alpha `$ candidates in the 6.8$`\times `$ 6.8 FORS field of view is maximized. Regions of the slitmasks not devoted to primary targets were used to obtain spectra of low metallicity emission line galaxies. The results of the observations of these low metallicity emission line galaxies will be presented and discussed in another paper (Maier et al. maier03 (2003), in preparation).

Depending on the distribution of candidates, we used two different grisms, the lower resolution 300 I grism, which gives a spectral resolution of about 1.2 nm, at 800 nm, for slitlets 1′′ wide, and the 600 RI grism, which gives a higher spectral resolution of about 0.8 nm, at 800 nm, for slitlets 1′′ wide. Three $`z5.7`$ Ly-$`\alpha `$ candidates, 01h-3238, 01h-28090, and 01h-27927, were observed for 215 min using the 300 I grism with FORS2, and one $`z4.8`$ Ly-$`\alpha `$ candidate, 23h-50707, was observed for a total of 150 min using the 600 RI grism with FORS2. Fluxes were calibrated using multiple observations of the spectrophotometric standard stars LTT 7379, EG 274 and LTT 7987 (Hamuy et al. hamuy92 (1992),hamuy94 (1994)). Seeing varied between 0′′.5 and 1′′.2.

For two of the four observed Lyman-$`\alpha `$ candidates we verified the emission line, for 01h-3238 at $`819.0\pm 0.3`$ nm, and for 23h-50707 at $`705.7\pm 0.3`$ nm (see Figure 2 and Figure 3). No continuum and no additional emission lines are seen on the VLT spectra, and, according to the CADIS measurements, these two objects satisfy the criteria (1)-(4) for Ly-$`\alpha `$ candidates. Therefore, we conclude that the two emission lines are very likely Ly-$`\alpha `$ lines of high redshift galaxies, 01h-3238 at $`z=5.735\pm 0.003`$ and 23h-50707 at $`z=4.803\pm 0.003`$, respectively.

The line of the other two observed Ly-$`\alpha `$ candidates, 01h-28090 and 01h-27927, is not seen on the VLT spectra. These are the Ly-$`\alpha `$ candidates with the faintest fluxes of our Ly-$`\alpha `$ candidates list. One explanation for the non-detection could be that the slit possibly missed (part of) the line emitting region of these two galaxies. This is indicated by some emission line objects on the same masks as the Ly-$`\alpha `$ candidates which show a spectroscopic determined flux smaller than the CADIS flux. Thus, the slit width of 1′′.0 used for the VLT observations was possibly too narrow, and we may have thus measured only a fraction of the emission line flux of the Ly-$`\alpha `$ candidates. The reason for this is that our astrometry delivers rather accurate coordinates (better than 0′′.2) for objects with multiple detections (FP and continuum), but can be off by more than 0′′.5 for objects detected only in the Fabry-Perot (Ly-$`\alpha `$ candidates). On the other hand, the two Ly-$`\alpha `$ candidates could be indeed spurious, since we expect a residual contamination of about 50% in our list of Ly-$`\alpha `$ candidates. The reason for the non-verification of the lines has to be established by extending the statistics of spectroscopic follow-up observations.

It should be noted that we could verify an emission line (at $`820.9\pm 0.3`$ nm) also for the galaxy 01h-4616. This object, which shows no flux in the B filter and in the veto filters, passed the criteria (1)-(3), but was classified as a probable $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ galaxy at $`z1.2`$ according to criterion (4). No other lines are seen on the VLT spectrum, and a continuum on both sides of the line is detected. Therefore, the emission line is indeed likely to be $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ at $`z=1.202\pm 0.001`$.

We calculate the absolute Ly-$`\alpha `$ luminosity of 01h-3238 and 23h-50707 from the observed fluxes measured in the FP scan (since the spectroscopic observations are subject to slit losses, see above). The Ly$`\alpha `$ flux is $`f`$(Ly$`\alpha `$) = $`(4.1\pm 0.8)\times 10^{20}`$ W m-2 for 01h-3238, and $`f`$(Ly$`\alpha `$) = $`(4.0\pm 0.8)\times 10^{20}`$ W m-2 for 23h-50707. Using a $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and H$`{}_{0}{}^{}=70`$ km s-1 Mpc-1 cosmology, we obtain an absolute Ly$`\alpha `$ luminosity of $`L`$(Ly$`\alpha `$) $`(14.5\pm 3.0)\times 10^{35}`$ W for 01h-3238, and $`L`$(Ly$`\alpha `$) $`(9.5\pm 1.9)\times 10^{35}`$ W for 23h-50707. Note that these Ly$`\alpha `$ luminosities are the highest found among any $`z5`$ galaxies (luminosities are estimated by using the same cosmology): e.g., $`3.3\times 10^{35}`$ W for the lensed galaxy HCM 6A at $`z=6.56`$ (Hu et al. 2002); $`6.1\times 10^{35}`$ W for SSA22-HCM1 at $`z=5.74`$ (Hu et al. 1999); $`3.4\times 10^{35}`$ W for HDF 4-473.0 at $`z=5.60`$ (Weymann et al. 1998). Taking the Ly-$`\alpha `$ to $`\mathrm{H}\alpha `$ ratio for Case B recombination and no dust (Brocklehurst brockl (1971)), together with the Kennicutt (kenn (1983)) conversion between $`\mathrm{H}\alpha `$ luminosity and star formation rate, the derived star formation rates for the two Ly-$`\alpha `$ galaxies are: $`(14.5\pm 3.0)M_{}\mathrm{yr}^1`$ for 01h-3238, and $`(9.5\pm 1.9)M_{}\mathrm{yr}^1`$ for 23h-50707. Kudritzki et al. (kudr (2000)) found that the derived star formation rates depend very strongly on the assumption whether we witness a burst or continuous star formation, the SFR being higher in the case of a starburst. Therefore, the derived SFRs using Kennicutt (kenn (1983)) empirical formula for the case of continuous star formation represent only lower limits for the SFRs. The true values of the SFRs in case of a starburst may be a factor of five higher.

## 4 Comparison between observed and theoretical abundances of Ly-$`\alpha `$ emitting primeval galaxies

The CADIS selection suppresses the contamination in the list of Ly-$`\alpha `$ candidates to less than 50% (see section 2). As a consequence, the list of CADIS Ly-$`\alpha `$ candidates allows us to set *stringent* upper limits for the density of Ly-$`\alpha `$-emitting galaxies at $`z>4.7`$. We can therefore discuss the evolution of the population of Ly-$`\alpha `$ emitting galaxies from $`z=3`$ to $`z=6`$ by comparing the maximum abundance derived from the list of CADIS Ly-$`\alpha `$ candidates with the results of other systematic surveys at $`3<z<6`$.

Figure 4 shows the cumulative number counts of Lyman-$`\alpha `$-galaxies at $`z=4.8`$ (left panel) and at $`z=5.7`$ (right panel), presented as the total number of galaxies, $`N`$ per deg2 per $`\mathrm{\Delta }z=0.1`$, which are brighter than a certain observed flux, $`F_{lim}`$. In the cumulative histogram of CADIS Ly-$`\alpha `$ galaxy candidates, the upper edge represents the total number of CADIS candidates per deg 2 and per $`\mathrm{\Delta }z=0.1`$, brighter than the respective $`F_{lim}`$, based on the list of Ly-$`\alpha `$ candidates from Table 2: e.g., at $`z=5.7`$, in four fields with a total $`\mathrm{\Delta }\mathrm{\Omega }_{total}=413`$ arcmin$`{}_{}{}^{2}1/9\mathrm{}^{}`$, we find one candidate brighter than $`F_{lim}=6.410^{20}\mathrm{W}\mathrm{m}^2`$, two candidates brighter than $`F_{lim}=5.710^{20}\mathrm{W}\mathrm{m}^2`$, etc. The filled squares additionally take into account the likely residual contamination by noise plus unidentified $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$-emitting galaxies, and therefore represent our best-guess upper limit to the abundance of Ly-$`\alpha `$ galaxies.

Additional abundance measurements, showed as filled triangles, are from Hu et al. (hu98 (1998)) at $`z=4.5`$, and from Hu et al. (hu99 (1999)) at $`z=5.7`$. Hu et al. detected two Ly-$`\alpha `$ galaxies at $`z=4.52`$ in a 24 arcmin2 field, and one Ly-$`\alpha `$ galaxy at $`z=5.74`$ using narrowband observations in a 30 arcmin2 field. Stars indicate preliminary results from the Large Area Lyman Alpha Survey (LALA) at $`z=4.6`$ (Malhotra & Rhoads malh (2002)) and at $`z=5.7`$ (Rhoads et al. rhoads02 (2002)). Malhotra & Rhoads (malh (2002)) found 157 Ly-$`\alpha `$*candidates* at $`4.37<z<4.57`$ in one field of 0.36 deg2. The arrow in Fig. 4 (left panel) indicates the upper limit of the number density of Ly-$`\alpha `$-emitters based on the list of these candidates. Rhoads et al. (rhoads02 (2002)) reported 18 Ly-$`\alpha `$*candidates* at $`z5.7`$ in one field of $`0.2`$ deg2, of which three (out of four observed) candidates have been verified spectroscopically. Rhoads et al. infer a number density $`39<N<54`$ of Ly-$`\alpha `$-emitters per deg2 per $`\mathrm{\Delta }z=0.1`$ down to their detection limit (see Fig. 4, right panel). For comparison, we show in Fig. 4 the measured abundance of Lyman-$`\alpha `$-emitting galaxies at $`z=3.5`$ through emission-line-surveys from Hu et al. (hu98 (1998), small circles), and Kudritzki et al. (kudr (2000), large circles; converted from $`z=3.1`$).

In order to parametrize the observed number counts at $`z=3.5`$, we use the model by Thommes & Meisenheimer (thommes03 (2003)): It basically assumes (i) that Ly-$`\alpha `$-galaxies mark the onset of star-formation in present-day spheroids, (ii) a rather short-lived Ly-$`\alpha `$-bright phase due to rapid dust formation ($`10^8`$ yrs), and (iii) a mass dependent formation history as given by the peak formalism. In order to fit the abundance at $`z=3.5`$ we used model parameters which set the peak of Ly-$`\alpha `$-emission for the bulge of the Milky Way around $`z6`$. The normalization at $`z=3.5`$ is set by the abundance of present day spheroids together with a free parameter $``$ the duration of the Ly-$`\alpha `$-bright phase. For the present work, we decided not to fine tune the model in a way that a consistent description of the abundances at all redshifts $`3.5z5.7`$ is obtained. This will be presented in a forthcoming paper. Here we rather would like to demonstrate a robust qualitative statement $``$ namely that Ly-$`\alpha `$-galaxies are less common at $`z=5.7`$ than at $`z=3.5`$. To this end we simply *shift* the model function for $`z=3.5`$ to higher redshift by taking into account both the higher luminosity distance (shift to the right), and the smaller comoving volume $`\mathrm{\Delta }V/\mathrm{\Delta }z`$ (shift to bottom). This “no-evolution model” is shown by the dotted line in Fig. 4.

It is obvious from Fig. 4 that the CADIS upper limits both at $`z=4.8`$ and $`z=5.7`$ fall short of such a non-evolution model, while at the faint end, F$`{}_{\mathrm{line}}{}^{}10^{20}\mathrm{Wm}^2`$, the present results could still be compatible with no evolution.

Since our statistics at $`z4.8`$ is extremely limited, we will focus any further discussion on redshifts around $`z=5.7`$.

While the number counts of Hu et al. at $`z=5.7`$ (a single galaxy found by Hu et al. hu99 (1999)) still seem compatible with the no-evolution model, the abundance of Ly-$`\alpha `$-galaxies from LALA and the limits given by the CADIS candidates fall significantly short of the no-evolution prediction. The under-abundance of Ly-$`\alpha `$ emitting galaxies at $`z=5.7`$ compared to $`z=3.5`$ becomes even more obvious if one considers that the noise distribution of our survey predicts residual contamination of about 50%. Nevertheless, before drawing definite conclusions, one should consider whether selection effects could lead to an under-estimation of the number of Ly-$`\alpha `$-emitting galaxies. Potential influences on the upper limits include: objects classified as probable $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ galaxies for which the emission line is actually Ly-$`\alpha `$, the fainter surface brightnesses of high redshift objects, and the large-scale structure. The weakly confined bright end of the luminosity function at $`z=3.5`$ may limit the comparison at the bright end. Therefore, we have to discuss how robust the CADIS upper limits are.

The distinction between a galaxy with $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ emission, but no flux in the CADIS B filter, and a Lyman-$`\alpha `$-galaxy is ambiguous (see step (4) in Section 2). Therefore, some of the probable $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$ objects may turn out to be Ly-$`\alpha `$ galaxies, and some presumed Ly-$`\alpha `$ galaxies may turn out to be $`[\mathrm{O}\mathrm{II}]\lambda \mathrm{\hspace{0.17em}3727}`$-emitting galaxies. However, any asymmetry in this mutual contamination has already been accounted for by assuming a conservative estimate of 50% contamination of the Ly-$`\alpha `$ candidates sample.

Surface brightness depends on the apparent size of an object at redshift $`z`$, which does not change much between $`z=3.5`$ and $`z=5.7`$, and on the luminosity distance. Since the model luminosity function is converted to $`z=5.7`$ taking the luminosity distance into consideration, most of the effect of the $`(1+z)^4`$ dimming towards higher redshift has already been accounted for.

Large-scale structure can influence galaxy counts even at high redshifts (Steidel et al. steidel98 (1998)). Since we derive the number counts of Ly-$`\alpha `$ galaxies at $`z5.7`$ by combining four CADIS fields (with a total area about 14 times the size of the field searched by Hu et al. hu98 (1998)), large scale structure should average out. Nevertheless, it should be kept in mind that results from small fields, like those searched by Hu et al. (hu98 (1998), hu99 (1999)), could be affected by large scale structure.

Two points should be noted about the bright end of the luminosity function at $`z3.5`$. First, the *observed* Ly-$`\alpha `$ line fluxes for the most luminous galaxies at $`z3.5`$ found by Hu et al. (1998) and Kudritzki et al. (2000) are on the order of $`10^{19}\mathrm{W}\mathrm{m}^2`$. Taking into account the higher luminosity distance at $`z=5.7`$ (compared to $`z=3.5`$) these galaxies would have apparent line fluxes of $`3\times 10^{20}\mathrm{W}\mathrm{m}^2`$ at $`z=5.7`$ (using a $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and H$`{}_{0}{}^{}=70`$km s-1 Mpc-1 cosmology); i.e., galaxies with intrinsically bright Ly-$`\alpha `$ emission, corresponding to observed line fluxes $`>3\times 10^{20}\mathrm{Wm}^2`$ at $`z5.7`$ $``$ as searched by CADIS $``$ have not yet been found in the small fields searched at $`z3.5`$. Thus, the shape of the bright end of the luminosity function relies mainly on our model. Second, the number counts at the bright end of the luminosity function at $`z=3.5`$, resulting from the brightest galaxy detected by Kudritzki et al. (2000, large circles), are higher than the results of Hu et al. (1998, small circles). Since large-scale structure can influence galaxy counts even at high redshifts (Steidel et al. 1998), the Kudritzki et al. number counts could be biased, because Kudritzki et al. searched only one field, whereas Hu et al. survey goes over two fields. Disregarding the brightest galaxy from Kudritzki et al. would result in a steeper luminosity function at the bright end at $`z=3.5`$ and thus lessen the discrepancy at $`z=5.7`$.

In summary, we conclude that the luminosity function of Ly-$`\alpha `$-bright galaxies declines between $`z=3.5`$ and $`z=5.7`$. This implies that the peak of the Ly-$`\alpha `$-bright phase $``$ i.e., the first formation of massive stars $``$ is reached after $`z=6`$. Although, in principle, both luminosity evolution or density evolution could account for this finding, we only can establish the decrease in density around $`3\times 10^{20}\mathrm{Wm}^2`$, as the bright end of the Ly-$`\alpha `$ luminosity function is still not determined well enough at $`z=3.5`$.

The onset of massive star formation should contribute a substantial fraction of the UV photons which reionize the universe. Moreover, both bright starbursts (SFR $`10\mathrm{M}_{}`$/yr) and the formation of supermassive black holes indicate a strong concentration of baryons within the centers of the most prominent density peaks. Therefore, our result seems to confirm additionaly the finding from the analysis of the HI absorption in the spectra of SDSS quasars at $`z>6`$ (Fan et al. fan02 (2002)), that the reionization of the universe did not happen long before $`z=6`$.


format: text

text:
# Recovering Isotropic Statistics in Turbulence Simulations: The Kolmogorov 4/5th-Law

## I Introduction

Both experimental and numerical studies of turbulence have attempted to observe the 1941 predictions of A.N. Kolmogorov K41b  for the statistics of isotropic, homogeneous, fully developed turbulence in the limit of infinite Reynolds number in an incompressible fluid. A main result of the theory is the so-called “4/5-law”

$`(\delta u_L(𝐫,𝐱))^3`$ $`=`$ $`{\displaystyle \frac{4}{5}}\epsilon r`$ (1)
$`\delta u_L(𝐫,𝐱)`$ $`=`$ $`[𝐮(𝐱+𝐫)𝐮(𝐱)]\widehat{𝐫}`$
$`\widehat{𝐫}`$ $`=`$ $`𝐫/r`$
where $``$ denotes ensemble averaging. The lefthand side of Eq. 1 is the well-known third-order longitudinal structure function. The length scale $`r`$ must lie in the inertial range $`\eta <<r<<L`$, sufficiently far from the large scales $`L`$ and the dissipation scales given by the Kolmogorov length $`\eta `$. The mean energy dissipation rate of the flow is given by $`\epsilon `$. The 4/5ths law is one of the few exact, non-trivial results known in the theory of statistical hydrodynamics. It may be reformulated in terms of other components of the structure function by using the incompressibility constraint MYII

$$\delta u_L(\delta u_T)^2=\frac{1}{6}\frac{}{r}(r(\delta u_L)^3)$$
(2)
where $`\delta u_T`$ is a velocity increment along a vector transverse to the separation vector $`𝐫`$. In Eq. 2 and henceforward, the vector argument $`𝐫`$ is implicit. This combined with the Eq.1 gives “4/15ths law” and the “4/3rds law”

$`\delta u_L(\delta u_T)^2={\displaystyle \frac{4}{15}}\epsilon r`$ (3)
$`\delta u_L|\delta 𝐮|^2={\displaystyle \frac{4}{3}}\epsilon r`$ (4)
where $`|\delta 𝐮|`$ denotes the total magnitude of the velocity difference across $`𝐫`$. We will refer to the three laws given by Eqs. (1), (3) and (4), and the related theory collectively as K$`41^3`$.

The K$`41^3`$ results have served as invaluable benchmarks for the empirical study of high-Reynolds number turbulence in both experiments and numerical simulations. Considered as exact results, they have allowed investigators to assess the degree to which homogeneous, isotropic, and high Reynolds number conditions have been attained. Furthermore, the derivation of the K$`41^3`$ relations requires a fundamental, unproven assumption, namely, that turbulent energy dissipation has a strictly positive limit as viscosity tends to zero. Hence, the validity of the K$`41^3`$ relations constitute an important test of this basic assumption. Experiments in high-Reynolds number turbulence performed over the past half-century do, by and large, support the linear scaling in $`r`$ of the third-order structure functions. The convergence to the 4/5ths coefficient is quite slow as Reynolds number increases for large scale anisotropic experiments SVBSCC96 ; BD01 , although there is empirical consensus that indeed this is asymptotically the correct coefficient. Recent numerical simulations GFN02  of isotropically forced Navier-Stokes also emphasize the slow approach to the 4/5ths law as the Reynolds numbers are pushed as high as computational power would allow. A key feature of both experimental and numerical endeavors, is the large volumes of data required – very long time-averages, extending over many integral length-scales or eddy-turnover times are needed to obtain adequate statistics and to observe the trend toward K$`41^3`$.

A modified version of the 4/5ths law, which does not assume isotropy of the flow, now exists. Nie and Tanveer NieTan99  proved that the 4/3rds and consequently, the 4/5ths laws can be recovered in homogeneous, but not necessarily isotropic flows,

$`(\delta u_L)^3`$ $`=`$ $`\underset{T\mathrm{}}{lim}{\displaystyle \frac{1}{T}}{\displaystyle _0^T}𝑑t{\displaystyle \frac{d\mathrm{\Omega }}{4\pi }\frac{d𝐱}{L^3}[\delta u_L(𝐫;𝐱,t)]^3}`$ (5)
$`=`$ $`{\displaystyle \frac{4}{5}}\epsilon r.`$
The angle integration $`d\mathrm{\Omega }`$ integrates in $`𝐫`$ over the sphere of radius $`r`$. For each point $`𝐱`$ the vector increment $`𝐫`$ is allowed to vary over all angles and the resulting longitudinal moments are integrated. The integration over $`𝐱`$ is over the entire flow domain. The integration over time $`t`$ extends over long-times, and the long time average is consistent with the ensemble average of the original K41 theory since ergodicity allows identification of ensemble-averages with time-averages Frisch95 . In Eq. 5, the integration over $`\mathrm{\Omega }`$ extracts the isotropic component of a generally anisotropic flow. This is fully consistent with recent experimental ADKLPS98 ; KLPS00  and numerical ABMP99 ; BLMT02  efforts to quantify anisotropic contributions by projecting the structure function onto a particular irreducible representation, labeled by $`j=0,1,\mathrm{}`$, of the SO(3) symmetry group. The angle-averaged Eq. 5 correspond to projecting onto the $`j=0`$ (isotropic) sector by integration over the sphere. The authors of NieTan99  do not perform numerically the average over the sphere. However, they do point out that the direction of $`𝐫`$ matters strongly. In their anisotropic DNS simulation at moderate Reynolds number, the result of taking $`𝐫`$ along a coordinate direction gave very poor agreement with the laws, whereas taking $`𝐫`$ along a body-diagonal as defined by the square-grid gave much better agreement.

A local version of the 4/3rds law was recently derived by J. Duchon and R. Robert DucRob99 . Subsequently, G.L. Eyink Eyink02  derived the corresponding version of the 4/5ths and 4/15ths laws. The statement is the following: Given $`any`$ local region $`B`$ of size $`R`$ of the flow, for $`r<<R`$, and in the limits $`\nu 0`$, next $`r0`$, and finally $`\delta 0`$,

$`(\delta u_L)^3_{(\mathrm{\Omega },B)}`$ $`=`$ $`\underset{\delta 0}{lim}{\displaystyle \frac{1}{\delta }}{\displaystyle _t^{t+\delta }}𝑑\tau {\displaystyle \frac{d\mathrm{\Omega }}{4\pi }_B\frac{d𝐱}{R^3}[\delta u_L(𝐫;𝐱,\tau )]^3}`$ (6)
$`=`$ $`{\displaystyle \frac{4}{5}}\epsilon _Br.`$
for almost every (Lebesgue) point $`t`$ in time, where $`\epsilon _B`$ is the instantaneous mean energy dissipation rate over the local region $`B`$. This version of the K$`41^3`$ results does not require isotropy or homogeneity of the flow. Long-time or ensemble averages are also not required as in the original K41 theory Frisch95 . The Duchon-Robert DucRob99  and Eyink Eyink02  versions of K$`41^3`$ are truly local in space and time.

We are motivated in the present work by the existence of isotropic statistics embedded in anisotropic data as suggested by previous work. It is clear that both experiments and simulations face intrinsic difficulties in achieving the high-Reynolds numbers and isotropic limit required by K41 theory. Both anisotropy and finite Reynolds numbers conspire to shorten the inertial range. Experiments have achieved Reynolds numbers several orders of magnitude higher than simulations. The indication is that at such high Reynolds numbers, the large scale anisotropies decay faster than the isotropic scales, allowing the latter to dominate at small scales KurSre00 . However, while the linear scaling in $`r`$ of the third-order structure function is fairly robust, the coefficient exhibits only a slow trend toward 4/5 as indicated by the numerical work of GFN02 . It is clear that for anisotropic forcing, some choices of directions for the vector increment $`𝐫`$ are more “isotropic” than others NieTan99 .

The concept of averaging over the sphere in order to extract the isotropic component of turbulence data has existed for some time ALP99 ; KSreview . Present high-Reynolds number experiments provide limited data – often only a few spatial points of data acquisition, with vector increment directions limited by the location of the probes and the implementation of a suitable space-time surrogation (Taylor’s hypothesis). Such configurations are not suitable to spherical averaging. Numerical data has, in principle, complete space-time information of the flow. However, the interpolation of square grid data over spherical shells has been deemed too expensive NieTan99 , or, when some such interpolation scheme is implemented, has not been used at sufficiently high Reynolds numbers as will allow for observations of the K41 type ABMP99 . The new angle-averaged and “local” laws of NieTan99 ; DucRob99 ; Eyink02  provide us with the theoretical impetus to investigate and extract the isotropic component of the flow in high-Reynolds number anisotropic turbulence. We use a novel means of taking the average over angles which avoids the expense and effort of interpolating the square-grid data over spherical shells.

In section II we discuss the numerical method and describe the stochastic and deterministic numerical forcing schemes used in the past and reimplemented by us. In section III we present an easily implemented scheme to average over a finite number of angles. Using the second and third-order isotropy relations, we demonstrate that this scheme is a good approximation to the true spherical average. We then present results for the angle-averaged third-order structure functions computed both from a snapshot of the flow at a single instant in time, and from time-averaging the data from several snapshots. A summary and concluding discussion of the results is given in section IV.

## II Numerical Simulations

The numerical simulation of the forced Navier-Stokes equation for an incompressible flow are given by

$`_t𝐮+\omega \times 𝐮+\varphi `$ $`=`$ $`\nu ^2𝐮+𝐟`$ (7)
$`𝐮`$ $`=`$ $`0`$ (8)
where the vorticity $`\omega =\times 𝐮`$ and $`\varphi `$ is determined so as to maintain $`𝐮=0`$. The domain is a periodic box of side $`L=2\pi `$ with $`N=512`$ grid points to a side. A standard Fourier pseudo-spectral method is used for the spatial discretization and the equations are integrated in time using a fourth-order Runge-Kutta scheme. Aliasing errors from the nonlinear term are effectively controlled by removing all coefficients with wave-number magnitude greater than $`k_{max}=\frac{\sqrt{2}}{3}N`$. The code is optimized for distributed memory parallel computers and uses MPI for inter-process communication. The runs were made using 256 processors of a Compaq ES45 cluster.

We make use of two different types of low wave number forcing. The first is modeled after the deterministic forcing schemes described in Kerr94 ; SVBSCC96 ; OvePop98 , where the energy in a few low wave numbers is relaxed back to a target spectrum. We refer to the output using this forcing as the deterministic dataset. The second forcing is the stochastic forcing used in GFN02 , where the Fourier coefficients of $`𝐟`$ are chosen randomly, and we refer to the data produces with this forcing as the stochastic dataset. Both forcings have advantages and disadvantages. The deterministic forcing equilibrates quickly and has less variance in time so that less data is needed for converged time averages. But there is an unavoidable anisotropy throughout the simulation if the forcing is restricted to the lowest wave numbers. The stochastic forcing has a larger variance in time so that data from more snapshots is needed to obtain converged time averages, but statistics from those snapshots are observed to be more isotropic. We perform both kinds of forcing in order to demonstrate the equivalence of the results when angle-averaging is applied to the data.

Parameters of interest for both simulations are given in Table 1. For the stochastic forcing, we have chosen parameters similar to those used for the $`512^3`$ simulations in GFN02 . The parameters for the deterministic case were chosen so that $`R_\lambda `$ would be similar in both cases.

### II.1 Deterministic forcing

We first define the energy in each spherical wave number $`k`$ in the usual way:

$$E(k)=\underset{k.5|𝐤|<k+.5}{}\frac{1}{2}|\stackrel{~}{𝐮}_𝐤|^2$$
(9)
where $`\stackrel{~}{𝐮}_𝐤`$ is the $`k`$th Fourier coefficient of $`𝐮`$. We then choose a target spectrum function given by $`F(k)`$, which we set to $`F(1)=F(2)=.5`$ and $`F(k)=0`$ for $`k>2`$. We generate a velocity field $`\stackrel{~}{𝐯}_𝐤`$ with energy $`F(k)`$ by setting

$$\stackrel{~}{𝐯}_𝐤=\sqrt{\frac{F(k)}{E(k)}}\stackrel{~}{𝐮}_𝐤$$

The Fourier coefficients of the forcing function are then given by

$$\stackrel{~}{𝐟}_𝐤=\{\begin{array}{cc}\tau \left(\stackrel{~}{𝐯}_𝐤\stackrel{~}{𝐮}_𝐤\right),\hfill & F(k)>E(k)\hfill \\ 0,\hfill & F(k)E(k)\hfill \end{array}$$

The relaxation parameter is chosen as $`\tau ^1=2|u|`$, a simplified version of the formula given in OvePop98 .

This forcing simply relaxes the amplitudes of the Fourier coefficients in the first two wave numbers so that the energy matches the target spectrum $`F(k)`$ in those wave numbers. It has no effect on the phase of the coefficients. The phases are observed to change very slowly, giving rise to persistent anisotropy in the large scales.

### II.2 Stochastic forcing

Our second type of forcing is modeled after the stochastic scheme of GFN02  in which the wave-numbers $`|𝐤|2.5`$ are stochastically forced. This ensures that the phase of each forced mode changes sufficiently rapidly so that the large scales will be statistically isotropic. At the beginning of each time-step, we choose a divergence-free forcing function,

$$𝐟=\times ^1𝐠$$
(10)
where the Fourier coefficients of $`𝐠`$, denoted by $`\stackrel{~}{𝐠}_𝐤`$, are chosen randomly with uniformly distributed phase and Gaussian distributed amplitude. The variance in the Gaussian distribution is chosen so that

$$\underset{.5|𝐤|<1.5}{}|\stackrel{~}{𝐠}_𝐤|^2=\underset{1.5|𝐤|<2.5}{}|\stackrel{~}{𝐠}_𝐤|^2=18$$
(11)
## III Angle-Averaging Technique

We would like to extract the isotropic component of a flow by a suitable average of the two-point structure functions over the solid angle $`\mathrm{\Omega }`$ as defined by Eqs. (56). We approximate the spherically averaged third-order longitudinal structure function by the following average over a finite number $`N_d`$ of directions:

$$(\delta u_L(r))^3=\frac{1}{N_d}\frac{1}{N^3}\underset{j=1}{\overset{N_d}{}}\underset{i=1}{\overset{N^3}{}}w_j\left[\delta u_L(𝐫_j;𝐱_i)\right]^3.$$
(12)
where $`𝐱_i`$ denotes grid-points, $`𝐫_j`$ denotes the increment vector in the $`j`$th direction, $`r=|𝐫_j|`$ is fixed, and the $`w_j`$ are quadrature weights. Here we are using the longitudinal structure function as an example. The procedure applies equally well to any two-point structure function.

The simulation is computed on a fixed uniform rectangular mesh. Thus we are faced with the difficulty of evaluating $`𝐮`$ at points $`(𝐱_i+𝐫_j)`$, most of which will not be grid points. The most straightforward approach would be to perform 3-dimensional interpolations at each of the points $`(𝐱_i+𝐫_j)`$. This requires $`N^3N_d`$ 3-D interpolations of the velocity-vector field for each separation distance $`r`$, which is prohibitively expensive NieTan99 .

We have developed a less expensive technique for angle-averaging which does not require any 3-D interpolations. We first choose the vectors $`𝐫_j`$ from among those natural to a square computational grid. We restrict ourselves to the set of all unique directions which can be expressed with integer components with length less than or equal to $`\sqrt{11}`$. Let $`j=1\mathrm{}N_d`$ be the index for this set. Each $`𝐫_j`$ is the minimum grid-point separation distance in the $`j`$th direction. This set is generated by the vectors (1,0,0), (1,1,0), (1,1,1), (2,1,0), (2,1,1), (2,2,1), (3,1,0) and (3,1,1) by taking all index and sign permutations of the three coordinates, and removing any vector which is a positive or negative multiple of any other vector in the set. This procedure generates a total of $`N_d=73`$ unique directions. The unit vectors associated with each direction are plotted as points on the sphere in Fig. 1. One can see that these points are well distributed over the sphere. Both the unit vectors $`\widehat{𝐫}_j`$ and $`\widehat{𝐫}_j`$ are plotted, but below we do not consider the $`\widehat{𝐫}_j`$ directions since they give the same contribution as $`\widehat{𝐫}_j`$ when averaged over the periodic computational domain.

For each of the $`N_d`$ directions, we form a set of $`\mathrm{}=1\mathrm{}N_r`$ separation vectors, $`𝐱_i+\mathrm{}𝐫_j`$. Since $`𝐫_j`$ is the minimum separation distance of grid-points in the $`j`$th direction and $`\mathrm{}`$ is an integer, all the $`𝐱_i+\mathrm{}𝐫_j`$ lie on our computational grid. This is illustrated (in two-dimensions) for four directions in Fig. 2, where the black dots represent the points $`𝐱_i+\mathrm{}𝐫_j`$ and $`𝐱_i`$ is shown at the origin. We can now efficiently compute structure functions in $`N_d`$ different directions, at $`N_r`$ separation distances for each direction, without any 3-d interpolations:

$$(\delta u_L(\mathrm{}𝐫_j))^3=\frac{1}{N^3}\underset{i=1}{\overset{N^3}{}}\left[\delta u_L(\mathrm{}𝐫_j;𝐱_i)\right]^3.$$
(13)
For each direction, we get a 1-dimensional curve as a function of $`\mathrm{}𝐫_j`$, as shown in Fig. 3.

In the figure, points represent structure function values at the separation distances $`\mathrm{}|𝐫_j|`$, and each line is a cubic-spline fit to the data at $`\mathrm{}=1\mathrm{}N_r`$ along each of the $`N_d`$ directions. One can see that only a few directions are computed at each of the separation distances, so we cannot directly take an angle-average from this data. But one can also see that the curves are quite smooth and the cubic-spline is an excellent interpolant. Thus we use cubic-spline interpolation to calculate the structure function in each of the $`N_d`$ directions, at separation vector $`r\widehat{𝐫}_j`$ of any desired length $`r`$.

Once the data for each direction has been interpolated to a common separation distance $`r`$, we can approximate the angle-average at $`r`$ by quadrature over the $`N_d`$ directions:

$$(\delta u_L(r))^3=\frac{1}{N_d}\underset{j=1}{\overset{N_d}{}}w_j(\delta u_L(r\widehat{𝐫}_j))^3$$
(14)
In order to determine the quadrature weights $`w_j`$, we use the software package Stripack STRIPACK  to compute the Voronoi tiling generated by the points $`\widehat{𝐫}_j`$ on the unit sphere centered at $`𝐱`$. The weight $`w_j`$ is the solid angle subtended by the Voronoi cell containing the point $`\widehat{𝐫}_j`$.

The angle-averaging procedure described can be implemented efficiently on parallel computers, requiring only the same type of parallel data transpose operator already used by a parallel pseudo-spectral code. The total cost of this angle-averaging procedure for one snapshot (73 directions and 100 different separation distances) is about the same as 150 timesteps of the Navier-Stokes code. Thus for a single eddy turnover time, where thousands of timesteps are required, the angle-averaging statistics can be computed during the computation with minimal impact on the total CPU time requirement.

### III.1 Extracting the isotropic component

We first present results demonstrating how well the angle-average procedure performs at extracting the isotropic component from our DNS data. We again follow GFN02  and examine the relations between the second and third order velocity structure functions:

$`(\delta u_T)^2=(1+{\displaystyle \frac{r}{2}}{\displaystyle \frac{d}{dr}})(\delta u_L)^2`$ (15)
$`\delta u_L(\delta u_T)^2={\displaystyle \frac{1}{6}}{\displaystyle \frac{d}{dr}}r(\delta u_L)^3`$ (16)
These equations require only isotropy and incompressibility. Thus in DNS data, where incompressibility is obtained to numerical round off error, deviations in the above relations are a measure of the anisotropy in the data. In GFN02 , the left and right sides of these equations are plotted after averaging in time. Excellent agreement is obtained in the inertial range, with some departure at larger scales.

In Fig. 4, we show the second order isotropy relation for our stochastic dataset, and in Fig. 5 we show the third order relation for the deterministic dataset. This data is computed by angle-averaging over a single snapshot of the flow. The agreement is excellent, both in the inertial range and at the largest scales. For comparison, the figures also show the same relations from the same snapshot but using only a single coordinate direction instead of angle-averaging. In that case, there are significant differences for scales well into the inertial range. Thus the angle-averaging technique appears to be extremely effective in extracting the isotropic component of anisotropic data, even at large scales where anisotropy remains after time averaging over many snapshots. Similar results were obtained for the second order isotropy relation from the deterministic dataset and for the third order isotropy relation from the stochastic dataset.

### III.2 Angle-averaging a single snapshot

We now present results from using angle-averaging to compute the third-order longitudinal structure function in the 4/5ths law. Figures 6 and 7 show the result of the angle-averaging procedure described above for single snapshots of the stochastic and deterministic datasets respectively. The snapshots are taken after the flow has had time to equilibrate. The value of the mean energy dissipation rate $`\epsilon `$ was calculated from the snapshot. This is to be contrasted with previous works in which $`\epsilon `$ is a long-time or ensemble average. We have therefore computed a version of the 4/5ths relation which is local in time. The dots represent the data from all 73 directions at all values of $`r`$ that were computed. The final weighted angle-average of Eq. (14) is given by the thick curves in both Figs. 6 and 7. One can see that the results from different directions are quite different, while the angle-averaged results are quite reasonable and similar to each other as well as similar to the results obtained from long time averaging of the coordinate directions presented in GFN02  and shown for our data in section III.4. Thus we conclude that angle-averaging the data from a single snapshot yields a very reasonable result. Similar results (not plotted) are obtained for the 4/3’rds and 4/15’ths laws.

### III.3 Temporal variance

To illustrate the variance in time of the third-order longitudinal structure function, with and without angle-averaging, we plot the peak value as a function of time for each dataset in Figures 8 and 9. The solid line is the angle-averaged value, and the dashed line is the value from a single coordinate direction. The angle-averaged value has a significantly reduced variance as compared to the single direction value, but one can see that there still is some variance from snapshot to snapshot. Thus in order to obtained fully converged statistics, some additional averaging is needed. In the next section we present results combining angle and time averaging.

Based on the local version of the K$`41^3`$ laws proved in DucRob99 ; Eyink02 , we expect that increasing the spatial resolution would allow us to obtain converged statistics from a single snapshot when used with angle-averaging. However, we could not expect such convergence without angle-averaging. This is because even in an isotropic flow, individual snapshots are not necessary isotropic - only the ensemble of all snapshots is guaranteed to be isotropic.

We note that the stochastic dataset (Fig. 9) shows a larger variation from snapshot to snapshot when compared to the deterministic dataset (Fig. 8). This is true for both the angle-averaged and single direction quantities shown in the figures, suggesting that the stochastic dataset produces data with a slightly larger variance in time, as expected.

### III.4 Time averaged results

We now look at the 4/5ths law using both angle-averaging (which extracts the isotropic component of the statistics) and time-averaging (to remove the variance observed from snapshot to snapshot). The time-average is taken from 60 snapshots taken over 6 eddy turnover times. The results are shown in Fig. 10. The two datasets produce nearly identical results at all scales, even though the large scale forcing is quite different. The peak value of the stochastic and deterministic datasets are .755 and .752, respectively.

Thus we conclude that flows with similar geometry and Reynolds number have the same underlying isotropic component at all scales, at least up to third order statistics.

## IV Conclusions

We have proposed a new, computationally efficient and easily implemented means of extracting isotropic statistics from an arbitrarily forced flow. As a first test of the method, we averaged the third-order structure functions over sufficiently many angles and discovered that the K$`41^3`$ relations are obtained, with tolerable variance, from a $`single`$ snapshot of homogeneous flow with either stochastic or deterministic forcing. This is a stronger result than was predicted by the original Kolmogorov ensemble approach or even the Nie-Tanveer version of NieTan99 . It appears that the results are, in fact, approaching the $`local`$ versions of K$`41^3`$ proposed in DucRob99 ; Eyink02 .

Using our procedure to extract the isotropic component, we are able to separate the effect of anisotropy from the effect of finite-Reynolds number on the statistics of the flow. This is an important point to make in the debate on how the two effects contaminate the inertial range. Once the anisotropy is eliminated, a more fruitful study of finite-Reynolds number effects can be made. It is clear from Fig. 10 that the Reynolds numbers are still not sufficient to give the wide inertial ranges that have been seen in high-Reynolds number experiments. However, it is also clear that angle-averaging has given a significant improvement in the results. With angle-averaging, less data is needed to obtain converged statistics, and deterministic forcings can be used without regard to the increased anisotropy they introduce.

The procedure we have described above can be used to investigate the isotropic component of higher-order structure functions or any other statistic as well. For example, the angle-averaged $`n`$th-order longitudinal structure functions may be measured in this way in order to determine scaling exponents which are truly independent of anisotropy. This method may also be used to isolate the anisotropic contributions themselves, as has been done in ABMP99 ; BLMT02 , by subtracting from the full structure function its angle-averaged value. Individual moments in a spherical harmonics expansion of structure functions can be computed by introducing the basis function of interest to the integrand in equation 14. In this way, the dominant scaling in anisotropic sectors can be determined, which is important to determine the rate of return to isotropy at small-scales. We plan to investigate such questions in future work.

###### Acknowledgements.

We thank Toshiyuki Gotoh for his assistance with the stochastic forcing procedure and for fruitful discussions.