format: text

text:
# A first principles simulation of rigid water

## I Introduction

The importance of water in many areas of science has motivated a large number of experimental and theoretical investigations of the liquid. However, it is only recently that x-ray and neutron diffraction measurements have come to an overall agreement for properties such as the structure of water at ambient conditions jsorenson00 ; asoper00b . In addition, several details of the water structure remain the subject of debate, and many dynamical properties of the liquid are not yet well understood.

Theoretical models have played an important role in the interpretation of experimental measurements and in understanding the physical properties of water bguillot02 . Over the last thirty years, rather accurate empirical force fields have been developed, which can reproduce not only the structure but also many dynamical properties of the liquid fstillinger74 ; hberendsen81 ; mmahoney00 . Although empirical models work well for pure water under ambient conditions, they are usually difficult to generalize to complex solutions or thermodynamic states far from ambient conditions. For example, the majority of empirical water models that are in current use employ potentials that do not change depending on the environment.

In recent years, it has become possible to simulate the properties of a liquid entirely from first principles, without having to resort to fitted potentials. This is due in large part to the development of the Car-Parrinello (CP) method rcar85  along with the continual increase in high performance computing resources. Although rather accurate and with the potential of being a truly predictive tool, CP simulations are much more computationally intensive than classical simulations.

The case of water is particularly demanding for CP simulations. The ionic vibrational spectrum of the liquid exhibits high frequency modes, i.e. O-H stretch (3200 to 3600 cm-1) and H-O-H bending modes ($``$1600 cm-1). Therefore, in order to avoid a coupling between the ionic and electronic degrees of freedom, which could cause severe inaccuracies in a CP simulation, a relatively small fictitious electronic mass ( $`\mu `$ 400 a.u. for protonated water) needs to be used jgrossman03 . In turn, the use of small values of $`\mu `$, together with the high kinetic energy cutoff required to describe the oxygen pseudopotential in a plane wave description, necessitates the use of small integration time steps. The time step may need to be as small as 0.08 fs, which is approximately ten times smaller than what is often used in classical MD simulations with empirical inter-atomic potentials. This poses a severe restriction on the time scales that can be accessed in CP simulations of water.

We note that when using Born-Oppenheimer (BO) dynamics (where the total energy of the system is minimized at each ionic step), one can use larger time steps than in CP simulations, since electronic degrees of freedom are not propagated at the same time as ionic coordinates. However, the accuracy required to reduce systematic errors on the ionic forces so as to have conservative dynamics is such that large number of iterations are usually necessary to minimize the Kohn-Sham energy at each ionic step. Therefore, the gain in efficiency obtained with a large time step is more than counter-balanced by the computational time requirement for total energy minimizations.

It is interesting to note that the problem of integrating fast vibrational modes in simulations of liquid water has also been encountered in classical MD simulations, where the most common approach has been to completely eliminate the high-frequency intra-molecular motion by using bond length and angle constraints wgunsteren82 ; tschlick97 . Based on this approach, a variety of classical water potentials, e.g. the TIP series mmahoney00 , are capable of accurately reproducing many of the interesting properties of water. In particular, results obtained with the rigid water TIP5P potential are in very good agreement with a variety of experimental measurements such as the structure, the temperature of maximum density, diffusion, as well as dielectric properties mmahoney00 ; mmahoney01 ; mmahoney01b .

With the aim of investigating how to increase the integration time step in CP simulations of water and thus access larger time scales, we have carried out calculations using a rigid water approximation. In this paper, we present the results of these simulations and compare them to those obtained with flexible water molecules (i.e. without imposing any constraints on the geometry of the molecules in the liquid), and we discuss the effect of a rigid model on the structural properties of water at ambient and high pressure conditions. Our results show that an ab initio rigid water model yields faster diffusion and radial distribution functions which are less structured than those found with a flexible model. Overall, the properties computed with the rigid model are in better agreement with experiment than those determined with a flexible model. Possible reasons for this improved agreement are discussed. In addition, we present a localized orbital analysis of the trajectories obtained with both a rigid and flexible water model, and we demonstrate that the large dipole moment changes in going from the gas to the liquid phase are not significantly altered by the rigid water approximation. The use of a rigid water model in ab initio simulations opens the way to much longer simulations of solutions where chemical reactions and dissociation effects do not occur.

## II Methods

In order to examine how the structural and dynamical properties of water are altered by a rigid water approximation, we have performed a series of first principle molecular dynamics simulations gp  of water with and without intramolecular bond and angle constraints under ambient and high pressure and temperature conditions. The simulations consist of 54 water molecules in a periodically repeated cubic cell with a lateral dimension of either 11.74 Å or 10.10 Å, which correspond to densities of 1.00 g/cc and 1.57 g/cc, respectively. At each density, we have compared simulations where the intra-molecular geometry of the water molecules are rigid to those where the geometries of the water molecules are fully flexible.

In each simulation, the electronic structure was described within density functional theory (DFT) phohenberg64 ; wkohn65  with the PBE generalized gradient approximation jperdew98 . The valence wavefunctions and charge density were expanded in a plane wave basis, which was truncated in reciprocal space at 85 and 340 Ry, respectively. Norm-conserving pseudopotentials of the Hamman type were used to describe valence-core interactions dhamann89 ; lkleinman82 .

The simulations were performed with the CP technique rcar85 , which is based on the use of a Lagrangian that couples together the system’s electronic and ionic degrees of freedom,

$`L_{CP}`$ $`=`$ $`\mu {\displaystyle \underset{i}{}}f_i{\displaystyle 𝑑𝐫\left|\dot{\mathrm{\Psi }}_i\left(𝐫\right)\right|^2}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{I}{}}M_I\dot{𝐑}_I^2`$ (1)
$`E_{KS}[\left\{\mathrm{\Psi }_i\right\},𝐑_I]+`$
$`{\displaystyle \underset{ij}{}}\mathrm{\Lambda }_{ij}\left({\displaystyle 𝑑𝐫\mathrm{\Psi }_i^{}\left(𝐫\right)\mathrm{\Psi }_j\left(𝐫\right)}\delta _{ij}\right).`$
In Eq. 1, $`\mathrm{\Psi }_i\left(r\right)`$ are the Kohn-Sham orbitals, $`\mu `$ is a fictitious mass parameter used to evolve the electronic degrees of freedom in time, $`M_I`$ are ion masses, $`E_{KS}`$ is the Kohn-Sham energy, $`f_i`$ are occupation numbers, and $`\mathrm{\Lambda }_{ij}`$ are Lagrange multipliers, used to impose the orthonormality constraint $`\mathrm{\Psi }_i^{}\mathrm{\Psi }_j=\delta _{ij}`$. The equations of motion derived from Eq. 1 are,

$$\mu \ddot{\mathrm{\Psi }}_i(𝐫,t)=\frac{\delta E}{\delta \mathrm{\Psi }_i^{}(𝐫,t)}+\underset{j}{}\mathrm{\Lambda }_{ij}\mathrm{\Psi }_j(𝐫,t)$$
(2)
and

$$M_I\ddot{𝐑}_I=\frac{E}{𝐑_I\left(t\right)}.$$
(3)
Central to the CP method is the introduction of the fictitious mass parameter, $`\mu `$, which enables the electronic and ionic trajectories to be propagated simultaneously at each time step. The particular value of $`\mu `$ is chosen so that the dynamics of the electronic degrees of freedom occurs on a time scale that is much faster than, and thus decoupled from, the ion dynamics.

For a given system, an appropriate value of $`\mu `$ can be determined by approximating the dynamics of the orbitals generated by Eq. 2 as a superposition of oscillators with frequencies, gpastore91 ; ggalli91

$$\omega _{ij}=\left[\frac{2\left(ϵ_jϵ_i\right)}{\mu }\right]^{1/2},$$
(4)
where $`ϵ_i`$ and $`ϵ_j`$ are Kohn-Sham eigenvalues of occupied and unoccupied states, respectively. The lowest frequency obtained from Eq. 4 occurs when $`ϵ_i`$ is the highest occupied state (HOMO) and $`ϵ_j`$ is the lowest unoccupied state (LUMO). As discussed in Ref. gpastore91, , although Eq. 4 is a rather crude approximation of the true orbital dynamics, it still provides a useful estimate for selecting an appropriate value of $`\mu `$ jgrossman03 .

In the case of water, the highest ionic frequency is due to O-H stretching modes, which are at $``$3500 cm-1 within density functional theory. In addition, the HOMO-LUMO gap for water within density functional theory is approximately 4.6 eV klaasonen93 . As discussed in Ref. jgrossman03, , this means that values of $`\mu `$340 au are needed to ensure a clear separation between the ionic and electronic degrees of freedom in CP simulations of flexible, protonated water. However, when the rigid water approximation is used, much larger values of $`\mu `$ can be used because the bond distance and angle constraints suppress the high frequency ionic motions. The various values of $`\mu `$ that we have used for both the rigid and the flexible water simulations are listed in Table 1.

The initial starting configurations of the simulations were generated from previous simulations of flexible water eschwegler01  by scaling the intramolecular O-H distances and H-O-H angles, while leaving the spatial orientations of the molecules unchanged. The particular values for the geometry, $`0.9926`$ Å for the O-H distance and $`104.6^{}`$ for the H-O-H angle, were chosen to be the same as the peak values for the intra-molecular O-H distance and H-O-H angular distributions obtained in the previous flexible water simulations eschwegler01 . These bond distances and angles are similar to the ones used in the well-known TIP classical water models ($`0.9572`$ Å, $`104.52^{}`$) mmahoney00 . In order to maintain the distance and angle constraints at every iteration, the SHAKE algorithm was used jryckaert77 ; mpallen87 .

In simulations A and B (Table 1), the samples were initially heated to a temperature of 450 K and cooled to 300 K for a period of 2 ps. In simulation D, the system was heated to a temperature of 900 K and cooled to 600 K over 4 ps. The thermostats were then removed and the simulations were run under constant energy conditions for the times listed in Table 1.

Maximally localized Wannier functions (MLWF) nmarzari97 ; psilvestrelli98 , analogous to the orbitals obtained by the Boys localization procedure sboys60 , were used to determine how the rigid water approximation may influence the electronic structure of the water molecules. By using a recently proposed joint approximate diagonalization scheme, we were able to compute the MLWFs “on-the-fly” for a large number of configurations fgygi03 . Following the procedure proposed by Silvestreli et al. psilvestrelli99 , the centers of the MLFWs have been used to compute an approximate dipole moment of each of the water molecules in the simulation.

## III Results

The primary effect of the rigid water approximation is best illustrated by comparing the vibrational spectra obtained from a rigid and a flexible water simulation. In Fig. 1, the power spectra collected from simulations A and C are shown. The spectra were computed directly from the velocity autocorrelation function with the maximum entropy method. The two highest frequency modes in the flexible water simulation C are completely absent in simulation A. These modes at 1600 cm-1 and 3150 cm-1 correspond to intra-molecular H-O-H bending and O-H stretching, respectively. By suppressing the high frequency modes large values of $`\mu `$ can be used, which in a fully flexible simulation would normally lead to a coupling of the electronic and ionic degrees of freedom, and in turn to a breakdown in the adiabiticity of the CP dynamics jgrossman03 .

In addition to the obvious absence of high frequency modes in the rigid water simulation, there are small differences in the lower frequency range of the spectra below $``$1000 cm-1. For example, the librational modes are slightly shifted to lower frequencies by approximately 60 wavenumbers in the rigid water simulation A. Overall, the differences in the low frequency region of the rigid and flexible water simulations are small.

In Fig. 2, the oxygen-oxygen radial distribution functions (RDFs) obtained from simulations A, B and C are shown. Also displayed in Fig. 2, is the latest experimental oxygen-oxygen RDF obtained by an analysis of neutron diffraction data asoper00b , which is nearly identical to the RDF determined from a recent x-ray diffraction study of water jsorenson00 . The small differences between simulations A and B indicate that the choice of $`dt`$=0.07 or 0.24 fs has little effect on the oxygen-oxygen RDF when the rigid water approximation is used. The differences in the height of the first peak and the first minimum in the RDFs are most likely due to the $``$10 K higher average temperature of simulation A over simulation B (see Table 1). However, as discussed in Ref. jgrossman03, , for the simulation times used here (15-25 ps), differences of this magnitude are within the expected error bars. In Fig. 2 we also show the oxygen-oxygen RDF obtained from simulation C with flexible water molecules. The differences between the rigid and the flexible water simulations are significant. The first peak height is approximately 0.4 higher in simulation C than in simulations A and B, and is shifted inward by 0.03 Å. In general, the RDF peak heights and depths are decreased in the rigid water simulations, showing that the rigid water approximation causes an overall decrease in the structure of the liquid. The same trend in terms of the peak heights and depths has been observed in flexible and rigid water simulations by using the SPC classical interaction potential jlobaugh97 .

The tendency of the rigid water approximation to decrease the structure of water can also be seen in the oxygen-hydrogen and the hydrogen-hydrogen RDFs, as shown in Figs. 3 and 4. It is interesting to note that in simulation C the second peak in the oxygen-hydrogen RDF at r$``$2.8 Å is larger in height than the third peak at r$``$3.2 Å, whereas in simulations A and B, as well as in the experimental measurement asoper00b , the second peak is smaller than the corresponding third peak. Since the second oxygen-hydrogen peak coincides with pairs of hydrogen bonded water molecules in the liquid, the relative heights of the second and third peaks have been used as an indicator of the amount of hydrogen bonding in the liquid jlobaugh97 .

Overall, Figs. 2 to 4 suggest that the main effect of the rigid water approximation is to remove a large fraction of the over-structure that is characteristic of flexible water as described within DFT both by CP jgrossman03  and BO molecular dynamics dasthagiri03 , and in general, the rigid water distribution functions appear to be in better agreement with experiment. There are several possible reasons for this improved agreement. For example, it is possible that the rigid water approximation mimics some of the effects coming from the quantum nature of protons, which are explicitly ignored in simulations that allow for intramolecular flexibility with classical dynamics itironi96 ; jlobaugh97 ; bhess02 . More specifically, in Ref. itironi96,  it was pointed out that at a temperature of 300 K, $`k_BT200`$ cm-1, whereas the high frequency intramolecular modes in water range from $``$1000 to 3500 cm-1. In the quantum system the amount of thermal energy available to excite vibrational modes is much smaller than the lowest possible intramolecular vibrational excitations,

$$\mathrm{}\omega k_BT.$$
(5)
This indicates that the real quantum system will essentially be restricted to its vibrational ground state at 300 K. Therefore, the quantum system may be better described by a classical rigid water model than by a classical flexible model, despite the fact that in the quantum system the protons become delocalized compared to the classical case.

The effect of the proton quantum motion can be considered with path-integral (PI) methods rfeynman65 . To date, all PI simulations with empirical interaction potentials have found an overall softening in the RDFs of water at ambient conditions rkuharski85 ; gbuono91 ; jlobaugh97 ; hstern01 ; mmahoney01 . As pointed out in Ref. gbuono91,  these structural changes are quite similar to a 50 K increase in the simulation temperature. However, it is not clear if the lack of quantum effects can fully account for all of the overstructure that appears to be characteristic in the DFT/GGA treatment of water. In addition, a recent PI-DFT simulation bchen03  has found the surprising result that quantum effects may enhance hydrogen bonding in water.

In order to examine how the dynamical properties of water are affected by the rigid water approximation, we have estimated the diffusion coefficients, D, for the simulations of rigid and flexible water. To compute the self diffusion coefficient, the mean square displacement (MSD) of the oxygen atoms in the simulations were tracked as a function of the simulation time. The statistical sampling of the simulation data was improved by defining multiple starting configurations separated by 4 fs. The slope of the resulting MSD in the range of 1 to 10 ps was then determined and used to estimate D according to the Einstein relation mpallen87 ,

$$6D_1=\underset{t\mathrm{}}{lim}\frac{d}{dt}\left|𝐫_i\left(t\right)𝐫_i\left(0\right)\right|^2.$$
(6)
The diffusion coefficient can also be estimated via integration of the oxygen velocity autocorrelation function as

$$D_2=\frac{1}{3}_0^{\mathrm{}}𝑑t𝐯\left(t\right)𝐯\left(0\right),$$
(7)
which for the simulation times used here yields nearly identical results as Eq. 6. The measured diffusion coefficients for the different simulations as computed by Eqs. 6 and 7 are reported in Table 2. For the rigid water simulations A and B, the diffusion coefficients are approximately 6.2 and 3.2 times larger, respectively, than the flexible water simulation C. Although simulations A and B were both performed at higher temperatures than simulation C, the observed increases in the diffusion coefficients are outside the increase (1.4 to 1.7 fold) expected for a 15 to 26 increase in temperature over ambient, based on experimental data kkrynicki78 . Given that the rigid water approximation results in a less structured liquid than the flexible one, it is not surprising that it also leads to faster diffusion, which is in closer agreement with the experimental measurement of $`2.4\times 10^5cm^2/s`$ rmills73 ; kkrynicki78 .

In addition to possible changes in structure and dynamics, we have examined how the rigid molecule approximation changes the electronic properties of water. In order to do this we have performed a localized orbital analysis by computing the MLWFs for a series of well-separated snapshots from the rigid water simulation B and the flexible water simulation C. Within the pseudopotential approximation, there are four doubly occupied MLWFs around each of the water molecules in the simulations. Two of the MLWFs are localized on the oxygen-hydrogen covalent bonds, and the other two are localized on the lone-pair locations of the oxygen atoms. Given the large amount of data that the MLWFs represent, in the following, we only consider the centers of the MLWFs rather than the orbitals themselves. In Fig. 5, the oxygen-MLWF center RDFs for simulations B and C are shown. The RDFs consist of two distinct distributions centered at r$``$0.33 Å and r$``$0.49 Å, which correspond to lone-pair and covalent bond locations, respectively. For comparison, the dashed vertical lines in Fig. 5 represent the locations of the MLWF centers around the oxygen atom of an isolated gas phase water molecule. Surprisingly, the rigid water approximation does not significantly alter the large changes in the MLWF centers that are expected when going from an isolated water molecule to the liquid state. As can be seen in Fig. 5, the lone pair distributions are shifted away from the oxygen atoms by $``$0.03 Å, and the covalent bond distributions are shifted toward the oxygen atoms by $``$0.04 Å.

As proposed by Silvestrelli et al. psilvestrelli99 , an approximate dipole moment for each water molecule in the liquid can be defined by assigning the total charge of each MLWF to a point charge located at its corresponding center. Because the MLWF on neighboring water molecules in the liquid do not significantly overlap, this provides a less ambiguous definition of the molecular dipole moments than arbitrarily assigning electron density to individual water molecules. As pointed out in Ref. apasquarello03, , dipole moments computed in this manner from static configurations may not be representative of the experimentally measured dipole moments in the fluid. However, the MLWF dipole moments are useful for examining qualitative differences in the polarization of water as a function of different approximations or of solutes present in the liquid.

In Fig. 6, the probability distributions of molecular dipole moments for rigid and flexible water systems calculated from the MLWF centers are shown. Fig. 6 indicates that the rigid water approximation has a rather small effect on distribution of dipole moments in the liquid. In particular, the average moment in simulation B is shifted to 3.08 Debye as compared to 3.20 Debye in simulation C. Apparently, an explicit description of high frequency O-H stretch and H-O-H bending modes is not necessary to reproduce the broad range of moments that are characteristic of the liquid. It is also interesting to note that the latest experimental estimate based on an analysis of the x-ray structure factor of water indicates that the dipole moment of water in the liquid is 2.9 Debye ybadyal00 , which is in closer agreement with the rigid water model than the flexible simulation.

The decrease in the average dipole moment obtained from the rigid water approximation offers another explanation for the observed softening of the liquid structure. In addition to mimicking quantum effects, it is possible that the rigid water approximation to some extent corrects for the general tendency of simple GGA-based functionals to overestimate the polarizability of molecules smcdowell95 ; pcalaminici98 ; acohen99 ; ccaillie00 . For example, the static isotropic polarizability of an isolated water molecule is 10.74 au with the PBE functional as compared to the experimental value of 9.64 au ccaillie00 . It is interesting to note that hybrid DFT functionals, which include some amount of Hartree-Fock exchange, appear to significantly improve on the polarizability of water. In particular, the average polarizability of the water molecule is 9.78 au with hybrid PBE0 functional ccaillie00 .

In addition to water under ambient conditions, we have also examined how the rigid water approximation affects the properties of water under extreme temperatures and pressures. In particular, we have performed a simulation of rigid water (simulation D) and of flexible water (simulation E) at a density of 1.57 g/cc and an average temperature of $``$600 K. These high density and temperature conditions correspond to a regime where molecular dissociation is still considered a rare event eschwegler00 ; eschwegler01 . However, the pressure ($``$10 GPa) is high enough to cause a large increase in the nearest neighbor coordination of each water from 4.5 at ambient conditions to nearly 13 at high pressure eschwegler00 . The oxygen-oxygen, oxygen-hydrogen and hydrogen-hydrogen RDFs for simulations D and E are compared in Figs. 7 to 9. Except for the expected intramolecular differences due to the constraints, the RDFs obtained from simulations D and E are remarkably similar. In particular, both the large increase in the number of nearest neighbors as well as the stiffness of the first peak in gOO(r) as a function of compression are reproduced by the rigid water model eschwegler00 . As higher densities and temperatures are considered, intramolecular dissociation will become an common event in flexible water simulations eschwegler01  and the rigid water approximation is expected to be inappropriate for the description of the liquid.

## IV Conclusion

In summary, we have used a series of Car-Parrinello molecular dynamics simulations to examine how the rigid water approximation affects the computed properties of water in the liquid state at ambient conditions. In agreement with previous observations based on empirical interaction potentials jlobaugh97 , the rigid water approximation is found to cause an overall decrease in structure and an increase in diffusion of the liquid. These changes result in properties that are in better agreement with experimental measurements than the corresponding first principles simulations with flexible water molecules. At higher temperatures and densities in a regime where intramolecular dissociation is still a rare event, the differences between simulations where water molecules are either rigid or flexible become negligible.

In addition to an improved structural and dynamical description of water, the rigid water model enables the use of time steps as large as 0.24 fs within the Car-Parrinello scheme (i.e. $``$3 times larger than in a flexible water simulation). A similar conclusion was reached in Ref. etsuchida01,  for first-principle simulations of a cytosine molecule in the gas phase. This represents an important advantage for first-principle simulations of aqueous solutions where chemical reactions do not occur, and opens up the possibility of investigating phenomena that take place on a long timescale. For example, understanding how water orients around a hydrophobic solute may require simulations of the order of 100 to 200 ps. The rigid water approximation presented here may prove to be an accurate and efficient approach for describing the interaction between a hydrophobic solute and water within a first-principles context.

The authors would like to thank J. C. Grossman for many useful discussions. This work was performed under the auspices of the U.S. Dept. of Energy at the University of California/Lawrence Livermore National Laboratory under contract no. W-7405-Eng-48.


text:
# Reply to Comment on “Roughness of Interfacial Crack Fronts: Stress-Weighted Percolaton in the Damage Zone”

Alava and Zapperi (AZ) az03  question whether the fracture fronts we observe in shb03  are self affine. We use Family-Vicsek scaling to determine the two scaling exponents $`\alpha `$ and $`z`$. AZ claim that this is not enough to determine whether the front is self affine and they go on to point to the presence of overhangs as their evidence of fractality rather than self affinity. In Fig. 1, we show details of experimental fracture fronts from Fig. 1 of ms01 . There are significant overhangs, but these experimental fronts are self affine. In fact, as long as one allows a damage cloud to develop, overhangs are unavoidable. Family-Vicsek scaling implies non-isotropic scaling. A consequence is that the front width $`w`$ scales with the width of the system $`L_x`$. Non-isotropic scaling is the defining property of self-affine surfaces.

AZ claim that the source of the $`L_x`$-dependency of the damage length scale $`l_y`$, Eq. (9) in shb03 , comes from the rescaling of the model’s elastic constants when changing $`L_x`$. However, this rescaling is necessary to ensure that the elastic properties remain unchanged when the system size is changed under uniform loading conditions. It is caused by the non-local nature of the problem introduced by the Green function $`\overline{G}_{i,j}`$, Eq. (4) in shb03 . Under uniform loading condition, Eq. (5) in shb03  reads $`u_i=u=_j\overline{G}_{i,j}f_j=_j\overline{G}_{i,j}b^2\sigma `$. Here $`b`$ is the lattice constant. If $`u`$, local deformation, and $`\sigma `$, local stress, are to be independent of size, $`L_x`$ and $`L_y`$, we must have that

$$\underset{j}{}\overline{G}_{i,j}=\mathrm{constant}.$$
(1)
Since $`\overline{G}_{i,j}=\overline{G}_{ij}`$, we may for estetic reasons write $`_j\overline{G}_{i,j}=\mathrm{constant}=_{i,j}\overline{G}_{i,j}/(L_x\times L_y)`$, as there is no dependency on the index $`i`$ in Eq. (1). This was the way we chose to present it in shb03 . We point out again that both indices in the double sum $`_{i,j}`$ run over $`L_x\times L_y`$ sites. We demonstrate the correctness of the rescaling in Fig. 2, where we show the collapse of the loading curves obtained for different system sizes after using Eq. (1). Only when the elastic properties have been rescaled as just described, one may proceed to use finite size scaling as done in shb03 .

We end this Reply by pointing out that Ramanathan and Fisher rf98  measured numerically $`\nu =1.52\pm 0.02`$ using a very different model, which by construction cannot produce fractal fracture lines. This is in complete agreement with our model, which gave $`\nu =1.54`$. The fracture roughness exponents, however, are very different.


text:
# On the structure and spectrum of classical two-dimensional clusters with a logarithmic interaction potential

## I Introduction

In recent years, classical clusters with a finite number of particles moving in a two-dimensional plane were studied for different kinds of confinement and interaction potentials. Bedanov et al. bedanov ; schweigert  studied the classical system with a finite number of particles interacting through a repulsive $`1/r`$ potential and moving in a two-dimensional plane, confined by a parabolic potential. The $`1/r`$ potential is the Coulomb potential in a three-dimensional world. This theoretical system models experimental realizations, such as electrons on the surface of liquid helium helium , electrons in quantum dots in high magnetic fields reimann , colloidal suspensions colloidal  and confined plasma crystals plasma . A detailed investigation of the structure bedanov  and the spectral properties schweigert  of these 2D clusters was carried out. Recently, ground state and metastable configurations and/or spectral properties of 2D clusters with a parabolic confinement potential but different interaction potentials were investigated (see Ref. bolton ; bedanov ; schweigert ; minghuitransition ; koulakov ; lini  for a $`1/r`$ potential, Ref. lini ; lozovik  for a dipole interaction, Ref. lini ; lozovik  for a logarithmic potential, Ref. lini ; candido  for screened Coulomb interaction, …) and with a $`1/r`$ Coulomb interaction potential but with different confinement potentials (see Ref. farias  for a Coulomb confinement potential, Ref. minghui  for a hard wall potential).

Laughlin laughlin  has shown that for strongly interacting electrons in a magnetic field at certain fractional fillings $`(1/3,1/5,\mathrm{})`$ of the Landau levels, the maxima in the pair correlation function correspond to the equilibrium positions in a classical one dimensional plasma where the electrons interact via a repulsive logarithmic potential. The idea was extended to quantum dots manninen  with few electrons as well as bosons and to rings viefers  with a few electrons. It was shown that the electrons actually get localized at the positions where the pair correlation function peaks and the excitations of the system can be understood as rotational modes of the center of mass and vibrational modes. The eigenenergies of the system can be easily determined from the frequencies of the classical normal modes to a high degree of accuracy.

In this paper we consider systems with various confinement potentials containing particles interacting through a repulsive logarithmic interaction. As the solution of the 2D Poisson equation yields a logarithmic potential, one can say that the Coulomb potential in a 2D world is logarithmic. Vortices in a film of liquid Helium interact through a logarithmic potential vorthelium . Also vortices in a type II superconducting 2D film, for a low concentration of vortices, is expected to interact through a logarithmic potential sup . The logarithmic interaction between vortices was used to study the stable vortex configurations in a disk shaped superconductor supvort . Since the vortices are infinitely long in the $`z`$-direction, the symmetry makes it effectively 2D. Hence the model may have some experimental realizations and it may be possible to verify experimentally the predictions made here.

In Section II we describe our model systems. The results for the ground state configurations are shown in Section III, and for the spectra in Section IV. Section V tests the sensitivity of the results by looking at the spectrum of vortices in a type II superconductor, which interact through a logarithmic potential for a large penetration depth. Our conclusions are given in Section V.

## II Model and numerical approach

The Hamitonian of a 2D system of $`N`$ charged particles in a $`r^n`$ confinement potential and interacting through a repulsive logarithmic potential is given by

$$H=\underset{i=1}{\overset{N}{}}\frac{1}{2}m\omega _0^2R^2\left(\frac{r_i}{R}\right)^n+\underset{i>j}{\overset{N}{}}V(\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|),$$
(1)
where $`m`$ is the mass of the particle, $`\omega _0`$ the radial confinement frequency, and $`\stackrel{}{r}_i=(x_i,y_i)`$ the position of the $`i`$th particle with $`r_i|\stackrel{}{r}_i|`$. A hard wall confinement is obtained for $`n\mathrm{}`$ and in this case $`R`$ equals the radius of the hard wall. The interaction potential is taken to be logarithmic:

$$V(r)=\beta \mathrm{ln}(r/R).$$
(2)
We can write the Hamiltonian in a dimensionless form if we express the coordinates and energy in the following units: $`r^{}=\beta ^{1/n}\alpha ^{1/n}R^{(n2)/n},E^{}=\beta `$, with $`\alpha =\frac{1}{2}m\omega _0^2`$. The dimensionless Hamiltonian is given by

$$H=\underset{i=1}{\overset{N}{}}r_i^n\underset{i>j}{\overset{N}{}}\mathrm{ln}\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|.$$
(3)
In the limit of a hard wall confinement, the lengthunit becomes $`r^{}R`$.

We also present results for the actual interaction potential for vortices in a type II superconductor for $`\lambda \mathrm{}`$ as well as for smaller values of $`\lambda `$ (with $`\lambda `$ the penetration depth). Vortices in a type II superconductor interact with a potential tinkham

$$V_s(r)=\beta K_0\left(\frac{r}{\lambda }\right),$$
(4)
with $`K_0`$ the zero-order Hankel function with imaginary argument, $`\beta =\mathrm{\Phi }_0/(2\pi \lambda ^2)`$ and $`\mathrm{\Phi }_0=hc/2e`$ the flux quantum. In the limit of $`\lambda \mathrm{}`$, this potential becomes logarithmic:

$$V_s(r)\beta \mathrm{ln}r,\text{ for }\lambda \mathrm{}\text{ (up to a constant)}.$$
(5)
With the above units, the dimensionless form of this interaction potential is

$$V_s(r)=K_0\left(\frac{r}{\lambda }\right),$$
(6)
where now also $`\lambda `$ is expressed in the unit $`r^{}`$.

The numerical method to obtain the ground state configuration is based on the Monte Carlo simulation technique supplemented with the Newton method in order to increase the accuracy of the energy. The latter technique is outlined and compared with the Monte Carlo technique in Ref. schweigert . The eigenmode frequencies are obtained from the eigenvalues of the dynamical matrix

$$E_{\alpha \beta ,ij}=\frac{^2E}{r_{\alpha ,i}r_{\beta ,j}}|_{r_{\alpha ,i}=r_{\alpha ,i}^n},$$
(7)
where $`\left\{r_{\alpha ,i}^n\right\}`$ is the ground state configuration. The eigenvalues of the dynamical matrix are the squared eigenfrequencies of the system. The eigenfrequencies are expressed in the unit $`\omega ^{}=\sqrt{E^{}/r^2/m}`$.

## III Ground state configurations

In this section we will discuss the numerically obtained ground state configurations for different powers of the confinement potential and show how they can be understood from Earnshaw’s theorem. Fig. 1 shows as an example the ground state configurations for 40 particles for (a) a hard wall confinement (i.e. $`n\mathrm{}`$), (b) a confinement with $`n=3`$ and finally (c) a parabolic confinement (i.e. $`n=2`$). One can see that in the case of hard wall confinement all particles are situated on one ring at the boundary of the potential. This is true for any number of particles in the system, also for clusters with many particles. When the confinement differs from the hard wall confinement, particles can be situated in the central region. Notice that the density of particles for $`n=3`$ is much larger at the edge than in the central region, while for the parabolic confinement the density is uniform.

These observations can be understood from Earnshaw’s theorem earnshaw  in two dimensions. Since in the present case the inter-particle potential is logarithmic we can apply the 2D version of Gauss’s theorem

$$_C𝐄𝑑𝐫_S\rho (𝐫)𝑑S,$$
(8)
where $`S`$ is the surface enclosed by $`C`$, and $`\rho (𝐫)`$ is the enclosed charge. Now it is easy to show that charged particles on a 2D surface can not have a stable static equilibrium when they interact through a logarithmic potential (i.e. the 2D Coulomb potential). If we take a small circle anywhere in the 2D region with no enclosed charges, then, since the line integral of the field along the circle has to be zero, the number of outgoing lines of force are equal to the number of incoming ones. As long as there are outgoing lines of force, a charge left inside this circle can lower its energy by moving along it. Thus in case of hard wall confinement, particles can not be in a stable static equilibrium if they are not in contact with the hard wall potential. Our simulations indeed show that all particles are pushed on one ring situated at the boundary of the hard wall. This is different from the system with a hard wall confinement in which the particles interact through a $`1/r`$ potential. In Ref. bedanov  it was shown that stable multiple ring structures can be found if a sufficiently large number of particles are present in the 2D system.

As shown above the logarithmic interaction potential alone can not result in a stable multi-ring configuration. It is only possible if the hard wall confinement is changed into some soft wall $`r^n`$ confinement potential (as shown in Figs. 1(b) and (c)). In this case, when a particle in the central region is moved in the outward radial direction, the restoring force is given by the confinement potential alone. This must balance the nonrestoring forces provided by the Coulomb repulsion with the other particles. For the parabolic confinement case, this results in a uniform density as the restoring force is everywhere the same.

## IV The eigenmode spectrum

Also the eigenmode spectrum for systems with a logarithmic interaction potential differs from all the previously studied cases. Again the results can be understood from Earnshaw’s theorem. The eigenfrequencies of the system with logarithmic interaction for a $`n=3`$ confinement potential are, as an example for a soft wall confinement, shown in Fig. 2(a) for $`N=2`$ up to 50. Notice that there is a pronounced highest frequency branch. We checked the corresponding eigenmodes and they all show the same behaviour: all outer particles move in the radial direction and have a large amplitude. An example of these breathing-like modes is shown in Fig. 3(a) for $`N=40`$. It is interesting to look now at the eigenfrequencies in the case of a parabolic confinement potential, which is shown in Fig. 2(b). Again the mode corresponding to the highest frequency is the breathing mode (shown in Fig. 3(b) for $`N=40`$). But, as the frequency of the breathing mode in the parabolic case is independent of the number of particles (and given by $`\omega _{\text{breathing}}=\omega _{\text{max}}=2`$ thomson ), there is an upper bound for the frequencies.

For confinement potentials with $`n>2`$ the potential is steeper than the parabolic confinement at large distances from the center. So for the confinement with $`n>2`$ a mode with larger velocities at the boundary than at the center overtakes the breathing mode to be the highest mode. But in all these modes there is a simultaneous contraction of a large number of particles so we can speak of breathing-like modes.

Normal modes are effectively a motion of free particles moving in a resultant potential created by the other particles and the confinement potential. So whenever we try to move a particle from its static equilibrium position, there are restoring forces that bring it back. Then, addition of extra particles complicates the effective potential in which the individual particles are moving and creates higher normal modes. This potential can become steeper and steeper as more particles are added. This results in higher and higher normal modes. But this argument breaks down for the logarithmic interaction for which Gauss’s theorem can be applied: as already shown in previous section, the confinement potential is the only restoring force in the outward radial direction and forms an upper bound to the restoring force. Addition of new particles does not create a more complicated effective potential but it can only result in a slight increase in the frequency of the breathing-like mode due to the increase of the nonrestoring Coulomb force on the outer particles.

## V Vortices in superconductors

As mentioned in the Introduction, the logarithmic potential is sometimes a good model for the interaction between vortices in a type II superconductor. In order to check how important small deviations from this logarithmic potential are, we also studied the spectrum of vortices in a type II superconductor for a more generally valid interaction potential. In Ref. tinkham  it is shown that the interaction potential between vortices in a type II superconductor is given by Eq. (4) if the penetration length $`\lambda `$ the coherence length. In the limit of $`\lambda \mathrm{}`$, this potential reduces to the logarithmic potential studied in the previous sections.

The spectrum for the parabolically confined system is shown in Fig 4 for penetration lengths $`\lambda =1,10`$ and 50. One can clearly see that for $`\lambda =1`$ no upper bound frequency is found, while it is for $`\lambda =50`$. That the logarithmic limit is not at all reached for $`\lambda =1`$ is also seen in Fig. 5 which shows the eigenmode corresponding to the largest frequency for $`N=10`$. As shown above, this eigenmode corresponds with the breathing mode in case of a logarithmic interaction. This is not at all the case for $`\lambda =1`$ and 10 (see Figs.5(a) and (b)), but for $`\lambda =50`$ it does. Therefore, transition to the logarithmic-type of interaction takes place between $`\lambda =10`$ and 50.

It is also interesting to look at the eigenmodes corresponding to the highest frequency for $`N=20`$ (see Figs. 6(a-c)). Now the modes for $`\lambda =1,10`$ and 50 look all similar. However, they differ from the breathing mode: the middle and outer ring move in the opposite direction. From the spectrum one can observe that for $`N>10`$ there are two highest frequencies which lie very close together. These highest but one frequencies are therefore plot in Fig. 6(d-f). Notice that for $`\lambda =10`$ and 50 this mode indeed corresponds with the breathing mode.

## VI Conclusions

We investigated the ground state structure and the spectral properties of a classical 2D cluster with a finite number of particles interacting through a repulsive logarithmic potential. The logarithmic potential is the Coulomb potential in a 2D world. It was shown that for a hard wall confinement all particles are situated on one ring at the boundary of the potential, as a consequence of Earnshaw’s theorem in 2D. Multi-ring configurations exist for soft confinement potentials in which case only the external confinement delivers the restoring force. The particles are uniformly distributed for a parabolic confinement.

We also found that, independently of the number of particles in the cluster, all eigenfrequencies lie below the frequency of a particular mode, namely the breathing-like mode. This is again a consequence of Earnshaw’s theorem because adding extra particles does not create extra minima in the potential nor makes it more complex. In the case of a parabolic confinement the frequency of this breathing mode is independent from $`N`$, which results in an upper bound for the frequencies.

It is also shown that for vortices in a type II superconductor (for which the penetration length $`\lambda `$ coherence length) the transition to the logarithmic-type of interaction takes place between $`10<\lambda <50`$. It is not always the highest mode anymore which corresponds with the breathing mode, it can also be the highest but one mode.

###### Acknowledgements.

This work is supported by the Flemish Science Foundation.


text:
# A close look at Jeckelmann’s weak-coupling results

In his comment, Jeckelmann claims that our results  are inaccurate and our criticism of his work is groundless, but unfortunately, in contrast to his assertion, our conclusion is firmly grounded on a well-established fact, independent of our numerical calculation. The fact is that at the weak-coupling limit, $`U/V_c`$ must approach to 2. Jeckelmann’s results deviate from 2 in the weak-coupling limit, and therefore are troublesome. In his Comment, while avoiding this basic fact, Jeckelmann does not provide new information about his calculation but instead he tries to create more confusions by making the same mistakes as his previous paper did, which is unhelpful to the undergoing debate on the phase diagram.

On the other hand, our results agree with both the weak- and strong-coupling limit results, thus our conclusions in our paper  remain unchanged. In the following, we will investigate and reply to his comments in detail.

First, Jeckelmann completely ignored the first density-matrix renormalization group (DMRG) calculation  on the phase transition in the extended Hubbard model. As shown in our previous papers, since $`V_c`$ is so close to $`U/2`$, a $`V_c`$ versus $`U`$ plot or a variation of it is very insensitive to the phase transition boundary. Jeckelmann is not able to show a critical figure ($`U/V_c`$ versus $`U`$) that we plotted in our paper. This is the key mistake that he made before and now in this comment again. His Fig. 1 has the same problem as his old figure has, which hides the important difference between his results and the weak-coupling limit results though his Figure 1 already shows the deviation of his results (open diamonds) from the weak-coupling limit. This unfortunately misleads him to conclude his results are in quantitative agreement with other works, but it is interesting to note that his paper has already been commented by the same group of researchers that Jeckelmann claims his results agree with, (for details see Ref.). If he plotted his data on a $`U/V_c`$ versus $`U`$ plot, he would immediately see a deviation from the weak-coupling limit results. To be clear, in Fig. 1 we show such a deviation of his results, which we will come back to below.

Second, Jeckelmann repeatedly tries to claim his results agree quantitatively with the formal quantum Monte-Carlo results, which is untrue. The reason is that he has eight (8) sets of data while three-four (3-4) sets of data are available in the work by Sengupta. We pointed out his problem in our paper. We do not understand why the author could fail to see the difference between them. Instead of directly addressing our concerns, he questions our calculation.

Third, all the numerical methods have their limitations. We have been familiar with those references for a long time though Jeckelmann failed to notice our previous paper. The results by Cannon et al. were obtained in small clusters. Naturally, the finite-size effect is unavoidable. Hirsch’s MC results were obtained at one fixed temperature $`\tau `$. Cannon et al. showed that the MC results depend on $`\tau `$, thus it is necessary to perform $`\tau 0`$ extrapolation in order to get accurate results. Nakamura pointed out that the direct CDW-SDW level crossing point has a large size-effect for all regions. As already mentioned in our paper, our calculation, not different from all the other DMRG calculations, has several errors of different kinds: truncation error, finite-size effect, $`U`$-dependent error, numerical truncation error and many others. That is why we are very careful when we presented our results. Our results are obtained for system of 40 sites, truncation of 400 states, without any extrapolation and within the infinite-system algorithm. Our results are different from those above results since our calculations are done in longer chains, have smaller finite-size effects, do not have problems with temperature, and have a better control over state truncation. Therefore, the results tend to be more accurate. In particular, our results are fully consistent with both strong- and weak-coupling limit results. The reason why our present results are different from our previous results is simply because at present we used larger systems and higher truncation, then the results certainly become more accurate. In fact, in our paper we have systematically investigated this change (see Fig. 2 of our paper ), while there is no such result in Jeckelmann’s paper.

Fourth, Jeckelmann mentioned the infinite-system algorithm is known to give incorrect results occasionally, but this does not necessarily mean the results by the infinite-system algorithm are wrong. Without providing any evidence, he asserts the failure of the standard infinite-system DMRG algorithm for the present problem. He goes on to provide one example in his reference 9 of his comment to justify his argument. We looked at the reference and found that paper is about a Bosonic Hubbard model not about the extended Hubbard model. While it is true that the finite-system algorithm is more accurate, what is surprising is that his results from an algorithm that he claims ‘exact’ are not consistent with the weak-coupling limit results.

Fifth, the contradiction of his results with the weak-coupling limit is truly troublesome. Now he admits that his investigation focused on the intermediate coupling regime, but nowhere in his original letter mentions this, nor does this help him explain away the contradiction. As shown in our previous paper, the trend in his results is incorrect and deviates from $`U/V_c=2`$ substantially (see Fig. 1). Jeckelmann now bypasses this comparison by saying that the $`U0`$ of $`U/V_c`$ cannot be determined using the data obtained in his work. In order to justify his claim, he cited the paper by Nakamura that the weak-coupling limit is recovered only for $`U`$ much smaller than $`2t`$. We checked Nakamura’s paper. To our surprise, in sharp contrast to Jeckelmann’s claim, Nakamura showed that the intermediate-coupling is at $`U=4t`$ and smaller than $`4t`$ the results asymptotically approach to the weak coupling limit (see Fig. 12 in Nakamura’s paper). This further demonstrates Jeckelmann’s results are troubling. It is our hope that the author should directly and seriously address those issues since the evidence is very clear.

Finally, for some reason, Jeckelmann apparently misunderstood our statement in our paper. When we said that we could not find enough information, we mean that we could not find enough information in his paper not in the cited references. By contrast, we are pretty familiar with all those references for a long time. In his paper there is no information about system size, truncation, details about the extrapolation and how to determine the phase boundaries. That is why we were not able to compare our results with his results, which is a motivation of our study. Once again, his new Fig. 1 avoids the critical comparison among different calculations and does not provide a useful checkup of his data. Finally, we would like to take this opportunity to point out that in scientific community words like ‘ludicrous’ should be avoided and it is our belief that a little criticism on one’s work is very healthy for the author himself/herself as well as for the entire community.

In conclusion, we emphasize that none of our criticism on his paper as well as our main conclusion has been changed. We kindly suggest that the author check his data and address the criticisms more directly and throughly.


text:
# References

Cosmic Hydrogen Was Significantly Neutral

a Billion Years After the Big Bang

$``$ School of Physics, The University of Melbourne, Parkville, Vic 3010, Australia
$``$ Astronomy Dept., Harvard University, 60 Garden Street, Cambridge, MA 02138, USA

The ionization fraction of cosmic hydrogen, left over from the big bang, provides crucial fossil evidence for when the first stars and quasar black holes formed in the infant universe. Spectra of the two most distant quasars known1  show nearly complete absorption of photons with wavelengths shorter than the Ly$`\alpha `$ transition of neutral hydrogen, indicating that hydrogen in the intergalactic medium (IGM) had not been completely ionized at a redshift $`z6.3`$, about a billion years after the big bang. Here we show that the radii of influence of ionizing radiation from these quasars imply that the surrounding IGM had a neutral hydrogen fraction of tens of percent prior to the quasar activity, much higher than previous lower limits 1 ,2  of $`0.1\%`$. When combined with the recent inference of a large cumulative optical depth to electron scattering after cosmological recombination from the WMAP data3 , our result suggests the existence of a second peak in the mean ionization history, potentially due to an early formation episode of the first stars.

The detection of a Ly$`\alpha `$ emitting galaxy at $`z6.56`$ led to the claim that the IGM must already be highly ionized at that redshift4 . However, it was later shown5  that the presence of a small ionized region around the galaxy and broadening of the Ly$`\alpha `$ line, would allow detection of the line feature in a fully neutral IGM. Luminous quasars offer an alternative probe of the IGM and have higher luminosities, allowing acquisition of higher quality spectra. The spectra6 ,7  of SDSS J1030+0524 at $`z=6.28`$ and SDSS J1148+5251 at $`z=6.41`$ show transmitted flux down to wavelengths corresponding to the Ly$`\alpha `$ resonance of hydrogen at $`z6.20`$ and $`z6.32`$ respectively, even though no flux is detected at somewhat greater distances from the quasars7 . These redshift displacements correspond to velocity shifts of order $`3300\mathrm{km}\mathrm{s}^1`$, and imply ionized regions surrounding the quasars with observed physical radii of $`R_{\mathrm{obs}}=4.5`$ Mpc and $`4.7`$ Mpc, respectively. The inference is robust since high ionization broad-lines typically have velocity offsets from the true quasar redshift8  of only $`1000\mathrm{km}\mathrm{s}^1`$.

We model the evolution of a fully ionized region (the so-called Strömgren sphere) around a quasar that emits UV photons into a partially neutral IGM. Prior to the overlap phase of ionized hydrogen (H II ) regions that signals the end of the reionization epoch9 , galaxies generated their own small ionized regions filling a fraction of the IGM volume around the quasar. We define the neutral fraction, $`x_{\mathrm{HI}}`$, to be the mean fraction of hydrogen atoms (in a representative volume of the low-density IGM) that remain neutral prior to the quasar activity, and assume $`x_{\mathrm{HI}}`$ outside the Strömgren sphere to remain constant during its expansion. Since the luminosity of a bright quasar is much larger than the stellar luminosity of neighboring galaxies, a fully ionized quasar bubble expands into a partially ionized IGM composed of isolated H II regions.

Even for the deepest available spectra, with lower limits on the Ly$`\alpha `$ (Gunn-Peterson10 ) optical depth of $`\tau _{\mathrm{GP}}26`$, the IGM need only have a neutral fraction of $`10^3`$ in order to produce the observed absorption blueward of Ly$`\alpha `$1 . However, the value of $`R_{\mathrm{obs}}4.5`$Mpc may be used to constrain larger values of the neutral fraction in the range $`0.001x_{\mathrm{HI}}1.0`$. In the early expansion phase, the dependence of the physical radius of the Strömgren sphere, $`R_\mathrm{p}`$, on the quasar age1 ,11 , $`t_{\mathrm{age}}`$, may be crudely approximated as $`R_\mathrm{p}7x_{\mathrm{HI}}^{1/3}\left[t_{\mathrm{age}}/10^7\mathrm{yr}\right]^{1/3}`$Mpc, given the production rate of ionizing photons that is characteristic of SDSS J1148+5251 and SDSS J1030+0524. Clearly, the Strömgren sphere grows larger when embedded in an IGM with a small neutral fraction. Defining $`R_{\mathrm{max}}`$ to be the maximum radius achieved in a fully neutral IGM by the end of the quasar lifetime, $`t_{\mathrm{lt}}`$, the resulting a priori probability distribution is, $`dP/dR_\mathrm{p}|_{x_{\mathrm{HI}}}=3x_{\mathrm{HI}}R_\mathrm{p}^2/R_{\mathrm{max}}^3`$ for $`R_\mathrm{p}<R_{\mathrm{max}}x_{\mathrm{HI}}^{1/3}`$. In addition, the neutral fraction must be smaller than $`x_{\mathrm{max}}=\mathrm{min}(\left[t_{\mathrm{lt}}/10^7\mathrm{years}\right]\left[R_{\mathrm{obs}}/7\mathrm{M}\mathrm{p}\mathrm{c}\right]^3,1)`$ to allow the growth of the sphere to $`R_{\mathrm{obs}}`$ within $`t_{\mathrm{lt}}`$. Using a flat logarithmic prior for $`x_{\mathrm{HI}}`$, we find the cumulative a posteriori probability distribution, $`P(<x_{\mathrm{HI}})=\mathrm{min}(x_{\mathrm{HI}},x_{\mathrm{max}})/x_{\mathrm{max}}`$. Since estimates of quasar lifetimes12  are bracketed in the range $`10^610^8`$ years, we get $`x_{\mathrm{max}}1`$ and conclude that $`x_{\mathrm{HI}}0.1`$ ($`0.01`$) with 90% (99%) confidence.

To better quantify this simple argument, we write the general equation for the relativistic expansion of the co-moving radius \[$`R=(1+z)R_\mathrm{p}`$\] of the quasar H II region1  in an IGM with a neutral filling fraction $`x_{\mathrm{HI}}`$ (fixed by other ionizing sources),

$$\frac{dR}{dt}=c(1+z)\left[\frac{F_\gamma \dot{N}_{\mathrm{ion}}\alpha _\mathrm{B}CF_\mathrm{m}x_{\mathrm{HI}}\left(\overline{n}_0^\mathrm{H}\right)^2\left(1+z\right)^3\left(\frac{4\pi }{3}R^3\right)}{F_\gamma \dot{N}_{\mathrm{ion}}+4\pi R^2\left(1+z\right)cF_\mathrm{m}x_{\mathrm{HI}}\overline{n}_0^\mathrm{H}}\right],$$
(1)
where $`c`$ is the speed of light, $`\overline{n}_0^\mathrm{H}`$ is the mean number density of protons at $`z=0`$, $`\alpha _\mathrm{B}=2.6\times 10^{13}`$cm3s-1 is the case-B recombination coefficient at the characteristic temperature of $`10^4`$K, and $`\dot{N}_{\mathrm{ion}}`$ is the rate of ionizing photons crossing a shell at the radius of the H II region at time $`t`$. We use the distribution derived from numerical simulations for the over-densities in gas clumps13 , and calculate the mean free path $`d(\mathrm{\Delta }_\mathrm{c})`$ for ionizing photons13 ,14 . Following Barkana & Loeb14 , we then find the critical overdensity $`(\mathrm{\Delta }_\mathrm{c})`$ at which a fraction $`F_\gamma 0.5=\mathrm{exp}\left[R_\mathrm{p}/d(\mathrm{\Delta }_\mathrm{c})\right]`$, of photons emitted do not encounter an overdensity larger than $`\mathrm{\Delta }_\mathrm{c}`$ within the H II region. We also compute the mass fraction $`F_\mathrm{m}`$ ($`1`$) of gas within $`R_\mathrm{p}`$ that is at over-densities lower than $`\mathrm{\Delta }_\mathrm{c}`$. Finally, we calculate the clumping factor in the ionized regions, $`C(R)\mathrm{\Delta }^2/\mathrm{\Delta }^2`$, where the angular brackets denote an average over all regions with $`\mathrm{\Delta }<\mathrm{\Delta }_\mathrm{c}`$. For $`x_{\mathrm{HI}}=1`$ and $`dR_\mathrm{p}/dtc`$, equation (1) reduces to its well-known form11 ,15 -17 . The expansion of the quasar Strömgren sphere would change in an overdense region of the IGM. However the density contrast due to infall has been shown14  to be small ($`1\%`$) on scales of several Mpc around the massive hosts of the bright SDSS quasars at $`z6`$. Numerical experiments showed our results to be insensitive to infall. The emission rate of ionizing photons $`\dot{N}_{\mathrm{ion}}`$ in equation (1) is computed at $`t^{}=tt_{\mathrm{delay}}`$ to account for the finite light travel time between the source and the ionization front. The delay $`t_{\mathrm{delay}}`$ is derived from the relation $`R=_{t_Rt_{\mathrm{delay}}}^{t_R}c𝑑t^{}\left[1+z(t^{})\right]`$, where $`t_R`$ is the time when the photon crosses the co-moving radius $`R`$. The use of the Telfer et al.18  spectral template implies $`\dot{N}_{\mathrm{ion}}=6.5\times 10^{57}`$s-1 for SDSS J1030+0524 and $`\dot{N}_{\mathrm{ion}}=10.0\times 10^{57}`$s-1 for SDSS J1148+5251, while the template of Elvis et al.19  implies lower values of $`\dot{N}_{\mathrm{ion}}1.3\times 10^{57}`$s-1 and $`\dot{N}_{\mathrm{ion}}2.0\times 10^{57}`$s-1, respectively1 .

Our model for the evolution of the Strömgren sphere includes the quasar formation history. We compute the time dependent ionizing flux governed by black-hole growth through accretion and mergers. Following the observational inference that the relation between bulge velocity dispersion and black hole mass $`M_{\mathrm{bh}}`$ does not evolve with redshift20 , we extrapolate the present-day $`M_{\mathrm{bh}}`$$`M_{\mathrm{halo}}`$ relation21  using the dependence of virial velocity on redshift9  to obtain $`M_{\mathrm{halo}}=1.5\times 10^{12}M_{}\left(M_{\mathrm{bh}}/10^9M_{}\right)^{3/5}\left[(1+z)/7\right]^{3/2}`$. The luminosities of the SDSS quasars imply black-hole masses of $`2\times 10^9`$M accreting at their Eddington rate6 , leading to inferred host dark-matter halos of $`M_{\mathrm{halo}}2\times 10^{12}`$M. The inference of such a massive halo at this early epoch is supported by the signature of gas infall in the spectra of high redshift quasars14 , and the inference of a large molecular mass and velocity width in the host galaxy of SDSS J1148+52517 .

We begin with two dark-matter halos of comparable mass $`M_1`$ and $`M_2`$ that merge to form the host dark-matter halo with a mass $`M_{\mathrm{halo}}=(M_1+M_2)`$ at a redshift corresponding to one quasar lifetime prior the observed redshift, and construct merger trees for each of the progenitor halos using the algorithm described by Volonteri et al.22 . We use the $`M_{\mathrm{bh}}`$$`M_{\mathrm{halo}}`$ relation to find the mass of the black-hole at the center of each dark-matter halo in the merger tree. When there is a major merger (having a progenitor mass ratio $`<2`$), we assume that the black-holes merge and that the accretion of the mass deficit relative to the $`M_{\mathrm{bh}}`$$`M_{\mathrm{halo}}`$ relation produces the limiting Eddington luminosity with a radiative efficiency of $`ϵ=10\%`$. The resulting quasar lifetime is $`t_{\mathrm{lt}}=4\times 10^7(ϵ/0.1)\mathrm{ln}\left[M_{\mathrm{halo}}^{5/3}/(M_1^{5/3}+M_2^{5/3})\right]\text{ years}`$. The ionizing photon emission rate $`\dot{N}(t)`$ is taken to have an exponential time dependence, $`\dot{N}_{\mathrm{ion}}(t)=\dot{N}_0e^{t/t_{\mathrm{lt}}}`$, with $`\dot{N}_0M_{\mathrm{bh}}`$. We assume that all quasar episodes in the merger tree contribute to the volume of the observed Strömgren sphere, which is centered on the most massive halo in the tree. A similar prescription for quasar activity was shown to be consistent with the observed number counts of high redshift quasars22 . For a major merger resulting in the observed quasar activity, the lifetime given by this approach is $`10^7`$ years. This compares favorably with a variety observational estimates of quasar lifetimes23 -26 . To cover the full range of uncertainty we consider modifications to our fiducial model in which the quasar lifetime is multiplied by a factor $`f_{\mathrm{lt}}`$ between 0.1–10, resulting in lifetimes for the observed quasars of $`f_{\mathrm{lt}}t_{\mathrm{lt}}=10^6`$$`10^8`$ years.

We have produced 300 realizations of the merger tree, and computed the evolution of the Strömgren sphere in each case for the full range of neutral fractions $`0.001<x_{\mathrm{HI}}<\mathrm{\hspace{0.33em}1.0}`$, allowed by the Gunn-Peterson optical depth1 . We find that quasar activity associated with the hierarchical build up of the host galaxy produces an H II region with a typical size of $`2x_{\mathrm{HI}}^{1/3}`$Mpc, comparable to the observed radii if $`x_{\mathrm{HI}}`$ is of order unity. Figure 1 shows the conditional a priori probability distributions of $`R_\mathrm{p}`$ for quasars like SDSS J1030+0524, assuming $`x_{\mathrm{HI}}=1`$ and two different values of $`f_{\mathrm{lt}}`$. The observed radii are consistent with lifetimes close to the fiducial case of $`f_{\mathrm{lt}}=1`$. The plotted distribution of $`R_\mathrm{p}`$ implies that the hierarchical evolution of early quasars leads to a dearth of small Strömgren sphere radii for the observed SDSS quasars. This is in contrast to the simplified monolithic model of quasar formation16 ,1 , for which the H II region has zero volume at the time when the quasar turns on \[$`R(t_{\mathrm{age}}=0)=0`$\] and the distribution extends down to $`R_\mathrm{p}=0`$.

To constrain the neutral fraction, we have computed the likelihood of observing $`R_{\mathrm{obs}}=4.5`$Mpc around SDSS J1030+0524 and $`R_{\mathrm{obs}}=4.7`$Mpc around SDSS J1148+5251 as a function of $`x_{\mathrm{HI}}`$ and $`f_{\mathrm{lt}}`$. The results are plotted in Figure 2. The upper panels show the locus of most likely values (thick line) as well as likelihood contours at 0.1 of the peak value (dashed lines). The fiducial lifetime ($`f_{\mathrm{lt}}`$) favors $`x_{\mathrm{HI}}1`$, while with $`f_{\mathrm{lt}}1`$ the distributions for $`R_\mathrm{p}`$ lie substantially below the observed value of $`4.5`$Mpc, making smaller values of $`x_{\mathrm{HI}}`$ more likely. Extrapolation of the most likely contour for the Elvis et al.19  template yields $`f_{\mathrm{lt}}2x_{\mathrm{HI}}`$. This implies that $`x_{\mathrm{HI}}10^3`$ would require a lifetime as short as $`2\times 10^4`$ years, which is ruled out by variability properties of quasars in SDSS26 . The a posteriori probability distributions for $`x_{\mathrm{HI}}`$ given $`f_{\mathrm{lt}}`$ are plotted in the lower panels of Figure 2 and robustly yield the constraint $`x_{\mathrm{HI}}>\mathrm{\hspace{0.33em}0.01}`$. For $`f_{\mathrm{lt}}0.3`$, we find $`x_{\mathrm{HI}}>\mathrm{\hspace{0.33em}0.1}`$ and $`x_{\mathrm{HI}}>\mathrm{\hspace{0.33em}0.4}`$ at the 90% level for the Elvis et al.19  and Telfer et al.18  spectra, respectively. The fiducial model with $`f_{\mathrm{lt}}=1`$ yields corresponding constraints of $`x_{\mathrm{HI}}>\mathrm{\hspace{0.33em}0.3}`$ and $`x_{\mathrm{HI}}>\mathrm{\hspace{0.33em}0.6}`$.

The inference of a large neutral fraction at $`z6.3`$ presents a challenge to theories of cosmological reionization when combined with the large optical depth3  to electron scattering after cosmological recombination, $`\tau _{\mathrm{es}}=0.17\pm 0.04`$. Consider a toy model in which the universe was partially reionized at a high redshift $`z_{\mathrm{reion},\mathrm{I}}`$, leaving a neutral fraction $`x_{\mathrm{HI}}`$ until complete reionization was reached at $`z_{\mathrm{reion},\mathrm{II}}6.25`$. The optical depth is then $`\tau _{\mathrm{es}}=0.04+0.002(1x_{\mathrm{HI}})\left[(1+z_{\mathrm{reion},\mathrm{I}})^{3/2}19.5\right]`$. If the IGM ionization fraction increased monotonically, then given $`\tau _{\mathrm{es}}>0.13`$, the universe needed to be reionized earlier than $`z_{\mathrm{reion},\mathrm{I}}18`$ or $`z_{\mathrm{reion},\mathrm{I}}24`$ assuming the Elvis et al.19  and Telfer et al.18  spectra, respectively. If $`\tau _{\mathrm{es}}=0.17`$ and $`x_{\mathrm{HI}}>0.7`$, as implied by the Telfer et al.18  spectrum, then a monotonic reionization history requires significant reionization at an implausibly high redshift ($`z>\mathrm{\hspace{0.33em}30}`$). In this case, a more plausible alternative is a non-monotonic history with an early reionization peak, possibly due to the formation of massive population-III stars27 ,28 .

Correspondence and requests for materials to Abraham Loeb.

ACKNOWLEDGEMENTS. This work was supported in part by grants from ARC, NSF and NASA.