text:
# Self-Consistency and Calibration of Cluster Number Count Surveys for Dark Energy

## I Introduction

The number density of massive galaxy clusters is exponentially sensitive to the amplitude of initially Gaussian density fluctuations and has long been recognized as a sensitive cosmological probe Cluster . In particular, the evolution of the number counts above a given mass threshold can determine the properties of the dark energy that accelerates the expansion WanSte98 ; Dieetal01 ; HaiMohHol01 .

This potential can only be realized if, in addition to the cluster redshifts, the cluster masses themselves are known, at least statistically. Unfortunately the total mass is not a direct observable and must be inferred through scaling relations with, e.g. the Sunyaev-Zel’dovich (SZ) flux decrement, the $`X`$-ray flux or temperature, weak lensing shear, or optical velocity dispersion. Indeed the normalization of the mass-temperature relation is currently the leading source of ambiguity in interpreting the local cluster abundance (e.g. PieScoWhi01 ; Sel02 ). It compromises dark energy constraints if the cluster mass selection cannot be defined to a few percent in mass HaiMohHol01 .

The mass-observable relation can potentially be calibrated within a survey itself if its effect on the number counts is not degenerate with the cosmology. For example, with a single cut on the temperature, the survey itself can calibrate the mass-temperature normalization to better than the tens of percent that span the determinations in the current literature, if it does not evolve LevSchWhi02 . However this method fails if the relationship has uncertain evolution that mimics the cosmology MajMoh02 .

Because the number density of clusters as a function of mass has a fixed functional form given by cosmological simulations (e.g. Jenetal01 ; Reeetal03 ), cluster number counts as a function of the observable in principle have the ability to self-calibrate even an evolving mass-observable relation. Here we study the potential for internal calibration and consistency checks in idealized future cluster number count surveys.

## II Statistical Forecasts

The cosmological utility of cluster counts stems from the simulation-based prediction of their comoving differential number density as a function of mass Jenetal01 ,

$$\frac{d\overline{n}}{d\mathrm{ln}M}=0.3\frac{\rho _m}{M}\frac{d\mathrm{ln}\sigma ^1}{d\mathrm{ln}M}\mathrm{exp}[|\mathrm{ln}\sigma ^1+0.64|^{3.82}].$$
(1)
Here $`\sigma ^2(M,z)`$ is the variance in the linear density field smoothed with a top hat that encloses $`M`$ at the mean matter density today $`\rho _m`$. The cosmological sensitivity comes from these quantities and the comoving volume element in a redshift slice and solid angle.

To assess the sensitivity of counts to various cosmological and mass-observable scaling parameters $`p_\alpha `$, we employ the Fisher matrix technique HolHaiMoh01 ; HuKra02

$$F_{\alpha \beta }=\underset{ij}{}\frac{\overline{n}_i}{p_\alpha }(𝐂^1)_{ij}\frac{\overline{n}_j}{p_\beta },$$
(2)
where the covariance matrix is given by

$$C_{ij}=\left(n_in_j\overline{n}_i\overline{n}_j\right)+\delta _{ij}\overline{n}_i/V_i.$$
(3)
The first term represents sample covariance in the volume $`V_i`$ from large-scale structure, calculated as described in HuKra02 , and the second term shot variance. Here the number density in a bin $`n_i`$ is defined by the redshift interval $`\mathrm{\Delta }z`$ around $`z_i`$ and a selection based on some observable quantity $`f`$, such as the SZ flux decrement. Note that sampling errors for the selections in the same redshift bin completely covary.

The Fisher matrix is a local approximation to the covariance matrix of the parameters

$$𝐂_{\mathrm{tot}}=(𝐅+𝐂_{\mathrm{prior}}^1)^1,$$
(4)
where $`𝐂_{\mathrm{prior}}`$ is the covariance matrix from prior information.

The Fisher matrix is evaluated around a fiducial choice of parameters. For the cosmology, we take: the dark energy density $`\mathrm{\Omega }_{\mathrm{DE}}=0.65`$, equation of state $`w=1`$, physical matter density $`\mathrm{\Omega }_mh^2=0.148`$, physical baryon density $`\mathrm{\Omega }_bh^2=0.02`$, tilt $`n_s=1`$, and the initial normalization of the curvature fluctuations $`\delta _\zeta =4.79\times 10^5`$ at $`k=0.01`$ Mpc-1 Hu01c  (more conventionally: $`\delta _H=4.42\times 10^5`$, $`\sigma _8=0.92`$, or $`M_{}=1.2\times 10^{13}h^1M_{}`$). We will however work in a future context where the cosmic microwave background (CMB) has constrained the high redshift universe parameters with 1$`\sigma `$ errors of $`\sigma (\mathrm{ln}\mathrm{\Omega }_mh^2)=\sigma (\mathrm{ln}\mathrm{\Omega }_bh^2)=\sigma (n_s)=\sigma (\mathrm{ln}\delta _\zeta )=0.01`$ (e.g. Hu01c ) in a flat universe, leaving the dark energy parameters to be determined by the survey.

For the mass-observable relation, we take the form

$$\frac{M_f}{M_0}=e^{A(z)}\left(\frac{f}{f_0}\right)^{p(z)},$$
(5)
where $`M_0`$ and $`f_0`$ are dimensional constants and the normalization parameter $`A(z)`$ is dimensionless. We take the scaling parameters as the amplitudes of a piecewise constant form $`A(z_i)`$ and $`p(z_i)`$. Specifically, given a variation $`\delta A(z_i)`$ and $`\delta \mathrm{ln}p(z_i)`$ from the fiducial model, the true mass is related to the apparent mass as

$$\mathrm{ln}M=\mathrm{ln}M_f+\delta A(z_i)+\delta \mathrm{ln}p(z_i)\mathrm{ln}\left(\frac{M_f}{M_p}\right),$$
(6)
where $`M_p=M_0e^A`$ is the pivot mass scale that defines a normalization point for possible external information.

This generalizes previous treatments which have assumed constant normalization LevSchWhi02  and power law evolution MajMoh02 . In the limit of small redshift bins, these simpler cases can be recovered from our more general treatment by noting that under a re-parameterization of the space to the set $`\pi _\mu (p_\alpha )`$, possibly of lower dimension, the covariance matrix transforms as

$$C_{\mu \nu }=\underset{\alpha \beta }{}\frac{\pi _\mu }{p_\alpha }C_{\alpha \beta }\frac{\pi _\nu }{p_\beta }.$$
(7)
The Fisher matrix definition itself can be so viewed.

In the usual approach (e.g. HaiMohHol01 ), all of the clusters in a given redshift bin above a given observable threshold are simply binned together. Clearly, allowing for arbitrary variations $`\delta A(z_i)`$ no cosmological information can be extracted. However the data at a given redshift may be binned into several apparent mass bins based on the observable $`f`$. The additional information supplied by multiple bins allows for a breaking of the degeneracy.

Likewise, even if the normalization is fixed by detailed followup by $`X`$-ray or weak lensing measurements at some mass $`M_p`$ (e.g. MajMoh02 ; HutWhi02 ), an uncertain scaling index $`\delta p(z_i)`$ would again destroy the cosmological information with a single bin if $`M_pM_{f,\mathrm{min}}`$. Multiple bins again in principle allow the cosmology and the mass-observable relations to be jointly determined.

This treatment ignores intrinsic scatter in the mass-observable relation as well as measurement error in the observable. More realistically the selection function is not sharp in mass but assuming its functional form is well characterized this does not substantially compromise the cosmological information Holetal00 ; LevSchWhi02 .

For definiteness, let us take a fiducial cluster survey with specifications similar to the planned South Pole Telescope (SPT) Survey for clusters with the SZ effect: an area of 4000 deg2 survey and a sensitivity corresponding to a constant $`\mathrm{log}_{10}M_{f,\mathrm{min}}/h^1M_{}=14.2`$ (e.g. MajMoh02 ). We divide the number counts into bins of redshift $`\mathrm{\Delta }z=0.1`$ out to $`z_{\mathrm{max}}=3`$ and $`\mathrm{\Delta }\mathrm{log}_{10}M_f=0.43\mathrm{\Delta }\mathrm{ln}M_f=0.2`$. This crude binning is sufficient to retain the cosmological information. For illustration purposes, we take a pivot mass scale of $`\mathrm{log}_{10}M_p/h^1M_{}=14.7`$ to reflect potential mass measurement followup on the high mass end.

## III Consistency and Calibration

Let us first consider that simulations or mass followup with different observables have placed strong priors on the mass-observable normalization and scaling index $`\sigma (\delta A(z_i))=\sigma (\delta \mathrm{ln}p(z_i))=0`$. In this case, the division into mass bins yields nearly independent measurements of the dark energy for a consistency check.

In Fig. 1, we illustrate the consistency check in the $`(\mathrm{\Omega }_{\mathrm{DE}},w)`$ plane for individual mass bins of $`\mathrm{\Delta }\mathrm{log}_{10}M_f=0.2`$ from $`\mathrm{log}(M_{f,\mathrm{min}}/h^1M_{})=14.2`$. Although the strongest constraints do come from the lowest mass bin, as expected from the increasing rarity of massive clusters at high redshift, the first 6 bins show significant constraints on the dark energy. This consistency test demonstrates that there is sufficient information in the mass binning to attempt some degree of self-calibration of the mass-observable relations.

Let us next keep the scaling index fixed $`p(z_i)`$ but allow the normalization $`A(z_i)`$ to vary (see Fig. 2). Solid lines represent constraints with the 6 mass bins and dashed lines those with a single mass bin or mass threshold. We first allow only a redshift-independent normalization factor $`\delta A(z_i)=\delta A_0`$. As shown in LevSchWhi02 , even a single mass bin is sufficient to calibrate the relation to good enough accuracy for dark energy constraints (here $`\sigma (\delta A)=0.05`$). We next employ a power law evolution in the normalization $`\delta A(z_i)=\delta A_0+n_a\mathrm{ln}(1+z_i)`$. Here errors on the dark energy parameters degrade substantially with a single mass bin MajMoh02 , but much of the loss is recovered from the multiple binning. Finally, we take a fully arbitrary evolution in the redshift bins, i.e. no constraint on the form of $`\delta A(z_i)`$. As expected, there is no constraint on the dark energy with a single mass bin, but with multiple bins interesting constraints on the dark energy can still be extracted. The key here is that the dark energy evolution is assumed to be smooth and parameterized by two numbers and so cannot compensate an arbitrary variation in the normalization evolution.

Figure 3 (top) shows the errors on $`\delta A(z_i)`$ with roughly 10% constraints on the normalization in each redshift band out to $`z2`$. Note that the covariance between the bands is nearly complete so that the overall normalization is still only known to $`5\%`$ but that the evolution is highly constrained. Conversely, external information on the mass-observable normalization would have to be substantially better than $`10\%`$ to improve on the internal results. In the context of our fiducial cosmology, the critical prior assumption is that the initial amplitude of the fluctuations is fixed by the CMB. Otherwise what is constrained by the multiple mass bins is a degenerate combination of the initial amplitude and the mass-observable normalization.

We have so far assumed that the scaling index $`p(z)`$ is known a priori but of course that too is part of the uncertain mass-observable relation. Allowing for arbitrary variations $`\delta A(z_i)`$ and $`\delta \mathrm{ln}p(z_i)`$ destroys most of the information on the dark energy. However in the less drastic cases of adding a constant $`\delta \mathrm{ln}p`$ to a constant $`\delta A`$ for a two parameter model, the additional degradation $`\sigma (w)`$ is a negligible $`3\%`$ and of the addition of power law variation in $`p(z)`$, $`\delta \mathrm{ln}p=\delta \mathrm{ln}p_0+n_p\mathrm{ln}(1+z_i)`$ to power law variation in $`\delta A`$ for a 4 parameter model, a factor of $`1.76`$ for a total of $`\sigma (w)=0.12`$.

Another interesting case to consider is if some independent mass calibration, say from weak lensing or $`X`$-ray temperature followup, normalizes the mass-observable relation on the high mass end. Extrapolation down to the survey mass limit can be dangerous due to the uncertain physics of low mass clusters. For illustrative purposes, let us take this constraint as $`\sigma (A(z_i))=0.1`$ at the pivot mass of $`\mathrm{log}_{10}M_p/h^1M_{}=14.7`$. Despite the constraint on the normalization, if we allow for arbitrary variations in the scaling index and retain only a threshold at $`>14.2`$, no constraint on the dark energy can be extracted. With multiple bins, the errors on $`w`$ are only degraded by a factor of 1.16 to $`\sigma (w)=0.09`$ for a completely arbitrary $`\delta \mathrm{ln}p(z_i)`$. Furthermore the errors on $`\delta \mathrm{ln}p`$ shown in Fig. 3 (bottom) are at the several percent level and would be of interest in studying cluster physics.

## IV Discussion

As in the case of classical cosmological tests for the dark energy involving standardized candles and rulers, e.g. supernovae and the peaks in the CMB power spectrum, cluster number count tests require a standardized mass based on observable quantities. We have demonstrated that consistency with the well-determined shape of the mass function from cosmological simulations can in principle be used to calibrate the survey internally. Two examples are a normalization that has an arbitrary evolution in redshift and a scaling from the more easily calibrated high mass end that has an arbitrary evolution in the index.

Our study involves a number of idealizations that merit future study and so internal calibration is best viewed as a useful check on cross-calibration studies with detailed multi-wavelength followup (e.g. MajMoh02 ; VerHaiSpe02 ) and cosmological simulations. The crucial assumptions are that the survey be sufficiently deep to explore a substantial dynamic range in the mass function, the underlying mass function is known given a cosmology, the mass-observable selection function is sharp compared with the binning, the mean mass-observable relation is power law in form, the high redshift normalization of the fluctuations is fixed by future CMB data, and the dark energy equation of state is constant.

In the fiducial model, a depth comparable to that planned for the SPT SZ survey is sufficient to make self-calibration useful but note that the necessary mass limit scales roughly with the non-linear mass scale $`M_{}`$. That depends on the currently uncertain normalization and so if the normalization is lower than our fiducial choice it may be beneficial to sacrifice survey width for depth.

While the mass function form is based on simulations and currently its scaling with cosmology is only known to the $`10\%`$ level Jenetal01 , determination requires only the well-understood gravitational physics of the dark matter and is far more secure than simulation-based mass-observable relations. In any case, this knowledge is a prerequisite for any cluster number count study of cosmology.

We have assumed a deterministic mass-observable relation that makes the selection function sharp in apparent mass. Gaussian scatter in this relation does not seriously compromise the cosmological information Holetal00 ; LevSchWhi02  but long uncharacterized tails to low mass would, given the steep mass function. Likewise the cluster physics of preheating and cooling can enter into the low mass end to make the mass-observable relations deviate from a power law. Here the binning serves as a useful consistency test for the implicitly assumed cluster physics.

The required prior information on the initial normalization is within reach of the upcoming CMB satellite missions if the extent of reionization can be determined (e.g. Hu01c ; Kapetal02 ). If not, binning still serves to constrain the evolution of the mass-observable relation but leaves a degeneracy between the overall normalization of the relations and the fluctuations which can be fixed, e.g. by CMB lensing, Ly$`\alpha `$ forest clustering or even local cluster abundance studies. Conversely, external calibration of the mass-observable relation at the percent level can determine the fluctuation normalization at a comparable level. Binning can then provide tests of the dark energy model assumptions WelBatKne01  or other cosmological priors. It is especially valuable if redshifts are available only locally (e.g. $`z_{\mathrm{max}}0.6`$ from current optical surveys).

In summary, utilizing the extra information in counting clusters as a function of both the redshift and their observable properties allows for joint solutions to the mass-observable relations and cosmology, if they are both simple. If not, it offers valuable internal consistency checks against overly simplistic assumptions.

Acknowledgments: I thank Z. Haiman, A.V. Kravtsov, J.J. Mohr, C. Pryke and the Chicago Thunch group for useful conversations and J. Valdes for \*nix expertise in file recovery. WH is supported by NASA NAG5-10840 and the DOE OJI program.


format: text

text:
# Vortex crystallisation in classical field theory

## Abstract

We show that the formation of a vortex lattice in a weakly interacting Bose condensed gas can be modeled with the nonlinear Schrödinger equation for both $`T=0`$ and finite temperatures without the need for an explicit damping term. Applying a weak rotating anisotropic harmonic potential we find numerically that the turbulent dynamics of the field produces an effective dissipation of the vortex motion and leads to the formation of a lattice. For $`T=0`$ this turbulent dynamics is triggered by an already known rotational dynamic instability of the condensate. For finite temperatures, noise is present at the start of the simulation and allows the formation of a vortex lattice at a lower rotation frequency, the Landau frequency. These two regimes have different vortex dynamics. We show that the multimode interpretation of the classical field is essential.

Several groups have recently observed the formation of a vortex lattice in weakly interacting Bose gases ENS ; MIT ; JILA ; Foot . Theoretically there have been attempts to understand the formation process Feder ; Ueda ; Gardiner ; Lundh  with simulations of the Gross-Pitaevskii equation for the condensate wavefunction. These papers either do not consider the effect of the noncondensed modes Feder ; Lundh  or they approximate it with an added damping term. All of them stress the need for explicitly including this term since the vortices have to lose energy to be able to relax to a lattice configuration.

Here we consider this problem within the framework of the classical theory of a complex field whose exact equation of motion is the standard nonlinear Schrödinger equation (NLSE). Within this framework we can give the true numerical answer which can be used as a testing ground for approximations such as the single mode one of Ueda ; Gardiner  to study the formation of the vortex lattice. If the dynamics of the field were ergodic we would expect that, in a fast enough rotating external potential, the vortex lattice would be formed without the need for any dissipative term. The issue of relaxation to thermal equilibrium for purely Hamiltonian dynamics is analogous to the Fermi-Ulam-Pasta problem Fermi  and has been studied also in the context of Bose gases Svistunov ; Sachdev ; Burnett ; Castin . Here we study the formation of the lattice in three dimensions from an initially nonrotating Bose condensed gas both at $`T=0`$ and at finite temperatures using the NLSE. We find that, contrary to the prevalent view, the lattice is formed in both cases without any need for a damping term which suggests that thermalisation takes place in our system.

We start our simulations with the nonrotating classical field in thermal equilibrium. For $`T=0`$ the initial state is simply a pure condensate, that is, with a field proportional to the condensate wavefunction $`\varphi `$ given by the Gross-Pitaevskii equation in the absence of rotation, $`\psi =\sqrt{N_0}\varphi `$ where $`N_0`$ is the condensate atom number. For finite temperatures we sample the initial thermal equilibrium in the Bogoliubov approximation at a given temperature $`T`$ for a fixed number $`N_0`$ of condensate particles. In this approximation the classical field is given by $`\psi (𝐫,0)=\sqrt{N_0}\varphi (𝐫)+\psi _{}(𝐫)`$. The random field $`\psi _{}(𝐫)`$ orthogonal to $`\varphi `$ Cartago  representing the thermal noise is given by

$$\psi _{}(𝐫)=\underset{n}{}b_nu_n(𝐫)+b_n^{}v_n^{}(𝐫)$$
(1)
where the $`u_n`$ and $`v_n`$ are the Bogoliubov mode functions associated with $`\varphi `$ and the $`b_n`$ are independent random c-numbers taken from a Gaussian distribution that obeys the classical equipartition formula, $`b_n^{}b_n=k_BT/ϵ_n`$, $`ϵ_n`$ being the Bogoliubov energy of mode $`n`$. In practice, to sample this distribution we use the Brownian motion method described in Cartago . In our work, the field $`\psi `$ is to be interpreted not as the condensate wavefunction but as the whole matter field (unlike Feder ; Gardiner ; Ueda ; Lundh ). We present here results from single realisations of the field $`\psi `$ which experimentally correspond to single runs. We have checked that different realisations lead to similar results.

In our simulations we consider a Bose condensed gas initially trapped in a cigar-shaped harmonic potential with oscillation frequencies whose ratio is 1:1:0.25, with $`10^5`$ atoms of mass $`m`$ and a coupling constant $`g=0.0343`$ in units of $`\mathrm{}\omega a_0^3`$ where $`\omega `$ is the radial frequency and $`a_0=\sqrt{\mathrm{}/m\omega }`$ is the oscillator length. The corresponding chemical potential is $`\mu =8\mathrm{}\omega `$. We start each simulation with the gas in thermal equilibrium. We abruptly turn on the trap anisotropy which leads to a change in the radial frequencies: $`\omega _{x,y}^2=\omega ^2(1ϵ)`$ where $`ϵ=0.025`$. Then the rotation frequency $`\mathrm{\Omega }(t)`$ of this anisotropy is slowly increased from zero to a final value $`\mathrm{\Omega }_f`$ over $`500\omega ^1`$, to follow Procedure I in Subhasis . After that we let the gas evolve in the presence of the rotating anisotropy until, for most of our simulations, the angular momentum of the gas reached a steady state.

The calculation is performed in the rotating frame so that the NLSE takes the form

$$i\mathrm{}_t\psi =\left[\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U(𝐫)+g|\psi |^2\mathrm{\Omega }(t)L_z\right]\psi $$
(2)
where $`L_z`$ is the angular momentum operator along $`z`$ and $`U`$ is the anisotropic harmonic potential. The field $`\psi `$ is subject to periodic boundary conditions in the rotating frame. Our gridsize is $`32\times 32\times 128`$ corresponding to an energy cutoff of $`32\mathrm{}\omega `$ per spatial direction, although we have also run simulations on a $`64\times 64\times 256`$ grid (see below). We have checked that it is necessary to have time independent boundary conditions in the rotating frame: periodic boundary conditions in the lab frame arrest the rotation of the non-condensed gas. We have also checked that a pure condensate cannot be set into rotation by the effect of the boundary conditions only, since the condensate density is extremely weak at the borders of the grid for our choice of gridsize. The harmonic trap anisotropy is then a crucial element for the formation of the lattice at $`T=0`$ by triggering the dynamic instability of Subhasis .

Zero initial temperature: This set of simulations can be divided into two groups: those for which the final rotation frequency is $`\mathrm{\Omega }_f/\omega 0.7`$ and those with $`\mathrm{\Omega }_f/\omega 0.75`$. Between these two values lies the threshold for the dynamic instability of the condensate which changes the subsequent dynamics dramatically Subhasis . In the first group, as the rotation frequency gradually increases with time, the condensate adiabatically follows a steady state, apart from excitations of the surface modes leading to a very small periodic variation of the angular momentum (see curve for $`\mathrm{\Omega }_f=0.7\omega `$ in Fig.2a). With increasing $`\mathrm{\Omega }_f`$, the condensate’s final state becomes more and more elliptically deformed, surrounded by a ring of vortices which however never enter it. The second group shows completely different behaviour when $`\mathrm{\Omega }(t)/\omega 0.75`$ (see left column of Fig.1): the instability sets in, the condensate becomes slightly S-shaped at $`t450\omega ^1`$ before being highly deformed and undergoing very turbulent motion Feder . This is accompanied by a large increase in angular momentum of the gas from almost zero when $`\mathrm{\Omega }(t)<0.75\omega `$ to between $`5\mathrm{}`$-$`7\mathrm{}`$ per particle (see Fig.2a). At this point ($`t670\omega ^1`$) several vortices enter the high density region and, in less than 200$`\omega ^1`$, settle down to form a well-defined lattice. After this, a period of relaxation of around 800$`\omega ^1`$ begins with the initially rotating lattice finally stopping in the rotating frame. There remains a small random motion of the vortices around their equilibrium positions in the lattice together with density fluctuations in and around the condensate.

At the end of the simulation, damping of the vortex motion has occurred and the initial energy of the vortex motion has been transferred in an effectively irreversible way to other degrees of freedom of the field. A similar phenomenon has been observed for the relative motion of two condensates Castin . If we assume that the field has reached a thermal distribution, we can calculate the temperature of the system by taking the final state of the simulation and evolving it with the conjugate gradient method in a trap rotating at $`\mathrm{\Omega }_f`$. This reduces its energy and takes it to the local minimum associated with the vortex lattice. We then calculate the energy difference $`\mathrm{\Delta }E`$ between the final state of the simulation and the one at the minimum. Assuming that Bogoliubov theory is valid, $`\mathrm{\Delta }E`$ must correspond to the energy of a classical thermal distribution of weakly coupled harmonic oscillators of amplitude $`b_n`$ which obeys the equipartition formula $`b_n^{}b_nϵ_n=k_BT`$, with $`n`$ being the Bogoliubov mode number. So, if $`𝒩`$ is the number of modes in the system (and keeping in mind that we have to subtract the one corresponding to the condensate) then we have

$$\mathrm{\Delta }E=\underset{n}{}b_n^{}b_nϵ_n=(𝒩1)k_BT.$$
(3)
The final temperature is $`0.616\mathrm{}\omega `$ for $`\mathrm{\Omega }_f=0.75\omega `$ and $`0.754\mathrm{}\omega `$ for $`\mathrm{\Omega }_f=0.8\omega `$, in other words it is extremely small, less than a tenth of the chemical potential.

We have also carried out a simulation on a larger grid ($`64\times 64\times 256`$) to check the dependence on size. We chose $`\mathrm{\Omega }_f=0.8\omega `$ and compared it with the one on the $`32\times 32\times 128`$ grid. The vortex nucleation and crystallisation phases are very similar and occur at roughly the same times. At longer times two differences arise: first, there are large underdamped oscillations of the angular momentum (see Fig.2a). An analysis of the simulation suggests that these oscillations are those of the scissors mode. Second, the final temperature ($`0.094\mathrm{}\omega `$) differs by the ratio of the number of modes as expected: at time $`t=500\omega ^1`$ when $`\mathrm{\Omega }(t)=\mathrm{\Omega }_f`$, $`\psi `$ had not yet reached the boundary in the smaller grid case and so the evolution of $`\psi `$ on both grids was identical up to this time with the same total energy which was conserved at later times resulting in the same value of $`\mathrm{\Delta }E`$. This exemplifies the fact that, in classical field theories, the relationship between energy and temperature depends on the energy cutoff.

Since the thermal occupation of the modes is directly proportional to the temperature, we expect that all relaxation processes which involve scattering from or into those modes (such as Landau-Beliaev damping) will be reduced. We are thus led to the conclusion that, for our simulations starting at $`T=0`$, relaxation rates in the equilibration period after the formation of the lattice could depend on the size of the grid. However, with the present numerical results, we were not able to demonstrate this.

Finite initial temperature: We performed simulations starting with $`k_BT=4\mathrm{}\omega `$ and $`k_BT=8\mathrm{}\omega `$. Now not only the condensate but also other modes are occupied in the initial state, with a thermal distribution. For a final rotation frequency below that of the dynamic instability, the situation is quite different from that of the zero temperature case: the condensate is never deformed and the vortices do enter the condensate if $`\mathrm{\Omega }_f0.55\omega `$ for $`k_BT=4\mathrm{}\omega `$ and if $`\mathrm{\Omega }_f0.5\omega `$ for $`k_BT=8\mathrm{}\omega `$. We have checked numerically that $`0.51\omega `$ is the Landau critical frequency above which the vortex-free condensate is no longer a local minimum of the energy. During the real time evolution of the right column of Fig.1 corresponding to $`\mathrm{\Omega }_f=0.6\omega `$, we find that the vortices enter only one at a time. That is, as the angular momentum of the cloud increases, one vortex, out of the group of vortices that surrounds the condensate will enter it and spiral slowly clockwise towards the center on a time scale of hundreds of $`\omega ^1`$. After that vortex has reached the center, a second one enters slowly, repeating the trajectory of the first until it starts to interact with it and the two orbit around each other for a while after which a third will enter. At the end of the simulation, coinciding with the achievement of the plateau in angular momentum, the lattice becomes stationary in the rotating frame and no further vortex enters the condensate (see right column of Fig.1). For $`\mathrm{\Omega }_f=0.7\omega `$ we find that the condensate deforms itself elliptically after which three vortices enter at the same time and form a rotating lattice. After that, and spaced by several hundred $`\omega ^1`$, a fourth and then a fifth vortex enter. Finally, two further vortices enter simultaneously to form the final seven vortex lattice. At each intermediate stage there is always a well defined lattice present although it is not stationary in the rotating frame. We should contrast this with the scenario of Ueda ; Gardiner  where a large number of vortices enter all at once into the condensate in a ring configuration and then some of them form a lattice while others are shed and leave the condensate.

For $`\mathrm{\Omega }_f`$ above the dynamic instability frequency the situation is quite similar to the corresponding one at $`T=0`$. Once the instability has set in the lattice is formed for both temperatures in about $`200\omega ^1`$ as in the $`T=0`$ case (see Fig.2b,c-this weak temperature dependence was also found experimentally by Raman ). The main difference is that the relaxation time for the lattice to stop rotating is much shorter, on the order of a hundred $`\omega ^1`$, not eight hundred.

It is important to emphasize the multimode interpretation of the field. Transposing Penrose and Onsager’s definition to the classical field theory, the condensate wavefunction is defined as the eigenvector corresponding to the largest eigenvalue of the one-body density matrix $`\psi ^{}(𝐫^{})\psi (𝐫)`$ where the average is over an ensemble of initial states. If the system becomes turbulent because it encounters an instability, the trajectories of the neighboring realisations will diverge exponentially. However, after averaging, we believe that the condensate wavefunction will not be a turbulent function. For $`T=0`$ there is only one initial state and so we replace ensemble averaging by one over time in the steady state regime LL . In our simulations with $`\mathrm{\Omega }_f=0.8\omega `$ the system must therefore be understood as becoming intrinsically multimode even though we started at $`T=0`$ with a pure condensate.

Conclusions: We have identified two very different regimes for the crystallisation of the vortex lattice in the classical field theory. At $`T=0`$ when the dynamic instability sets in (for high enough $`\mathrm{\Omega }_f`$), the vortices enter the condensate abruptly and settle into a lattice even if there is no dissipative term in the equation of motion. This scenario is consistent with the one found experimentally at the ENS ENS2  with comparable timescales. At finite temperatures, where Ueda ; Gardiner  add an explicit damping term, we see that the thermal classical field can provide the damping of the vortex motion on its own which leads to very different behaviour from that observed by those authors in that the vortices can enter one by one into the condensate and settle into a lattice before the entry of the following one. So far there has been no experimental check of this scenario. Finally, as the $`T=0`$ case shows, any theoretical model which singles out the condensate mode for separate treatment with a Gross-Pitaevskii-type equation could run into trouble in turbulent situations since the separation between condensed and non-condensed modes would be hard to keep.

We have been informed that crystallisation of the vortex lattice has also been observed in a simulation without a damping term by the group of N. Bigelow. We thank B. Durin, L. Carr, I. Carusotto, G. Shlyapnikov and J. Dalibard for useful contributions. C. L. acknowledges a fellowship from the Fundação para a Ciência e Tecnologia of Portugal. LKB is a unit of ENS and of Université Paris 6 associated to CNRS.


text:
# On reducing Terrorism Power: A Hint from Physics

## 1 Acknowledgments

We would like to thank Dietrich Stauffer for very helpful comment on making the manuscript clearer.

## 2 References

1. S. Galam, Eur. Phys. J. B 26, Rapid Note, 269-272 (2002).
2. D. Stauffer and A. Aharony, Introduction to Percolation theory, Taylor and Francis, London (1994).
3. M. Barth l my and L. A. N. Amaral, Phys. Rev. Lett. 82, 3180-3183 (1999).
4. D. Helbing, I. Farkas and T. Vicsek, Nature 407, 487-490 (2000).
5. S. Solomon, G. Weisbuch, L. de Arcangelis, N. Jan and D. Stauffer, Physica A277 (1-2), 239-247 (2000).
6. F. Lilieros, C. R. Edling, L. A. Nunes Amaral, H. E. Stanley and Y. Aberg, Nature 411, 907-908 (2001).
7. V. E. Krebs, Connections 24(3), 43-52 (2002).
8. S. Galam and A. Mauger, Phys. Rev. E53, 2177 (1996); Phys. Rev. E56, 1 (1997).
9. D. Stauffer and R. M. Ziff, Int. J. Mod. Phys. C 11, 205 (2000); A. Ticona and D. Stauffer, Physica A 290, 1 (2001).

## 3 Figure captions

Fig. 1: A schematic view of 64 persons on a square earth surface with 20 passive supporters in black circles, others 44 being in gray circle. People are aggregated into plaquettes of 4 nearest neighbors as shown for the upper left 16 ones.
Fig. 2: Nearest neighbor bonds on the earth surface before the establishment of the representational social space
Fig. 3: Construction of the third dimension associated to the 8 square lattice sites of Figure (2). For each plaquette of 4 sites, the associated vertical line contains 4 sites to rank the averaged degree of identification to the terrorist goal.
Fig. 4: The representational social space (upper part) with only passive supporters occupying their respective sites. Other sites are empty. Below is the two-dimensional earth space.
Fig. 5: Earth configuration after the establishment of the representational social space. Two more additional links connecting respectively passive supporters (a)-(f) and (f)-(h) have been created. Their respective range goes beyond the original nearest neighbor range.
Fig. 6: Variations of the percolation thresholds as function of the connectivity q from q=5 to q=30 at respectively d=2, d=4 and d=10. The horizontal line pc=0.1 shows when percolation occurs for density p=0.10.
Fig. 7: The three-dimensional plane of the percolation threshold as function of both q and d for respectively q=5 to q=20 and d=2 to d=10. The intersection with the plane at pc=0.10 is shown. Above the plane, terrorism is localized geographically while below it is long ranged.


text:
# Experimental Realization of Entanglement Concentration and A Quantum Repeater

## Acknowledgements

This work was supported by the National Natural Science Foundation of China, the Chinese Academy of Sciences and the National Fundamental Research Program (under Grant No. 2001CB309303).

* To whom correspondence should be addressed. E-mail:pan@ap.univie.ac.at

Figure Captions:

Figure captions:

Figure $`1`$: Schematic drawings showing the principles of entanglement concentration (a) and quantum repeater (b). a, The sources S emit two entangled photon pairs into the modes (1,2) and (3,4). The two pairs are originally in the same unknown non-maximally entangled state. One member of each pair is then sent to Alice, and the other one to Bob. After photon 4 passes through the half wave plate R90, Bob superposes the photons at the polarizing beam splitter PBS. By performing a polarization measurement on each of the photons in mode 3 and $`4^{}`$ in the $`+/`$ basis, Alice and Bob can project the photons 1 and $`2^{}`$ onto a maximally entangled state. b, Alice and Charlie share the pair (1,2) and, Bob and Charlie share the pair (3,4). Interfering the photons 2 and 4 at the PBS and performing a subsequent polarization measurement on each of the modes $`2^{}`$ and $`4^{}`$ in the $`+/`$ basis, Charlie can thus generate a maximally entangled state between Alice and Bob (see text).

Figure $`2`$: The experimental setups for entanglement concentration (a) and quantum repeater (b). a, A UV pulse passes through a BBO twice to generate two entangled photon pairs. By inserting one, two or four pieces of Brewster’s windows (BWs) into each of the modes 1 and 3, we can prepare three different nonmaximally entangled states. After the photons 2 and 4 pass through the PBS, conditional on a two-fold coincidence detection behind the two $`45^{}`$ polarizers Pol in the modes 3 and $`4^{}`$, the photons $`1`$ and $`2^{}`$ will be in a maximally entangled state. b, If the pairs (1,2) and (3,4) are in the same nonmaximally entangled state, then detecting a two-fold coincidence behind the two $`45^{}`$ polarizers Pol in the modes $`2^{}`$ and $`4^{}`$ is sufficient to project the photons 1 and 3 onto the desired maximally entangled state. Otherwise, additional half wave plates and BWs are needed in the modes 1 and 3 to generate the desired entangled state between Alice and Bob.

Figure $`3`$: Experimental results demonstrating the entanglement concentration. a, b, and c show the fractions of the nonmaximally entangled state in the $`H/V`$ basis, when inserting one, two and four BWs into the modes 1 and 3. d, e, and f show the corresponding fractions of the photons 1 and $`2^{}`$ in the $`H/V`$ basis after entanglement concentration. The data clearly confirm that our concentration protocol significantly improved the relative intensity between $`|H|V`$ and $`|H|V`$ components.

Figure $`4`$: Experimental results confirming the entanglement between the photons 1 and $`2^{}`$. The two-fold coincidence count rate $`|+|+`$ and $`||+`$ is shown as a function of the pump delay mirror position under the condition that a $`|+|+`$ coincidence is registered in the modes $`3`$ and $`4^{}`$.


text:
# Mean-field dynamical density functional theory

## References


text:
# Dynamic structure factor of the Ising model with purely relaxational dynamics.

## I Introduction

The dynamic structure factor $`C(k,\omega )`$ is a physically interesting quantity that can be directly measured in scattering experiments. Indeed, in neutron-scattering experiments and in Born approximation, $`C(k,\omega )`$ is proportional to the cross section for inelastic scattering with momentum transfer $`k`$ and energy transfer $`\omega `$. At a continuous phase transition the structure factor shows a universal scaling behavior that depends on the dynamic universality class of the system. In this paper we consider the dynamic universality class of the three-dimensional Ising model with purely relaxational dynamics without conservation laws, which is also known as model A, see Ref. . This dynamic universality class should be appropriate to describe the dynamic critical properties of uniaxial magnetic systems in which the energy is not conserved due to the coupling to phonons and of alloys such as $`\beta `$-brass at the order-disorder transition; see, e.g., Ref. . Note that this universality class does not describe the dynamic behavior of simple fluids and mixtures at the liquid-vapor or mixing transitions because of additional conservation laws . The model-A dynamics for the Ising universality class may also be relevant for the dynamics of quarks and gluons at finite temperature and finite baryon-number chemical potential $`\mu `$. Indeed, using quantum chromodynamics, which is the current theory of strong interactions, one can argue that in the $`T\mu `$ plane there exists an Ising-like continuous transition at the endpoint of a first-order phase transition line . Model A (or model C if one takes into account the baryon-number conservation ) should describe the critical dynamics at this critical point (see also Ref. ).

In this paper we compute the structure factor $`C(k,\omega )`$ for the three-dimensional Ising universality class with purely relaxational dynamics (model A) in equilibrium. We consider the $`ϵ`$-expansion and compute $`C(k,\omega )`$ to two loops in the high-temperature phase and to one loop in the whole temperature magnetic-field plane. In the high-temperature phase we also perform a Monte Carlo simulation, using the standard Ising model and the Metropolis dynamics. We find that, for moderately large $`k`$ and $`\omega `$, $`C(k,\omega )`$ is very well approximated by its mean-field (Gaussian) expression. In the high-temperature phase the field-theoretical analysis and the simulation show that corrections to the mean-field behavior are less than 1% for $`k\xi 5`$ and $`\omega \tau 10`$, where $`\xi `$ is the correlation length and $`\tau `$ is the zero-momentum autocorrelation time. In the low-temperature phase, on the basis of a one-loop field-theoretical analysis, we expect slightly larger corrections, but still of the order of a few percent. Note that our study concerns the scaling behavior of $`C(k,\omega )`$ in equilibrium, but it should be observed that similar conclusions have been obtained for the nonequilibrium dynamics in which one quenches a disordered system at $`T_c`$ .

The paper is organized as follows. In Sec. II we define the quantities that are computed in the following sections. We report a list of definitions together with some properties that are used in the calculation. In Sec. III we present our field-theoretical results, obtained using the general formalism of Refs. . Sec. IV is devoted to the presentation of the Monte Carlo results. In the Appendix we report some technical details.

## II Definitions and basic observables

In this paper we consider the equilibrium dynamics for an Ising-like theory with scalar order parameter $`\phi (r,t)`$ at temperature $`T`$ in the presence of a time- and space-independent external (magnetic) field $`H`$. We consider the connected two-point correlation function of the order parameter

$$G(r,t_1t_2)\phi (r,t_1)\phi (0,t_2)_{\mathrm{conn}},$$
(1)
where we have assumed to be in equilibrium, so that the correlation function depends only on the difference $`t_1t_2`$. Then, we define its Fourier transform $`\stackrel{~}{G}(k,t)`$ with respect to $`r`$,

$$\stackrel{~}{G}(k,t)=d^dre^{ikr}G(r,t),$$
(2)
and the structure factor $`C(k,\omega )`$

$$C(k,\omega )=𝑑te^{i\omega t}\stackrel{~}{G}(k,t).$$
(3)
Here, we do not write explicitly the dependence on $`T`$ and $`H`$ that is always understood in the notation. Near the critical point correlations develop both in space and time. They can be characterized in terms of the second-moment correlation length $`\xi `$ and of the zero-momentum integrated autocorrelation time $`\tau `$ defined by

$`\xi ^2`$ $``$ $`{\displaystyle \frac{1}{2d\chi }}{\displaystyle d^dr|r|^2G(r,0)}={\displaystyle \frac{1}{\chi }}{\displaystyle \frac{\stackrel{~}{G}(k,0)}{k^2}}|_{k^2=0},`$ (4)
$`\tau `$ $``$ $`{\displaystyle \frac{1}{2\chi }}{\displaystyle 𝑑t\stackrel{~}{G}(0,t)}={\displaystyle \frac{1}{2\chi }}C(0,0),`$ (5)
where $`\chi \stackrel{~}{G}(0,0)`$ is the static magnetic susceptibility. As is well known, for $`TT_c`$ ($`T_c`$ is the critical temperature) and $`H0`$, $`\xi `$ and $`\tau `$ diverge. In the absence of a magnetic field,

$$\xi |TT_c|^\nu ,\tau |TT_c|^{z\nu }\xi ^z,$$
(6)
where $`\nu `$ is the usual static exponent and $`z`$ a dynamic exponent that depends on the considered dynamics. The static exponents for the three-dimensional Ising universality class are very well determined , see Ref.  for a review. Present-day lattice studies give estimates that can be summarized as follows : $`\gamma =1.2372(5)`$, $`\nu =0.6301(4)`$, $`\eta =0.0364(5)`$, $`\alpha =0.110(1)`$. The exponent $`z`$ depends on the dynamics. For model-A dynamics, estimates of $`z`$ in three dimensions have been obtained by employing several methods. There exist field-theoretical perturbative calculations in different schemes  and Monte Carlo analyses that determine $`z`$ by studying the equilibrium dynamics at $`T_c`$ in finite volume , damage spreading , the critical relaxation from an ordered state , hysteresis scaling , and the short-time critical dynamics . For experimental determinations see, e.g., Refs. . The exponent $`z`$ turns out to be slightly larger than two. For instance, $`z2.017`$ from the fixed-dimension field-theoretical expansion , $`z2.02`$ from an analysis interpolating the $`4ϵ`$ and $`1+ϵ`$ expansions , and $`z2.04`$ from Monte Carlo simulations.

Near the critical point correlation functions show a scaling behavior. For the static structure factor, neglecting scaling corrections, we have

$$\frac{\chi }{\stackrel{~}{G}(k,0)}g_{\mathrm{stat}}(y;x),$$
(7)
where $`xa_0(TT_c)M^{1/\beta }`$, $`yk^2\xi ^2`$, $`M\phi `$ is the time-independent (we only consider the equilibrium dynamics) static magnetization, and $`a_0`$ is a normalization factor that is fixed by requiring that $`x=1`$ corresponds to the coexistence line. The magnetization $`M`$ is related to $`T`$ and $`H`$ by the equation of state, which, in the critical limit, can be written in the scaling form

$$H=b_0M^\delta f(x),$$
(8)
where $`b_0`$ is a nonuniversal constant which is fixed by the condition $`f(0)=1`$.

The function $`g_{\mathrm{stat}}(y;x)`$ has been extensively studied, both in the high-temperature  and in the low-temperature phase ; see Ref.  for an extensive review. In the high-temperature phase, the static function $`g_{\mathrm{stat}}^+(y)`$ is known to $`O(ϵ^3)`$ , and satisfies $`g_{\mathrm{stat}}^+(y)=1+y+O(ϵ^2y^2)`$. Its small-momentum expansion in three dimensions has been accurately determined using high-temperature expansion techniques, see, e.g., Refs. , finding

$$g_{\mathrm{stat}}^+(y)=1+y0.000390(6)y^2+0.0000088(1)y^3+O(y^4).$$
(9)
There are also precise estimates of the equation-of-state scaling function $`f(x)`$ .

Eq. (7) can be extended to finite values of $`t`$. In the critical limit we can write

$$\frac{\chi }{\stackrel{~}{G}(k,t)}g(y,s;x)$$
(10)
with $`st/\tau `$. We can also define a scaling function for the structure factor:

$$\frac{C(k,\omega )}{2\tau \chi }𝒞(y,w;x),$$
(11)
where $`w\omega \tau `$ and

$$𝒞(y,w;x)=\frac{1}{2}𝑑se^{iws}[g(y,s;x)]^1.$$
(12)
We also define an integrated autocorrelation time at momentum $`k`$,

$$\tau (k)\frac{1}{2}𝑑t\frac{\stackrel{~}{G}(k,t)}{\stackrel{~}{G}(k,0)}=\frac{C(k,0)}{2\stackrel{~}{G}(k,0)},$$
(13)
and an exponential autocorrelation time

$$\tau _{\mathrm{exp}}(k)\underset{|t|\mathrm{}}{lim}\frac{|t|}{\mathrm{ln}\stackrel{~}{G}(k,t)},$$
(14)
which controls the large-$`t`$ behavior of $`\stackrel{~}{G}(k,t)`$: $`\stackrel{~}{G}(k,t)\mathrm{exp}[|t|/\tau _{\mathrm{exp}}(k)]`$ for $`|t|\mathrm{}`$. In the scaling limit, neglecting scaling corrections,

$`{\displaystyle \frac{\tau (k)}{\tau }}`$ $``$ $`𝒯(y;x)=𝒞(y,0;x)g_{\mathrm{stat}}(y),`$ (15)
$`{\displaystyle \frac{\tau _{\mathrm{exp}}(k)}{\tau (k)}}`$ $``$ $`𝒯_{\mathrm{exp}}(y;x)={\displaystyle \frac{1}{|w_0(y;x)|𝒞(y,0;x)g_{\mathrm{stat}}(y)}},`$ (16)
where $`\pm iw_0(y;x)`$ are the zeros of $`[𝒞(y,w;x)]^1`$ at fixed $`y`$ and $`x`$ on the imaginary $`w`$-axis that are nearest to the origin $`w=0`$.

For a Gaussian free theory we have

$$C(k,\omega )|_{\mathrm{Gaussian}}=\frac{2\mathrm{\Omega }\chi m^2}{\mathrm{\Omega }^2(m^2+k^2)^2+\omega ^2},$$
(17)
where $`\mathrm{\Omega }`$ is an Onsager transport coefficient and $`m1/\xi `$. It follows

$`[𝒞(y,w;x)]^1=(1+y)^2+w^2,`$ (18)
$`[𝒯(y;x)]^1=1+y,`$ (19)
$`𝒯_{\mathrm{exp}}(y;x)=1.`$ (20)
For $`y0`$ and $`w0`$ the above-defined scaling functions have a regular behavior and one can write

$`[𝒞(y,w;x)]^1`$ $`=`$ $`(1+y)^2+w^2+{\displaystyle \underset{m,n=0}{}}c_{n,m}(x)y^nw^{2m},`$ (21)
$`[𝒯(y;x)]^1`$ $`=`$ $`1+y+{\displaystyle \underset{n=0}{}}t_n(x)y^n,`$ (22)
$`𝒯_{\mathrm{exp}}(y;x)`$ $`=`$ $`1+{\displaystyle \underset{n=0}{}}t_{\mathrm{exp},n}(x)y^n,`$ (23)
with $`c_{0,0}(x)=0`$ because of the definition of $`\tau `$. The expansion coefficients $`c_{n,m}(x)`$, $`t_n(x)`$, and $`t_{\mathrm{exp},n}(x)`$ parametrize the deviations from the Gaussian behavior (20) in the low-frequency and low-momentum regime.

At the critical point, $`T=T_c`$, $`H=0`$, the structure factor obeys the scaling law

$$C(k,\omega )=\frac{1}{\omega ^{(2\eta +z)/z}}f_C(\omega k^z),$$
(24)
with $`f_C(\mathrm{})`$ finite, which implies that, for $`y\mathrm{}`$, $`w\mathrm{}`$ keeping $`uwy^{z/2}=u_0\omega k^z`$ fixed, we have

$$𝒞(y,w;x)=\frac{f_0}{w^{(2\eta +z)/z}}f_C(u/u_0),$$
(25)
where $`f_0`$ is a normalization constant.

For large $`w`$ at $`y`$ and $`x`$ fixed we have

$$𝒞(y,w;x)c_{\mathrm{}}(y;x)w^{(2\eta +z)/z},$$
(26)
where $`c_{\mathrm{}}(y;x)`$ is finite and $`x`$-independent for $`y\mathrm{}`$. The large-frequency behavior of the structure factor allows us to compute the nonanalytic small-$`s`$ behavior of $`g(y,s;x)`$ at $`y`$ and $`x`$ fixed. We obtain for $`s0^+`$

$$[g(y,s;x)]_{\mathrm{nonanalytic}}^1=g_0(y;x)s^{(2\eta )/z},$$
(27)
where

$$g_0(y;x)=\frac{1}{\pi }c_{\mathrm{}}(y;x)\mathrm{sin}(\pi \rho /2)\mathrm{\Gamma }(1\rho ),$$
(28)
with $`\rho 1+(2\eta )/z`$. Notice that, since $`(2\eta )/z0.96`$, the nonanalytic small-$`t`$ behavior of $`\stackrel{~}{G}(k,t)`$ turns out to be practically indistinguishable from the analytic background.

## III Field-theoretical results

### A Field-theoretical approach

In order to determine the critical behavior of a purely relaxational dynamics without conservation laws, the so-called model-A dynamics, one may start from the the stochastic Langevin equation

$$\frac{\phi (r,t)}{t}=\mathrm{\Omega }\frac{\delta (\phi )}{\delta \phi (r,t)}+\rho (r,t),$$
(29)
where $`\phi (r,t)`$ is the order parameter, $`(\phi )`$ is the Landau-Ginzburg-Wilson Hamiltonian

$$(\phi )=d^dx\left[\frac{1}{2}(\phi )^2+\frac{1}{2}r\phi ^2+\frac{1}{4!}u\phi ^4H\phi \right],$$
(30)
$`\mathrm{\Omega }`$ a transport coefficient (cf. Eq. (17)), and $`\rho (t)`$ a Gaussian white noise with correlations

$$\rho (r,t)=0,\rho (r_1,t_1)\rho (r_2,t_2)=\mathrm{\Omega }\delta (r_1r_2)\delta (t_1t_2).$$
(31)
The correlation functions generated by the Langevin equation (29) and averaged over the noise $`\rho `$, can be obtained starting from the field-theoretical action

$$S(\phi ,\widehat{\phi })=𝑑td^dx\left[\widehat{\phi }\frac{\phi }{t}+\mathrm{\Omega }\widehat{\phi }\frac{\delta (\phi )}{\delta \phi }\mathrm{\Omega }\widehat{\phi }^{\mathrm{\hspace{0.17em}2}}\mathrm{ln}J(\phi )\right].$$
(32)
The last term in the action is an appropriate Jacobian term that compensates the contributions of self-loops of response propagators .

In order to perform the field-theoretical calculation it is useful to introduce the response function $`Y(r,t)`$—it gives the linear response to an external magnetic field—defined by

$$Y(r,t_1t_2)=\widehat{\phi }(r,t_1)\phi (0,t_2),$$
(33)
(again we have assumed to be in equilibrium so that time-translation invariance holds), its Fourier transform $`\stackrel{~}{Y}(k,t)`$ with respect to $`r`$, and its double Fourier transform $`R(k,\omega )`$ with respect to $`r`$ and $`t`$, defined as $`C(k,\omega )`$ in Eq. (3). The response function and the two-point correlation function are strictly related. First, the zero-frequency response functions are related to the static correlation functions,

$$\stackrel{~}{G}(k,0)=\mathrm{\Omega }R(k,0).$$
(34)
Moreover, because of the fluctuation-dissipation theorem that holds for the equilibrium dynamics, we have

$$\omega C(k,\omega )=2\mathrm{\Omega }\mathrm{Im}R(k,\omega ).$$
(35)
Also the response function $`R(k,\omega )`$ shows a scaling behavior and one can write

$$\frac{\chi }{\mathrm{\Omega }R(k,\omega )}r(y,w;x),$$
(36)
neglecting scaling corrections. The function $`r(y,w;x)`$ is such that

$`r(y,0;x)=1+y+O(y^2),`$ (37)
$`r(0,w;x)=1iw+O(w^2),`$ (38)
$`[r(y,w;x)]^{}=r(y,w;x).`$ (39)
Then, it is easy to show by using Eqs. (34) and (35) that

$`r(y,0;x)`$ $`=`$ $`g_{\mathrm{stat}}(y;x),`$ (40)
$`𝒞(y,w;x)`$ $`=`$ $`{\displaystyle \frac{\mathrm{Im}r(y,w;x)}{w|r(y,w;x)|^2}}.`$ (41)
For a Gaussian theory

$$r(y,w;x)=1+yiw.$$
(42)
The behavior of $`r(y,w;x)`$ for small $`w`$ and large $`w`$ is similar to that of $`𝒞(y,w;x)`$. For small frequencies and momenta, the scaling function has a regular expansion in powers of $`w`$ and $`y`$:

$$r(y,w;x)=g_{\mathrm{stat}}(y;x)iw\left[1+\underset{n,m}{}r_{n,m}(x)(iw)^my^n\right],$$
(43)
where the coefficients $`r_{n,m}(x)`$ are real and parametrize the $`w`$-dependent deviations from the Gaussian behavior (42). For $`w\mathrm{}`$ at fixed $`y`$ we have

$$r(y,w;x)r_{\mathrm{}}^+(y;x)(iw)^{(2\eta )/z}.$$
(44)
### B Correlation functions in the disordered phase

In this Section we consider the equilibrium dynamics in the high-temperature phase $`H=0`$, $`T>T_c`$. In order to determine the two-point correlation function, we have computed the scaling function $`r^+(y,w)`$ (here and in the following we will not indicate $`x`$ and add instead a superscript “+” to remind the reader that we refer to the high-temperature phase) and we have then used Eq. (41).

A two-loop calculation in the framework of the $`ϵ`$ expansion gives

$$r^+(y,w)=g_{\mathrm{stat}}^+(y)iw\left[1+ϵ^2A(y,w)+O(ϵ^3)\right],$$
(45)
where $`A(y,w)`$ is reported in App. A. Note that $`A(0,0)=0`$ and $`A(y,w)^{}=A(y,w)`$, as expected from Eq. (39). The static function $`g_{\mathrm{stat}}^+(y)`$ is known to $`O(ϵ^3)`$ , and at order $`ϵ^2`$ it reads

$$g_{\mathrm{stat}}^+(y)=1+y+10^3ϵ^2[3.76012y^2+0.095966y^30.00407101y^4+O(y^5)]+O(ϵ^3).$$
(46)
Expanding $`A(y,w)`$ in powers of $`y`$ and $`w`$ one obtains the coefficients $`r_{n,m}^+`$ defined in Eq. (43). We have $`r_{n,m}^+=ϵ^2\overline{r}_{n,m}^+`$, where $`\overline{r}_{n,m}^+`$ are reported in Table I for $`n,m3`$. The coefficients $`\overline{r}_{n,m}^+`$ are rather small, the largest ones being of order $`10^3`$, and decrease quite rapidly. The analysis of the coefficients of the expansion of $`A(k^2,\omega )`$ in powers of $`k^2`$ (at fixed $`\omega `$) shows the presence of a singularity for $`w=3i`$. Therefore, we expect asymptotically

$$\overline{r}_{n,m}^+\frac{1}{3}\overline{r}_{n,m1}^+.$$
(47)
We have verified numerically this relation, although quantitative agreement is observed only for quite large values of $`m`$: for $`n=0`$, this relation is satisfied at the 10% level only for $`m41`$. Analogously, the coefficients of the expansion of $`A(k^2,\omega )`$ in powers of $`\omega `$ become singular for $`k^2=9`$, so that asymptotically

$$\overline{r}_{n,m}^+\frac{1}{9}\overline{r}_{n1,m}^+.$$
(48)
The behaviors (47) and (48) can be interpreted in terms of the analytic structure of $`R^+(k,\omega )`$. If one considers the structure factor, it is well known  that the nearest singularity  appearing in $`[\stackrel{~}{G}(k,0)]^1`$ is the three-particle cut at $`k=\pm 3im_{\mathrm{exp}}`$, where $`m_{\mathrm{exp}}`$ is the mass gap of the theory. Since in the critical limit $`m_{\mathrm{exp}}\xi 1`$  with very small corrections (more precisely $`m_{\mathrm{exp}}\xi 1=2.00(3)10^4`$, see Ref. ), the nearest singularity to the origin appearing in $`g_{\mathrm{stat}}^+(y)`$ corresponds to $`y9`$. In view of relation (48), it is natural to conjecture that the same behavior holds for $`R^+(k,\omega )`$, so that Eq. (48) should approximately hold for the three-dimensional coefficients $`r_{n,m}^+`$ and not only for their two-loop approximation.

Relation (47) is consistent with the idea that the three-particle cut also controls the small-$`w`$ behavior. In this case it is natural to conjecture that the coefficients of the expansion of $`[R^+(k,\omega )]^1`$ in powers of $`k^2`$ have a singularity for $`\omega =3i/\tau _{\mathrm{exp}}(0)`$. Thus, turning to the scaling function $`r^+(y;w)`$, we expect a singularity at $`w=3i\tau /\tau _{\mathrm{exp}}(0)3i`$, since, as we shall see, in the critical limit $`\tau /\tau _{\mathrm{exp}}(0)1`$. Therefore, we expect relation (47) to be a general property of the three-dimensional coefficients $`r_{n,m}^+`$.

This discussion indicates that $`𝒞(y,0)`$ and its $`w`$-derivatives at $`w=0`$ have a convergent expansion in $`y`$ for $`|y|9`$ and analogously that $`𝒞(0,w)`$ and its $`y`$-derivatives at $`y=0`$ have a convergent expansion for $`|w|3`$. Mathematically, this does not tell us much about the convergence of the double expansion which requires to know the singularity structure for both $`y,w0`$. At two loops, one can easily verify from the exact expression that $`𝒞(y,w)`$ has a convergent double expansion in the whole region $`|w|<3`$, $`|y|<9`$, and it is sensible to conjecture that the same is true for the exact expansion. From the results of Table I, one sees quite clearly that the response function $`R^+(k,\omega )`$ is well described by the Gaussian approximation for $`|w|3`$ and $`|y|9`$. Deviations should be smaller than 1% in this region. This result is very similar to that obtained for the static structure factor: in that case high-temperature expansions and Monte Carlo simulations  show that the deviations from the Gaussian behavior are less that 0.3% for $`y9`$.

We now consider the large-frequency behavior. At order $`ϵ^2`$ the function $`r_{\mathrm{}}^+(y)`$ defined in Eq. (44) turns out to be constant and given by

$$r_{\mathrm{}}^+(y)=1+0.00538992ϵ^2+O(ϵ^3).$$
(49)
Again the correction term is quite small.

Using the fluctuation-dissipation theorem (41), we obtain for the scaling function $`𝒞^+(y,w)`$:

$$[𝒞^+(y,w)]^1=g_{\mathrm{stat}}(y)^2+w^2+ϵ^2E(y,w)+O(ϵ^3),$$
(50)
where

$$E(y,w)=2w(1+y)\mathrm{Im}A(y,w)+\left[w^2(1+y)^2\right]\mathrm{Re}A(y,w).$$
(51)
We can then obtain the small-$`w`$ and small-$`y`$ behavior. For the coefficients $`c_{n,m}^+`$, see Eq. (23), we obtain $`c_{n,m}^+=ϵ^2\overline{c}_{n,m}^+`$, where the constants $`\overline{c}_{n,m}^+`$ are reported in Table II for $`n,m3`$.

Again, we should note that the coefficients $`\overline{c}_{n,m}^+`$ are very small and show the same pattern observed for $`\overline{r}_{n,m}^+`$. We expect that $`𝒞^+(y,w)`$ has singularities at $`y=9`$ and $`w=\pm 3i`$, so that $`|c_{n,m}/c_{n+1,m}||c_{n,m}/c_{n,m+1}|9`$. Thus, in complete analogy with what observed for the static structure factor and $`R^+(k,\omega )`$, the dynamic $`C^+(k,\omega )`$ is essentially Gaussian in the region $`y\mathrm{}<9`$ and $`|w|\mathrm{}<3`$.

We also compute the large-frequency behavior. For the coefficients $`c_{\mathrm{}}^+(y)`$ and $`g_0^+(y)`$, see Eqs. (26) and (27), we obtain

$`c_{\mathrm{}}^+(y)`$ $`=`$ $`10.00538992ϵ^2+O(ϵ^3),`$ (52)
$`2g_0^+(y)`$ $`=`$ $`1+0.00136716ϵ^2+O(ϵ^3),`$ (53)
where, since the corrections are very small, one may simply set $`ϵ=1`$ to obtain a three-dimensional numerical estimate. Therefore, for large $`w`$ we predict $`𝒞^+(y,w)0.995/w^{1.95}`$, which is not very different from the purely Gaussian behavior $`𝒞_{\mathrm{Gauss}}^+(y,w)1/w^2`$. Thus, the Gaussian approximation should be a reasonably good approximation even outside the small-$`w`$ region, $`w\mathrm{}<3`$, discussed above. Trusting the above estimate of $`c_{\mathrm{}}^+(y)`$ we find that $`𝒞^+(y,w)/𝒞_{\mathrm{Gauss}}^+(y,w)=1.12`$, 1.25, 1.41 respectively for $`w=10`$, 100, 1000. Thus, quite large values of $`w`$ are needed in order to observe a significant difference.

Finally we compute the scaling function $`𝒯_{\mathrm{exp}}^+(y)`$ defined in Eq. (16). For this purpose we need to compute $`\tau _{\mathrm{exp}}(k)`$ and therefore the large-$`t`$ behavior of $`\stackrel{~}{G}(k,t)`$. Because of the fluctuation-dissipation theorem, it is equivalent to consider $`\stackrel{~}{Y}(k,t)`$. For $`y<3`$ we obtain

$$\stackrel{~}{Y}(k,t)e^{s(1+y)}\{1ϵ^2sA[y,i(1+y)]\},$$
(54)
where $`st/\tau `$, while for $`y>3`$ we have

$$\stackrel{~}{Y}(k,t)e^{s(1+y)}+ϵ^2\frac{27}{8}\frac{e^{s(y/3+3)}}{s^4(y3)^2(y+9)^3}[1+O(s^{1/2})]+O(ϵ^3)$$
(55)
For $`y<3`$ the correction term exponentiates as expected, and as a consequence we obtain

$$𝒞^+(y,0)[g_{\mathrm{stat}}^+(y)]^2𝒯_{\mathrm{exp}}^+(y)=1+ϵ^2A[y,i(1+y)]+O(ϵ^3).$$
(56)
On the other hand, for $`y>3`$ the correction term decreases with a different exponential factor which dominates for large values of $`t`$, suggesting that, at leading order in $`ϵ`$, $`\tau _{\mathrm{exp}}(k)/\tau =(y+9)/3`$. In other words, the interaction turns on a new singularity (a three-particle cut) that becomes the leading one for $`y`$ large enough. However, this is not the end of the story. Indeed, by considering graphs in which one recursively replaces each line with a two-loop watermelon graph one obtains contributions to $`\stackrel{~}{Y}(k,t)`$ decreasing as $`\mathrm{exp}[s(3^ny+3^n)]`$ ($`3^n`$-particle cut) which would be more important for $`y`$ large enough. These singularities will not probably be the only ones, since we also expect a 5-particle cut, a 7-particle cut, etc. On the basis of these results, we expect $`𝒯_{\mathrm{exp}}^+(y)`$ to have several singularities on the positive real $`y`$ axis and to become eventually infinite as $`y\mathrm{}`$. This is not unexpected since, for $`y\mathrm{}`$, $`R^+(k,\omega )`$ behaves as $`\omega ^{(2\eta )/z}`$ and therefore has a branch cut starting at $`\omega =0`$.

For $`y<3`$, we can use Eq. (56) to compute the coefficients $`t_{\mathrm{exp},n}^+`$ defined in Eq. (23). We obtain, at order $`ϵ^2`$, $`t_{\mathrm{exp},0}^+=0.00110075ϵ^2`$, $`t_{\mathrm{exp},1}^+=0.00337789ϵ^2`$, $`t_{\mathrm{exp},2}^+=0.000217173ϵ^2`$, etc. The coefficients decrease as $`t_{\mathrm{exp},\mathrm{n}}^+/t_{\mathrm{exp},\mathrm{n}+1}^+3`$, which reflects the presence of a singularity at $`y=3`$. Again, for $`y<3`$ the deviations from a purely Gaussian behavior are very small.

### C Correlation function in the $`(t,H)`$ plane

In the presence of an external magnetic field $`H`$, a one-loop calculation gives

$`r(y,w;x)=g_{\mathrm{stat}}(y;x)iw\left[1+ϵ{\displaystyle \frac{2}{3+x}}B(y,w)+O(ϵ^2)\right],`$ (57)
where $`B(y,w)`$ is defined in App. A and

$`g_{\mathrm{stat}}(y;x)=1+y+{\displaystyle \frac{2ϵ}{3+x}}\left[1{\displaystyle \frac{y}{12}}+{\displaystyle \frac{\sqrt{4+y}}{2\sqrt{y}}}\mathrm{ln}{\displaystyle \frac{\sqrt{4+y}+\sqrt{y}}{\sqrt{4+y}\sqrt{y}}}\right]+O(ϵ^2).`$ (58)
Note that the $`O(ϵ)`$ correction vanishes for $`x\mathrm{}`$ in agreement with the results of the previous Section. Moreover, the $`x`$-dependence is very simple and in Eqs. (57) and (58) is always given by the prefactor $`2/(3+x)`$ that becomes 1 on the coexistence curve $`x=1`$. As a consequence, such a prefactor will always appear in this Section, multiplying the low-temperature results that will be specified by adding a superscript “$``$” to all definitions. Of course, such a simple $`x`$-dependence does not hold at higher loops, as it can be seen for instance from the two-loop results of Ref.  for the static structure factor.

One can easily derive the small-momentum and small-frequency behavior by expanding the function $`B(y,w)`$. The coefficients $`r_{n,m}(x)`$, see Eq. (43), are given by

$$r_{n,m}(x)=2ϵ\frac{\overline{r}_{n,m}^{}}{3+x}+O(ϵ^2),$$
(59)
where $`\overline{r}_{n,m}^{}`$ are given in Table III for $`m3`$ and $`n2`$.

Again, we note that the corrections to the Gaussian behavior are small, although a factor-of-ten larger than the corresponding high-temperature ones. For instance, $`\overline{r}_{0,1}^{}0.02`$ to be compared with $`\overline{r}_{0,1}^+0.002`$. Moreover, the coefficients decrease slower with $`n`$ and $`m`$. This fact can be understood in terms of the singularities of the function $`B(y,w)`$. A simple analysis shows the presence of singularities for $`w=2i`$ and $`y=4`$, so that asymptotically

$$\overline{r}_{n,m}^{}\frac{1}{2}\overline{r}_{n,m1}^{},\overline{r}_{n,m}^{}\frac{1}{4}\overline{r}_{n1,m}^{}.$$
(60)
This behavior can be understood on general grounds. Considering the static structure factor, it is known that the nearest singularity in the low-temperature phase is the two-particle cut $`k=\pm 2im_{\mathrm{exp}}`$, so that $`g_{\mathrm{stat}}^{}(y)`$ has a singularity for $`y=4(m_{\mathrm{exp}}\xi )^24`$, where we have used the fact that in the critical limit $`m_{\mathrm{exp}}\xi 1`$ (more precisely, $`m_{\mathrm{exp}}\xi 0.96(1)`$ ). As we did for the high-temperature phase we can thus conjecture that also the singularities of the dynamic functions are controlled by the two-particle cut. Therefore, we expect singularities for $`y=4(m_{\mathrm{exp}}\xi )^24`$ and $`w=2i\tau /\tau _{\mathrm{exp}}(0)2i`$, where we have used the fact that $`\tau /\tau _{\mathrm{exp}}(0)1`$, with corrections of order a few percent as discussed below, in the critical limit. Therefore, Eq. (60) should also approximately apply to the three-dimensional coefficients $`r_{n,m}^{}`$.

The above-reported discussion shows that in the region $`|y|4`$, $`|w|2`$ the response function can be reasonably approximated by a Gaussian form. Note however that, while in the high-temperature phase corrections are expected to be less than 1%, here deviations should be larger.

We have also studied the large-frequency behavior. The coefficient $`r_{\mathrm{}}(y;x)`$ turns out to be $`y`$-independent at one loop:

$$r_{\mathrm{}}(y;x)=1\frac{ϵ}{3+x}+O(ϵ^2).$$
(61)
Note that the correction is quite large and thus significant deviations for the Gaussian behavior should be observed as soon as $`w`$ is large.

Using the fluctuation-dissipation theorem, we can compute at one loop the scaling function $`𝒞(y,w;x)`$. For the small-$`w`$, small-$`y`$ coefficients, we obtain

$`c_{n,m}(x)={\displaystyle \frac{2ϵ}{3+x}}\overline{c}_{m,n}^{}+O(ϵ^2),`$ (62)
$`t_n(x)={\displaystyle \frac{2ϵ}{3+x}}\overline{t}_n^{}+O(ϵ^2).`$ (63)
The coefficients $`\overline{c}_{m,n}^{}`$ and $`\overline{t}_n^{}`$ are reported in Table III for $`n2`$ and $`m3`$.

We have also investigated the large-frequency behavior. It is very simple to show, using the above-reported formulas, that at this order $`c_{\mathrm{}}(y;x)=1/r_{\mathrm{}}(y;x)`$ and $`g_0(y;x)=c_{\mathrm{}}(y;x)/2`$.

Finally, we consider $`𝒯_{\mathrm{exp}}(y;x)`$. For this purpose we need to compute the large-$`t`$ behavior of $`\stackrel{~}{Y}(k,t)`$. We observe a behavior analogous to that observed in the high-temperature phase. For $`y<2`$,

$$𝒞(y,0;x)[g_{\mathrm{stat}}(y;x)]^2𝒯_{\mathrm{exp}}(y;x)=1+\frac{2ϵ}{3+x}B[y,i(1+y)]+O(ϵ^2),$$
(64)
while for $`y>2`$ the two-particle cut contribution dominates so that $`\tau _{\mathrm{exp}}(k)/\tau =2+y/2`$. The discussion reported in Sec. III B can be repeated also here. One can easily identify diagrams that decrease as $`\mathrm{exp}[s(2^ny+2^n)]`$, indicating that $`𝒯_{\mathrm{exp}}(y;x)`$ has an infinite number of singularities on the $`y`$ axis and that it diverges for $`y\mathrm{}`$. For small $`y`$, we can use Eq. (64) to compute the small-$`y`$ expansion coefficients $`t_{\mathrm{exp},n}(x)`$. We have

$$t_{\mathrm{exp},n}(x)=\frac{2ϵ}{3+x}\overline{t}_{\mathrm{exp},n}^{}+O(ϵ^2).$$
(65)
Numerical values are reported in Table III. Note that $`\tau _{\mathrm{exp}}(0)/\tau =1+t_{\mathrm{exp},0}^{}1+0.0284ϵ`$, and thus we expect this ratio to be 1 with corrections of order of a few percent.

## IV Monte Carlo results

We determine the dynamic structure factor $`C(k,\omega )`$ and the scaling function $`\stackrel{~}{G}(k,t)`$ in the high-temperature phase $`H=0`$, $`T>T_c`$ for small values of $`k`$—as we shall see, we are able to reach $`k10/\xi `$—by means of a large-scale Monte Carlo simulation. We consider the Ising model on a cubic lattice, i.e. the Hamiltonian

$$=\beta \underset{ij}{}\sigma _i\sigma _j,$$
(66)
where $`\beta 1/T`$, $`\sigma _i=\pm 1`$, and the summation is over all nearest-neighbor pairs $`ij`$. We measure the correlation function

$$\stackrel{~}{G}(k,t)=\frac{1}{3}\underset{x,y,z}{}(e^{iqx}+e^{iqy}+e^{iqz})\sigma _{0,0,0}(t=0)\sigma _{x,y,z}(t),$$
(67)
for four different values of $`L`$ and $`\beta `$: (a) $`L=64`$, $`\beta =0.215`$; (b) $`L=64`$, $`\beta =0.219`$; (c) $`L=128`$, $`\beta =0.2204`$; (d) $`L=128`$, $`\beta =0.221`$. Of course, in Eq. (67) $`q=2\pi n/L`$ where $`n`$ is an integer. For each $`\beta `$ and $`L`$ we first reached equilibrium by running 20000 Swendsen-Wang iterations, then we collected $`N_{\mathrm{it}}`$ iterations using the Metropolis algorithm . The results of the simulations are reported in Table IV. There we report the number of iterations $`N_{\mathrm{it}}`$, the second-moment correlation length $`\xi `$ (for the $`L=128`$ lattices we report the more precise results of Ref. ) and the autocorrelation time $`\tau `$. Note that all lattices have $`L/\xi 6`$, a condition that usually ensures that finite-size effects are reasonably small (for static quantities corrections are less than 1%).

The correlation length $`\xi `$ has been determined by using a discretized form of Eq. (4):

$$\xi ^2=\frac{\chi /F1}{4\mathrm{sin}^2(\pi /L)},$$
(68)
where $`F=\stackrel{~}{G}(k,0)`$ with $`k=(2\pi /L,0,0)`$. The integrated autocorrelation time $`\tau `$, and also the autocorrelation times $`\tau (k)`$ considered below, have been determined using the self-consistent method of Ref. :

$$\tau (k)=\frac{1}{2}+\underset{t=1}{\overset{M(k)}{}}\frac{\stackrel{~}{G}(k,t)}{\stackrel{~}{G}(k,0)},$$
(69)
where $`t`$ is the Monte Carlo time in sweeps and the cutoff $`M(k)`$ is chosen self-consistently so that $`6\tau (k)<M(k)6\tau (k)+1`$. Since $`\stackrel{~}{G}(k,t)`$ decays exponentially, this choice makes the systematic error due to the truncation small, keeping the statistical variance small at the same time; see Ref.  for a discussion.

First, we check that $`\tau \xi ^z|TT_c|^{z\nu }`$. Using the precise estimate $`\beta _c=0.22165459(10)`$ of Ref. , we obtain from a least-square fit $`z=2.10(2)`$ including all data and $`z=2.11(5)`$ discarding the estimate of $`\tau `$ for lattice (a). This result is in reasonable agreement with the estimates reported in Sec. II, if we take into account that we quote here only the statistical error. The systematic error due to corrections to scaling and to neglected finite-size effects is probably larger.

Then, we determine the correlation function $`\stackrel{~}{G}(k,t)`$. In Fig. 1 we report the function

$$f^+(y,s)\frac{g^+(y,0)}{g^+(y,s)}\frac{\stackrel{~}{G}(k,t)}{\stackrel{~}{G}(k,0)},$$
(70)
for three different values of $`yk^2\xi ^2`$, $`y=0,4,16`$, as computed from lattices (a), (b), and (c). We have not included the results for lattice (d), because they have much larger errors. In order to obtain $`\stackrel{~}{G}(k,t)`$ for a given $`k2\pi n/L`$, we have performed a linear interpolation, using two nearby values of $`k`$. First, we observe reasonable scaling: corrections due to the finite values of $`\xi `$ and $`L`$ are under control, although they increase as $`y`$ increases. For $`y=0`$ the results for the three different lattices agree within a few percent, while for larger values of $`y`$ we observe larger discrepancies. In particular, for $`y=4`$ and $`y=16`$, the estimates of $`f^+(y,s)`$ obtained from lattice (b) are always larger than those obtained from (a) and (c), the discrepancy being of order 20% when $`f^+(y,s)10^1`$ and 80% when $`f^+(y,s)10^2`$. These differences are probably finite-size effects, since (a) and (c) have $`L/\xi 10`$, while $`L/\xi 8`$ for (b).

It is also remarkable that the plot of $`\mathrm{ln}f^+(y,s)`$ is a straight line, indicating that $`f^+(y,s)`$ is quite precisely a pure exponential. No deviations can be observed in Fig. 1. Therefore,

$$\stackrel{~}{G}(k,t)\stackrel{~}{G}(k,0)\mathrm{exp}\left[t/\tau _{\mathrm{exp}}(k)\right],$$
(71)
within the precision of our results. This behavior appears to be well satisfied in the region that we can safely investigate, i.e. $`1/10t/\tau (k)4`$ and $`k\xi 5`$. Therefore, the dynamic structure factor is well approximated by a Lorentzian in the region of not too large frequencies, i.e. for $`\omega \tau (k)\mathrm{}<10`$.

Then, we consider the scaling function $`𝒯(y)`$ that encodes the $`k`$-dependence of $`\tau (k)`$. In Fig. 2 we report our numerical results. Again, we observe good scaling up to quite large values of $`y`$. In the figure, we also report the Gaussian prediction $`𝒯(y)=1/(1+y)`$. It can be seen that the Gaussian approximation describes very well the numerical data. This result should have been expected on the basis of the results of Sec. III where we showed that the deviations from a Gaussian behavior are very small in the small-$`y`$ regime $`y\mathrm{}<9`$, and should remain small even for larger $`y`$. For instance, using the data with largest $`y`$ reported in Fig. 2, we estimate $`𝒯(y)=0.0053(3)`$ for $`y=181`$, to be compared with the Gaussian prediction $`0.0055`$. Thus, in the range $`y200`$ the discrepancy should be at most 4-10%.

Finally, we consider the function $`𝒯_{\mathrm{exp}}^+(y)`$. In order to compute $`\tau _{\mathrm{exp}}(k)`$ we define an effective quantity

$$\tau _{\mathrm{eff}}(t;k)\left[\mathrm{ln}\frac{\stackrel{~}{G}(k,t+1)}{\stackrel{~}{G}(k,t)}\right]^1.$$
(72)
The exponential autocorrelation time $`\tau _{\mathrm{exp}}(k)`$ is obtained from $`\tau _{\mathrm{eff}}(t;k)`$ by letting $`t`$ go to infinity. In practice, we can only compute $`\tau _{\mathrm{eff}}(t;k)`$ up to $`t`$ of the order of (1-2)$`\times \tau (k)`$ since errors increase rapidly. In Fig. 3 we report the ratio $`\tau _{\mathrm{eff}}(t;k)/\tau (k)`$ for $`t=\tau (k)`$ for lattice (c) which is the only one that allows us to reach large values of $`y`$. We observe $`\tau _{\mathrm{eff}}(t;k)\tau (k)`$ within the precision of our results. It is tempting to conclude that $`\tau _{\mathrm{exp}}(k)\tau (k)`$ for $`y<50`$, but this is in contrast with the theoretical results of Sec. III B. Indeed, we showed there that $`\tau _{\mathrm{exp}}(k)\tau (k)`$ with very small corrections for $`y<3`$, but we noticed that this relation breaks down for larger values of $`y`$. For instance, for $`y>3`$, our two-loop calculation gives $`\tau _{\mathrm{exp}}(k)/\tau (k)=(3+3y)/(9+y)`$ which is significantly larger than 1 for $`y>3`$. As we already discussed this prediction should not be taken seriously, unless $`y`$ is close to 3, since other singularities should be present, and indeed we expect $`\tau _{\mathrm{exp}}(k)/\tau (k)`$ to diverge as $`k\mathrm{}`$. Therefore, our numerical data show that the asymptotic large-$`t`$ behavior sets in only for large values of $`t`$, i.e. for $`t\tau (k)`$, where the correlation function $`\stackrel{~}{G}(k,t)`$ is very small. Therefore, even if Eq. (71) breaks down for $`y3`$ and $`t`$ large, it still represents a very good approximation (even for $`y50`$) for the values of $`t`$ for which $`\stackrel{~}{G}(k,t)`$ is sizeable.

## Acknowledgements

V.M.-M. is a Ramon y Cajal research fellow and is partly supported by MCyT (Spain), project Nos. FPA2001-1813 and FPA2000-0956.

## A Integrals entering the field-theoretical calculations

In this appendix we report some integrals that enter the perturbative field-theoretical calculations.

In the two-loop calculation of the response scaling function in the high-temperature phase, cf. Eq. (45), one needs to compute the function

$$A(k^2,\omega )=\frac{2}{27}N_d^2\left[I(k,\omega )I(0,0)\right],$$
(A1)
with (dimensional regularization near four dimensions is understood)

$$I(k,\omega )=_0^{\mathrm{}}𝑑te^{i\omega t}\frac{d^dp_1}{(2\pi )^d}\frac{d^dp_2}{(2\pi )^d}\underset{i=1}{\overset{3}{}}\frac{1}{p_i^2+1}e^{t_i(p_i^2+1)},$$
(A2)
where $`p_3=kp_1p_2`$, and $`N_d=2/[(4\pi )^{d/2}\mathrm{\Gamma }(d/2)]`$. The integral $`I(k,\omega )`$ can be written in the form

$$N_d^2I(k,\omega )=\frac{1}{3}\frac{1}{(4\pi )^dN_d^2}_0^1t^2𝑑t_0^{\mathrm{}}𝑑se^ss^{3d}e^{i\omega s\frac{1t}{3}}_0^1u𝑑u𝑑v\frac{e^{s\frac{Q}{\mathrm{\Delta }}k^2}}{\mathrm{\Delta }^{d/2}},$$
(A3)
where

$`\mathrm{\Delta }=t^2u[1u+uv(1v)]+{\displaystyle \frac{1}{3}}(1t^2),`$ (A4)
$`Q=t^3u^2(1u)v(1v)+{\displaystyle \frac{1t}{3}}t^2u[1u+uv(1v)]+{\displaystyle \frac{(1t)^2}{9}}t+{\displaystyle \frac{(1t)^3}{27}}.`$ (A5)
We will also need the singularity structure of the $`A(k^2,\omega )`$. For this purpose, we will determine the large-$`t`$ behavior of $`\stackrel{~}{A}(k^2,t)`$ which is the Fourier transform with respect to $`\omega `$ of $`A(k^2,\omega )`$. This behavior can easily be derived from Eq. (A2). Setting $`p_1=k/3+q_1/\sqrt{t}`$, $`p_2=k/3+q_2/\sqrt{t}`$, $`p_3=k/3+q_3/\sqrt{t}`$, we obtain that for $`t\mathrm{}`$,

$`\stackrel{~}{A}(k^2,t)`$ $``$ $`{\displaystyle \frac{2}{27}}N_d^2t^de^{k^2t/33t}{\displaystyle \frac{1}{(k^2/9+1)^3}}{\displaystyle \frac{d^dq_1}{(2\pi )^d}\frac{d^dq_2}{(2\pi )^d}e^{(q_1^2+q_2^2+q_3^2)}[1+O(t^{1/2})]}`$ (A6)
$`=`$ $`{\displaystyle \frac{1}{(k^2/9+1)^3}}{\displaystyle \frac{(1/3)^{d/2}}{54}}\mathrm{\Gamma }^2(d/2)t^de^{k^2t/33t}[1+O(t^{1/2})]`$ (A7)
$``$ $`{\displaystyle \frac{1}{(k^2+9)^3}}{\displaystyle \frac{3e^{k^2t/33t}}{2t^4}}[1+O(t^{1/2})]+O(ϵ)`$ (A8)
This result implies the presence of a branching cut in $`A(k^2,\omega )`$ starting at $`\omega =i(3+k^2/3)`$.

The one-loop expression of the response function in the ordered phase, cf. Eq. (57), is written in terms of the function

$$B(k^2,\omega )=N_d^1\left[J(k,\omega )J(0,0)\right],$$
(A9)
where

$$J(k,\omega )=_0^{\mathrm{}}𝑑te^{i\omega t}\frac{d^dp_1}{(2\pi )^d}\underset{i=1}{\overset{2}{}}\frac{1}{p_i^2+1}e^{t_i^2(p_i^2+1)},$$
(A10)
with $`p_2=kp_1`$. The function $`B(k^2,\omega )`$ can be written in the form

$$B(k^2,\omega )=\frac{1}{2}_0^1t𝑑t𝑑u\frac{i2\omega (1t)[1t^2+4t^2u(1u)]k^2}{4+[1t^2+4t^2u(1u)]k^22i\omega (1t)}.$$
(A11)
Such an integral can be computed exactly obtaining

$`B(k^2,\omega )`$ $`=`$ $`{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{k^2}}\mathrm{ln}{\displaystyle \frac{4+2i\omega +k^2}{42i\omega +k^2}}+{\displaystyle \frac{2i}{\omega }}\sqrt{{\displaystyle \frac{4+k^2}{k^2}}}\mathrm{ln}{\displaystyle \frac{\sqrt{4+k^2}\sqrt{k^2}}{\sqrt{4+k^2}+\sqrt{k^2}}}`$ (A13)
$`{\displaystyle \frac{1}{k^2}}\mathrm{ln}{\displaystyle \frac{4\omega ^2+(4+k^2)^2}{16}}+{\displaystyle \frac{2}{\omega k^2}}F\left(\mathrm{ln}{\displaystyle \frac{F\omega }{F+\omega }}+\mathrm{ln}{\displaystyle \frac{F+\omega +ik^2}{F\omega ik^2}}\right)`$
with $`F\sqrt{\omega ^2+2i\omega k^2k^2(4+k^2)}`$. It is easy to see using this exact expression or repeating the argument presented for $`A(k^2,\omega )`$ that $`B(k^2,\omega )`$ is singular for $`\omega =i(2+k^2/2)`$.


text:
# Very High Redshift X-ray Selected AGN in the Chandra Deep Field-North

## 1. Introduction

Most of our current information on very high redshift active galactic nuclei (AGN) relates to the very high luminosity tail of the luminosity function (LF). At $`z>5`$ the Sloan Digital Sky Survey (SDSS) probes to an $`1450`$ Å absolute magnitude of about $`26`$ (Fan et al. 2001a, b), much brighter than the average luminosity of AGN at lower redshifts (Pei 1995). Calculations of the contributions of high-redshift AGN to the ionizing flux in the Universe must therefore rely on large extrapolations (e.g., Madau, Haardt, & Rees 1999), and deeper observations of the LF are needed to constrain this quantity directly. Such observations can also test models of supermassive black hole formation (e.g., Haiman & Loeb 1998, 1999; Haehnelt, Natarajan, & Rees 1998).

We can efficiently search for faint, high-redshift AGN using combined ultradeep X-ray and optical imaging. Since at high redshifts observed-frame X-ray bands correspond to very high rest-frame energies, we can obtain a relatively complete sample of AGN, including sources surrounded by very high column densities of gas and dust. Some Compton-thick sources may be omitted from the sample but would not be expected to contribute much to the ionizing light. As there are only a handful of stars in X-ray samples, there is little problem in separating red stars from high-redshift AGN. Here we use deep broad-band optical imaging and follow-up optical spectroscopy to measure the $`z>5`$ AGN population in the Chandra Deep Field-North (CDF-N) 2 Ms exposure. We assume $`\mathrm{\Omega }_M=1/3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=2/3`$ and $`H_o=65`$ km s-1 Mpc-1.

## 2. Data

The deepest X-ray image is the $`2`$ Ms CDF-N exposure centered on the Hubble Deep Field-North. D. M. Alexander et al. (in preparation) merged samples detected in seven X-ray bands to form a catalog of 503 point sources over an area of 460 arcmin2. Near the aim point the limiting fluxes are $`1.5\times 10^{17}`$ erg cm-2 s-1 ($`0.52`$ keV) and $`10^{16}`$ erg cm-2 s-1 ($`28`$ keV). The sensitivity degrades with off-axis angle beyond $`6^{}`$, rising to about $`10^{15}`$ erg cm-2 s-1 ($`0.52`$ keV) and $`10^{14}`$ erg cm-2 s-1 ($`28`$ keV) at the maximum angle of $`14^{}`$. The $`28`$ keV number counts and completeness levels in the 1 Ms exposure of the CDF-N are described in Cowie et al. (2002a).

Johnson $`B`$ and $`V`$, Cousins $`R`$ and $`I`$, and Sloan $`z^{}`$ images were obtained on the Subaru111The Subaru telescope is operated by the National Astronomical Observatory of Japan. 8.2 m telescope’s Prime Focus Camera (Suprime-Cam; Miyazaki et al. 2002). The observations and reductions are summarized in P. Capak et al. (in preparation). The exposure times are $`z^{}`$ (3.9 hrs), $`I`$ (2.9 hrs), $`R`$ (5.3 hrs), $`V`$ (8.4 hrs), and $`B`$ (1.7 hrs). The image quality is $`0.8^{\prime \prime }`$, except in the $`V`$ ($`1.0^{\prime \prime }`$) and $`R`$ ($`1.1^{\prime \prime }`$) bands. Magnitudes and images are given in Barger et al. (2002; 1 Ms sample) and A. J. Barger et al. (in preparation; 2 Ms). We chose the optical position to be the highest surface brightness peak in the $`z^{}`$ image within $`1.5^{\prime \prime }`$ of the X-ray position. For most sources the optical counterpart is unambiguous, but in some cases (about $`7\%`$), there is more than one possible counterpart (see, e.g., Fig. 2 of Barger et al. 2002). We discuss alternative counterparts in § 3.2. For off-axis sources in extended galaxies, and optically faint sources, we used the X-ray position.

We measured magnitudes in $`2^{\prime \prime }`$ diameter apertures on images smoothed to the lowest resolution image ($`R`$) to provide the most accurate colors. We estimated total magnitudes by correcting the $`2^{\prime \prime }`$ diameter magnitudes to $`6^{\prime \prime }`$ diameter using a median offset ($`0.8`$ for each filter) measured on the smoothed images. Three of the sources lie too close to bright stars for accurate photometry; we excluded these objects. We use the AB magnitude system with $`m_{AB}=2.5\mathrm{log}f_\nu 48.60`$, where $`f_\nu `$ is the flux of the source in units of erg cm-2 s-1 Hz-1. The $`1\sigma `$ limits on the total magnitudes are 28.5 ($`B`$), 28.4 ($`V`$), 28.2 ($`R`$), 27.3 ($`I`$), and 26.9 ($`z^{}`$) based on measuring randomly positioned apertures in the field. These are slightly more conservative (brighter) than those used in Barger et al. (2002), where the limits were computed for $`3^{\prime \prime }`$ diameter aperture magnitudes using the distribution of sky backgrounds in the individual pixels. The present limits deal better with correlated noise.

Figure 2 shows redshift versus $`0.52`$ keV flux. Barger et al. (2002) and A. J. Barger et al. (in preparation) give optical spectroscopic identifications for 249 of the X-ray sources, including 16 off-axis sources where the X-ray position is within the envelope of a bright galaxy but not within $`1.5^{\prime \prime }`$ of the nucleus (Hornschemeier et al. 2002). We analyze the colors for the 423 (out of the 500 with accurate photometry) X-ray sources with optical counterparts brighter than the $`5\sigma `$ limit of $`z^{}=25.2`$. In Fig. 2 we split the spectroscopically unidentified counterparts into $`z^{}<25.2`$ (open squares at $`z=0.6`$) and $`z^{}>25.2`$ (plus signs at $`z=0.3`$).

The solid curves in Fig. 2 are for fixed rest-frame $`28`$ keV luminosities ranging from $`10^{42}`$ erg s-1 (lowest curve) to $`10^{45}`$ erg s-1 (highest curve). The open diamonds show the $`0.52`$ keV fluxes of the highest redshift sources from the SDSS sample (Brandt et al. 2002; Mathur et al. 2002). These form an optically complete sample for $`z>5.7`$ and $`z^{}<20`$ over an area of 1500 deg2 (Fan et al. 2001b). The present X-ray sample probes almost two orders of magnitude fainter than the typical X-ray flux seen in the SDSS sample; however, the area of sky covered is four orders of magnitude smaller. If the rest-frame optical to rest-frame X-ray flux ratios of high-redshift, low-luminosity sources are similar to those of the high-redshift, high-luminosity sources in the SDSS sample, then the magnitudes of the optical counterparts to the low-luminosity sources would go as faint as $`z^{}=25`$, which is well-matched to our $`z^{}=25.2`$ color-selection limit. However, if instead high-redshift, low-luminosity sources have lower rest-frame optical to X-ray flux ratios (Green et al. 1995; Vignali, Brandt, & Schneider 2003), then the counterparts could be fainter than $`z^{}=25`$.

## 3. High Redshift AGN

### 3.1. Optically Bright Sample

High-redshift sources in the X-ray sample are easily identified because the increasing line blanketing of the Lyman forest at high redshifts produces extremely large breaks across redshifted Ly$`\alpha `$ (e.g., Fan et al. 2001b). Songaila & Cowie (2002) measured the average forest transmission in the redshift range $`z=3.86.2`$. The magnitude break across Ly$`\alpha `$ is

$$\mathrm{\Delta }m=3.8+20.3\mathrm{log}_{10}\left(\frac{1+z}{7}\right).$$
(1)
At $`z=5`$, where Ly$`\alpha `$ lies between $`R`$ and $`I`$, the break is 2.4 magnitudes, while at $`z=6`$, where Ly$`\alpha `$ lies between $`I`$ and $`z^{}`$, it is 3.8 magnitudes. If the epoch of reionization is at $`z6.1`$ (Becker et al. 2001), then the break is extremely large at redshifts higher than this. The Ly$`\alpha `$ line moves through our longest wavelength ($`z^{}`$) band at $`z6.5`$, which is the largest redshift where a source would be seen in our optical data.

Figure 3.1 shows $`Iz^{}`$ versus redshift (spectroscopically unidentified sources are at $`z=0.3`$) for our $`z^{}<25.2`$ sample. We see from Fig. 3.1 that none of the sources lie at $`z>6`$, since such sources would be redder than $`Iz^{}>1.5`$. In Fig. 3.1 we show $`VI`$ versus $`Iz^{}`$. The dashed line is the track with redshift of the composite SDSS quasar spectrum of Vanden Berk et al. (2001), extrapolated with a $`f_\nu \nu ^{0.46}`$ power-law below Ly$`\alpha `$, modulated by the Songaila & Cowie (2002) forest transmissions, and convolved through the Suprime-Cam filters. The extrapolation below the Lyman continuum limit is not critical since the high incidence of Lyman limit systems at these redshifts truncates the spectrum at shorter wavelengths. The dot-dash line shows the same calculation for a $`f_\nu \nu ^0`$ galaxy with no emission lines and a cut-off at $`912`$ Å. The deep forest absorption moves the tracks far from the galaxy populations. Apart from the one $`z=5.19`$ AGN (large solid square) with $`z^{}=23.9`$, there are no other $`z>5`$ candidates. Given the possibility that in a small number of cases we may have chosen the wrong optical counterpart or the colors are contaminated by a nearby object, we also visually searched differenced images for objects with red $`VI`$ colors. No further candidates were found.

### 3.2. Optically Faint Sample

Of the 500 X-ray sources in our sample, 77 have $`z^{}>25.2`$. Two of these have spectroscopic redshifts ($`z=3.40`$ and $`z=4.14`$) from Barger et al. (2002). Of the remainder, 44 are detected above the $`2\sigma `$ level in $`V`$ or $`B`$ and cannot lie at $`z>5`$, leaving 31 faint sources that could be at $`z>5`$. Many of these faint sources may lie at lower redshifts (e.g., Yan et al. 2002), and some may be false X-ray detections; about half lie within a factor of 2 of the flux limit at their off-axis radius, so this is an upper bound. Within a $`6^{}`$ off-axis radius where the sensitivity is relatively uniform, 9 of these optically faint sources could have rest-frame $`28`$ keV luminosities in excess of $`10^{43}`$ erg s-1, if they were placed in the $`z=56.5`$ redshift range. Only 7 of these are not detected in $`R`$ or $`I`$ and could lie beyond $`z=6`$. The cumulative surface density of field sources to $`z^{}=25.2`$ is $`1.2\times 10^5`$ deg-2, so we expect that about 30 of our sources with bright optical identifications may be chance projections, with the true optical counterparts being fainter. If we assume that these 30 sources are similar to the 77 isolated faint sources discussed above, then we can include these sources in our candidate $`z>5`$ population with a scaling factor of 1.4.

Haiman & Loeb (1999) predicted a surface density of $`50`$ $`z>5`$ sources in a $`6^{}`$ radius circle above a flux of $`2\times 10^{16}`$ erg cm-2 s-1 (observed $`0.46`$ keV). We find less than four such $`z>5`$ objects, or a maximum of 5.6 if we use the 1.4 scale factor discussed above. The surface density is at least an order of magnitude lower than the Haiman & Loeb estimate, strengthening the similar conclusion reached by Alexander et al. (2001) and Hasinger (2002) by about a factor of three.

## 4. Discussion

We compute the number density of $`z=56.5`$ sources in the rest-frame $`28`$ keV luminosity range $`10^{43}`$ to $`10^{44}`$ erg s-1 following the procedures described in Cowie et al. (2002b). In Fig. 4a we show the number density obtained using only the one $`z>5`$ spectroscopic identification (solid diamond). We also computed an upper limit to this number density by assigning all the optically faint sources that could lie in the redshift interval redshifts at the center of the interval; the horizontal bar corresponds to those with $`L_x`$ in the specified luminosity range. In Fig. 4a we compare our $`z=56.5`$ number density with the lower redshift number densities from Cowie et al. (2002b) for both the $`10^{43}`$ to $`10^{44}`$ erg s-1 (solid circles) and $`10^{44}`$ to $`10^{45}`$ erg s-1 ranges (open circles). Again, the symbols denote the number densities obtained from only the spectroscopically identified AGN, while the horizontal bars show the upper limits. The bars are not consistent with one another because all the unidentified sources are included in each redshift bin if they lie in the specified luminosity range. To compare with the optically-selected samples of Boyle et al. (2000) and Fan et al. (2001b), in Fig. 4b we show the results for an $`\mathrm{\Omega }_M=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$ cosmology with $`H_o=50`$ km s-1 Mpc-1. The small number uncertainties are large in Fig. 4, but the number densities in the $`10^{43}`$ to $`10^{44}`$ erg s-1 range (AGN whose luminosities would classify them as Seyferts) show a slow decrease with increasing redshift from $`z=0`$ to $`z=6`$. This result holds even within the systematic uncertainties shown by the the horizontal bars. In contrast, the number densities in the $`10^{44}`$ to $`10^{45}`$ erg s-1 range (AGN with quasar luminosities) peak at $`z=1.53`$ and closely match in shape the evolution of the optically-selected samples. Figures 4a and b show that this conclusion does not depend on the cosmological geometry.

The number density at $`z=56.5`$ for $`L_x10^{43}`$ to $`10^{44}`$ erg s-1 is about three orders of magnitude higher than the optically-selected SDSS sample at $`L_x10^{45}`$ erg s-1. While this conclusion is based on the single $`z>5`$ object, the fact that other high-redshift AGN have also been found in the relatively small Chandra fields (Silverman et al. 2002; $`z=4.93`$) and in the small (74 arcmin2) optical field of Stern et al. (2000; $`z=5.5`$) gives us confidence that the CDF-N field is not anomalous. If we use our spectroscopically identified point, the slope $`\beta `$ of the LF ($`\varphi (L)dL=L^\beta dL`$) would be about $`2.6`$, which is similar to the bright end slope at lower redshifts in both the optical (Pei 1995) and X-ray (Miyaji et al. 2000). If we include the unidentified objects, the maximum value of $`\beta `$ would be about 3. A detailed LF determination for this redshift interval will require much larger area X-ray samples.

Fan et al. (2001b) investigated the AGN ionizing flux by fitting a variety of power-laws to extrapolate the bright SDSS data. A power-law with our observed slope of $`\beta =2.6`$ matched to the SDSS data fails by more than an order of magnitude to ionize the intergalactic medium at these high redshifts, even extrapolating to faint magnitudes (their Fig. 10). We may derive this result more directly by calculating the number of ionizing photons ($`n_i`$) per baryon ($`n_b`$) produced in the redshift interval by the observed AGN. We first sum the $`K`$-corrected fluxes of the sources to determine the rest-frame $`1450`$ Å background light they produce and then convert this to the number density of ionizing photons using the form of the near-ultraviolet AGN spectrum given in Madau, Haardt, & Rees (1999), which gives a ratio $`\eta =4.7\times 10^{25}`$ ionizing photons per erg Hz-1 at $`1450`$ Å. ($`\eta `$ would be less for a galaxy shaped spectrum.) We therefore have

$$\frac{n_i}{n_b}=\frac{4\pi \eta }{cAn_b}\underset{i}{}f_i.$$
(2)
The $`f_i`$ are the fluxes corresponding to the $`z^{}`$ magnitudes for $`z5.7`$, and $`A`$ is the observed solid angle in steradians; we set $`n_b=2.0\times 10^7`$ cm-3. In the central $`6^{}`$ of the CDF-N field, which has the deepest uniform sensitivity and where the $`z=5.19`$ source is found, we obtain $`n_i/n_b0.07`$. There is almost no change if we include all the possible $`z^{}>25.2`$ sources in the region because these objects are extremely faint in $`z^{}`$ and hence contribute little to the sum in Eq. 2. For a Poisson distribution there is less than a 2% probability of seeing 0 or 1 sources when the true value is greater than 6, so at $`98\%`$ confidence we can conclude that $`n_i/n_b<0.4`$. Thus, there are too few moderate luminosity AGN at $`z=56.5`$ to ionize the intergalactic medium.

Support came from the University of Wisconsin Research Committee (A.J.B.), the Alfred P. Sloan Foundation (A.J.B), NSF grants AST-0084847 (A.J.B.), AST-0084816 (L.L.C.), and AST-9983783 (D.M.A., W.N.B.), NASA grants DF1-2001X (L.L.C.) and NAS 8-01128 (G.P.G.), and CXC grant G02-3187A (D.M.A., F.E.B., W.N.B.).