text:
# Simulating IGM Reionization

## 1 Introduction

Recent discoveries of quasars at $`z>5.8`$ (e.g. Fan et al. 2000, 2001) are finally allowing quantitative studies of the high-redshift InterGalactic Medium (IGM) and its reionization history. In particular, the detection of a Gunn-Peterson trough (Gunn & Peterson 1965) in the Keck (Becker et al. 2001) and VLT (Pentericci et al. 2002) spectra of the Sloan Digital Sky Survey quasar SDSS 1030-0524 at $`z=6.28`$ is a clear indication that the universe is approaching the reionization epoch at $`z6`$.

Whatever the exact value of the reionization redshift, $`z_{ion}`$, it is clear that the hydrogen in the IGM is in a highly ionized state at $`z\stackrel{<}{}6`$. The nature of the responsible sources is the subject of a lively debate. Several authors (Madau, Haardt & Rees 1999 and references therein) have claimed that the known populations of quasars and galaxies provide $`10`$ times fewer ionizing photons than are necessary to explain the observed IGM ionization level. Thus, additional sources of ionizing photons are required at high redshift, the most promising being early galaxies and quasars. While no evidence for the existence of high-redshift quasars has yet been found, recent observations of temperature and metal abundance in the lowest density regions of the IGM suggest the existence of an early population of pregalactic stellar objects, which may have contributed to the reionization and metal enrichment of the IGM (Cowie & Songaila 1998; Ellison et al. 2000; Schaye et al. 2000; Madau, Ferrara & Rees 2001). For this reason, most theoretical work on IGM reionization has assumed stellar sources.

The study of IGM reionization by primeval stellar sources has been tackled by several authors, both via semi-analytic (e.g. Haiman & Loeb 1997; Valageas & Silk 1999; Miralda-Escudé, Haehnelt & Rees 2000; Cojazzi et al. 2000) and numerical (e.g. Gnedin & Ostriker 1997; Ciardi et al. 2000; Chiu & Ostriker 2000; Gnedin 2000; Razoumov et al. 2002) approaches. Two main ingredients are required for a proper treatment of the reionization process: i) a model of galaxy formation and ii) a reliable treatment of the radiative transfer of ionizing photons.

The commonly accepted scenario for structure formation is of a universe dominated by Cold Dark Matter (CDM). The evolution of dark matter structures in such a universe is now well understood, while the treatment of the physical processes that govern the formation and evolution of luminous objects (e.g. heating/cooling of the gas, star formation, feedback), is still quite uncertain. Despite many applications of hydrodynamical simulations to the structure formation process (e.g. Cen & Ostriker 1993, 2000; Katz, Weinberg & Hernquist 1996; Katz, Hernquist & Weinberg 1999; Blanton et al. 1999; Pearce et al. 2001; Weinberg, Hernquist & Katz 2002), much of our current understanding has come from semi-analytic models, based on simplified physical assumptions. These have proved successful in reproducing many (but not all) observed galaxy properties (e.g. Lacey et al. 1993; Kauffmann, White & Guiderdoni 1993; Heyl et al. 1995; Baugh et al. 1998; Mo, Mao & White 1998; Devriendt et al. 1998; Cole et al. 2000; Somerville, Primack & Faber 2001). Although these models do not provide as detailed information as numerical simulations, they have the advantage of allowing a larger dynamic range and of being fast enough to explore a vast range of parameters. Their main disadvantage is that no information on the spatial distribution of structures is provided. For this reason a new method has recently been developed combining N-body numerical simulations of the dark component of matter, with semi-analytic models which ‘attach’ to it the galaxies and their properties (e.g. Kauffmann, Nusser & Steinmetz 1997; Kauffmann et al. 1999; Benson et al. 2000; Diaferio et al. 2001; Springel et al. 2001, SWTK). The advantage of this method is that, once the cosmology is fixed, the N-body simulation, which provides the overall structure, has to be performed only once, while the semi-analytic model can be run many times to explore different physical parameters describing the galaxy formation process.

A limitation of the above method is related to the resolution of the N-body simulation, and turns out to be critical in simulations aimed at the study of the reionization process. The resolution must be high enough to follow the formation and evolution of the objects responsible for producing the bulk of the ionizing radiation. At the same time, a large simulation volume is required to have a region with “representative” properties and to avoid biases due to cosmic variance on small scales. So far, simulations of reionization have been run in fairly small volumes ($`L7h^1`$ Mpc comoving; Ciardi et al. 2000; Gnedin 2000; Razoumov et al. 2002). These fail to represent correctly the expected large-scale distribution either of the ionizing sources or of the material to be ionized. The minimum mass of objects that contribute substantially to the reionization process must be determined and then, simulations run in the largest possible volume compatible with resolving this mass. This procedure ensures the best possible treatment of the ionizing photon production (Ciardi 2002).

In addition to a reliable model of galaxy formation, an accurate treatment of the propagation of ionizing photons is required to study IGM reionization. The full solution of the seven dimension radiative transfer equation is still well beyond our computational capabilities, and although in some specific cases it is possible to reduce its dimensionality, for the reionization process no spatial symmetry can be invoked. Several authors (e.g. Umemura, Nakamoto & Susa 1999; Razoumov & Scott 1999; Abel, Norman & Madau 1999; Gnedin 2000; Ciardi et al. 2001, CFMR; Gnedin & Abel 2001; Cen 2002; Maselli, Ferrara & Ciardi 2003) have recently devoted their efforts to the development of radiative transfer codes based on a variety of approaches (e.g. ray tracing, Monte Carlo, local optical depth approximation). Given the intrinsic difficulty in solving a high dimensionality equation, most codes aim to solve simplified problems, and cannot be applied for more general use. Moreover, these codes face a speed problem. Most of them are prohibitively slow when applied in a cosmological contest. For this reason, approximate methods have been developed to allow faster calculations (e.g. Gnedin & Abel 2001; Abel & Wandelt 2002).

In this paper we study the IGM reionization process through a combination of high-resolution N-body simulations (to describe the distribution of dark matter and diffuse gas), a semi-analytic model of galaxy formation (to track the sources of ionization) and the Monte Carlo radiative transfer code CRASH (to follow the propagation of ionizing photons in the IGM; CFMR). In Section 2 we describe the numerical simulations of structure formation and in Section 3 the radiative transfer code. In Section 4 we present the results of the calculation and in Section 5 and 6 our discussion and conclusions.

## 2 Numerical Simulations of Structure Formation

As discussed in the Introduction, the mass resolution of the underlying N-body simulation plays a crucial role in modeling reionization. The objects that produce the bulk of the ionizing photons must be resolved. While at $`z>15`$ the main contributors are small mass objects ($`M10^7`$ M), at lower redshift, the radiation from such objects is negligible compared with that from galaxies with total (dark halo and galaxy) masses $`M>10^9`$ M. Thus, our numerical simulations should be able to resolve objects with $`M\stackrel{<}{}10^9`$ M and, at the same time, have the largest possible volume compatible with such mass resolution. Smaller mass objects remain important, however, for the small scale clumping (Ciardi 2002; Ciardi et al. 2000).

We started out with a simulation from the VIRGO consortium ($`479h^1`$ Mpc comoving on a side; Yoshida, Sheth & Diaferio 2001) based on a $`\mathrm{\Lambda }`$CDM cosmology with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`h=0.7`$, $`n=1`$, $`\sigma _8=0.9`$ and $`\mathrm{\Omega }_b=0.04`$. We then selected a “typical” spherical region of diameter $`50h^1`$ Mpc and resimulated it four times with higher mass resolutions within the region and lower resolution far outside it. Particle masses in the high resolution region are $`M_p=6.8\times 10^{10}h^1`$ M, $`4.8\times 10^9h^1`$ M, $`9.5\times 10^8h^1`$ M and $`1.7\times 10^8h^1`$ M for the ‘M0’, ‘M1’, ‘M2’ and ‘M3’ simulations, respectively. A discussion of the effect of mass resolution on our computations is presented in Section 5. All four simulations were prepared with the techniques of SWTK and were carried out with the N-body code GADGET (Springel, Yoshida & White 2001). For the highest mass resolution $`7\times 10^7`$ particles had to be followed. Finally, we have extracted a box of comoving side $`L=20h^1`$ Mpc, which will be used to study the reionization process. A larger box ($`30h^1`$ Mpc) could have been extracted, but would have required a prohibitively long radiative transfer. The choice of a $`20h^1`$ Mpc box allows a reasonably fast calculation in a region of the universe with “representative” properties. The location and mass of the dark matter haloes in the simulation have been determined by means of a friends-of-friends algorithm, while we identify gravitationally bound substructures with SUBFIND (SWTK) and build the merging tree for haloes and subhaloes following the prescription of SWTK. The smallest haloes have masses of $`M10^9`$ M, consistent with the requirements described above.

We then model the galaxy population with the semi-analytic technique of Kauffmann et al. (1999) in the implementation of SWTK. This procedure follows the merging tree of the dark matter haloes extracted from the simulation, and models galaxy formation in the haloes and subhaloes using simple recipes for gas cooling, star formation, feedback from supernovae and galaxy-galaxy merging. We adopt the same parameter values in these recipes as SWTK. At the end of this process we obtain a catalogue of galaxies for each of the 51 simulation outputs, containing for each galaxy, among other quantities, its position, stellar mass and star formation rate (SFR). The model reproduces reasonably well the extinction and incompleteness corrected SFR data of Somerville, Faber & Primack (2001) up to a redshift of 3. For higher redshifts the situation is less clear because of the lack of observational data, but we expect the model to be accurate to within a factor of two. The model is fully consistent with the upper bounds from the maximal extinction correction estimates of Somerville, Faber & Primack (2001).

In order to infer the emission properties of our model galaxies, we have assumed a Salpeter IMF and a time-dependent spectrum typical of metal-free stars (CFMR). This choice is further discussed in Section 5. Of the emitted ionizing photons, only a fraction $`f_{esc}`$ will actually be able to escape into the IGM. Throughout the calculation we adopt $`f_{esc}=5\%`$. This choice assumes that $`f_{esc}`$ is independent of other physical quantities and is the same for every galaxy. Actually, $`f_{esc}`$ may well vary with, e.g., redshift, mass and structure of a galaxy, as well as with the ionizing photon production rate (e.g. Wood & Loeb 1999; Ricotti & Shull 2000; Ciardi, Bianchi & Ferrara 2002). Given the large uncertainty in $`f_{esc}`$, we consider $`5\%`$ merely as a reference value. This choice is further discussed in Section 5.

In order to determine the effects of environment on the reionization process, we will also study how reionization proceeds in a proto-cluster rather than in a field region. For this purpose we use the cluster simulations performed by SWTK. Also in this case numerical simulations are available at different mass resolution, ‘S1’-‘S4’. Runs labelled with the same numbers in the ‘M’ and ‘S’ series have comparable mass resolution.

### 2.1 Simulation Convergence

We now use our series of simulations run to check for convergence. In Fig. 1 the total star formation rate, $`\dot{M}_{}`$, normalized to the total baryonic mass, $`M_b`$, is shown as a function of redshift for the field region simulations (upper panel). The ‘M2’ simulation converges to ‘M3’ at $`z_{conv}5.5`$, while the ‘M0’ simulation converges to higher resolution results only at $`z_{conv}0`$. The lower panel of Fig. 1 shows a similar comparison for the cluster simulations. In comparison to the field region, convergence is reached at earlier times. This is due to the fact that structure growth is accelerated in the proto-cluster region so that its star formation becomes dominated by relatively massive galaxies much earlier than in the field.

Thus, the ‘S2’ and ‘S1’ simulations converge to ‘S3’ already at $`z_{conv}9`$ and $`z_{conv}8`$, respectively. Comparing the ‘S3’ and the ‘S4’ simulations we conclude that the ‘S3’, which has a mass resolution of $`M_p=2.4\times 10^8h^1`$ M, somewhat lower than ‘M3’, already accounts for all significant star formation. A simulation of the field region with mass resolution comparable to ‘S4’ is not available. We therefore estimate the redshift after which the ‘M3’ run accounts for all star formation by deriving the contribution to the total star formation from objects with masses in different ranges.

In Fig. 2 we plot $`\dot{M}_{}`$ as a function of the halo mass, $`M`$ at different redshift. The bins are 0.5 wide, in units of log $`M`$ and they have been chosen so that the lower limit is equivalent to our mass resolution. Before redshift $`10`$ the main contribution comes from the smallest mass objects, but larger and larger objects become important as the redshift decreases. This can be seen more clearly in Fig. 3 where the redshift evolution of the total star formation rate is plotted (solid line) together with the contribution from haloes with masses in different ranges. As expected, at the higher redshifts the major contribution comes from objects with mass of few $`10^9`$ M. At $`z11`$ their contribution drops below half of the total. We have plotted the contribution of haloes with mass up to $`10^{11}`$ M, as more massive objects are important only at $`z<5`$, and thus are not relevant for our study. From this analysis we can conclude that the ‘M3’ simulation accounts for most star formation at $`z\stackrel{<}{}11`$.

## 3 Radiative Transfer

To follow the propagation of ionizing radiation produced by the sources through the given IGM density distribution, we use the Monte Carlo (MC) radiative transfer code CRASH (Cosmological RAdiative transfer Scheme for Hydrodynamics) described in CFMR. It should be noted that in the present study a wider range of densities is encountered, reaching values of the optical depth of a few hundreds in the highest density regions. For this reason, additional tests, that we do not report here, have been run to check its reliability, finding that it reaches the same accuracy providing a right number of photon packets is emitted (see below). For clarity, we briefly summarize the main features of the numerical scheme relevant to the present study.

In the application of a MC scheme to radiative transfer problems, the radiation intensity is discretized into a representative number of monochromatic photon packets. Then, all the processes involved (e.g. packet emission and absorption) are treated statistically by randomly sampling the appropriate distribution function and solving the time-dependent ionization equation. Differently from CFMR, here we deal with more than one ionizing source. In this case, a better solution of the discretized time-dependent ionization equation is obtained if the time step $`\mathrm{\Delta }t`$ used in their eq. 10 is not treated statistically. Thus, here $`\mathrm{\Delta }t`$ is the time elapsed since a photon packet has gone through the cell for which the equation is being solved. The only inputs needed for the calculation are the gas density field and ionization state, as well as the source positions and emission properties. Then, CRASH provides the final gas ionization state, while the gas density field remains constant throughout the radiative transfer simulation. The above procedure is repeated for each output of the N-body simulation described in the previous Section. The inputs and the CRASH running time, $`t_{rt}`$, are changed accordingly; in particular, $`t_{rt}=t_{out}`$, the time between two outputs of the simulation. In the following, we discuss how the above inputs are derived and updated.

For each output of the N-body simulation, the dark matter density distribution, $`n_{dm}`$, is tabulated on a mesh using a Triangular Shaped Cloud (TSC) interpolation (Hockney & Eastwood 1981). We discuss the effect of neglecting density variations on sub-grid scales in Section 5. Assuming that the gas distribution follows that of the dark matter, we set the gas density in each grid cell by requiring the adimensional baryon density to be $`\mathrm{\Omega }_b=0.04`$. A number $`N_c=128^3`$ of cells have been used (see Section 5 for a discussion on this choice). As long as the sources are of stellar type, the presence of helium does not sensibly affect the H reionization process. For this reason, we consider a gas of pure hydrogen. Moreover, a simulation with an H/He gas, although feasible with the updated version of the code CRASH (Maselli, Ferrara & Ciardi 2003), would be extremely time consuming. The other quantities, such as temperature, $`T`$, and ionization fraction, $`x`$, are initialized as in CFMR and updated as described in the following paragraphs.

The source positions are obtained directly from the outputs of the N-body simulations, as described in the previous Section. In order to reduce the computational cost of the calculation, we group all sources in a grid cell into a single source placed at the cell center. Mass and luminosity conservation are assured. To further reduce the computational cost of the radiative transfer simulation, each source is modeled with a number of photon packets proportional to its luminosity, i.e. its star formation rate. In practice, this means that we randomly sample the source luminosity distribution function. For the range of densities considered here, we find that a number of photon packets equal to $`𝒩_p=\left[5\times 10^7\left(E/10^{60}\mathrm{erg}\right)\right]`$, where $`E`$ is the total ionizing energy emitted at each output of the simulation, is needed to reach numerical convergence in the volume averaged ionization fraction to within few percent. If we were rather interested, e.g., in resolving the ionization fronts, a much higher number of photons packets would have been required (Maselli, Ferrara & Ciardi 2003). Of the emitted packets, only $`f_{esc}𝒩_p`$ will actually be able to escape from the galaxies and break into the IGM. The packets have the same energy, $`_p=E/𝒩_p`$, but contain a different number of monochromatic photons. We determine $`E`$ as follows. As ionizing photon production decreases rapidly as stars age, the bulk of the ionizing radiation is produced by newly formed stars. We thus assume that a mass of stars $`M_{}=t_{out}\times _{i=1}^{n_s}\dot{M}_{,i}`$, is responsible for the all emission, neglecting the contribution from stars formed in previous outputs. This is a good approximation as long as $`t_{out}\stackrel{>}{}10^7`$ yr, the mean lifetime of an OB star. $`\dot{M}_{,i}`$ is the star formation rate of the i-th galaxy and $`n_s`$ is the number of galaxies at the output under consideration. We derive $`E`$ integrating over $`t_{out}`$ the time-dependent SED of a Simple Stellar Population of metal-free stars, with a total mass $`M_{}`$, distributed according to a Salpeter IMF (see Section 2).

Given the above density distribution, source positions and emission properties, we run the code CRASH for the time $`t_{out,k}`$, corresponding to the output, $`k`$, under consideration. Before running CRASH for the next output, $`k+1`$, we update the relevant physical quantities as it follows. Let’s assume that the i-th cell has an ionization fraction $`x_i`$; we thus assign to each particle in the cell, $`p_i`$, the value $`x_i`$. We follow the same procedure for all the cells. At the output $`k+1`$ the $`p_i`$ particles may have moved into another cell. We thus reconstruct the new distribution of particles, together with their ionization fraction. Then, when we grid the box of the simulation at the $`k+1`$ output, we assign to each cell an ionization fraction averaged over those of the particles that end up in that cell. Particles flowing in from the external boundaries inherit the properties of the particles in the cell. This procedure is performed for all the quantities relevant to the calculation, with the exception of the gas density, which is updated from the simulations as described above. In this way, we can follow the reionization process in an evolving universe.

Finally, to mimic the presence of an unknown background, photons escaping from a side of the box re-enter from the opposite side. At the early times, when HII regions are still confined around their sources, few photons exit the box. Their number increases rapidly as the regions start to overlap. This is consistent with the expected build up of an ionizing background (e.g. Gnedin 2000).

## 4 Results

As discussed in the previous Sections, to study the effect of the environment on the reionization process, we have performed the same radiative transfer simulation on a 20 $`h^1`$ Mpc comoving box positioned in a field region of the universe and on a 10 $`h^1`$ Mpc comoving box centered on a clustered region. It should be noticed that in the redshift range where the radiative transfer simulations are applied, the cluster has not yet collapsed. The application of the CRASH code in its first version, that does not follow the detailed temperature evolution, is thus appropriate, since less than 1% of the gas is in regions with $`T>>10^4`$ K (see also Mo & White 2002). In the following we present the results of our simulations.

In Fig. 4 we show illustrative slices cut through the simulation boxes. The six panels show the neutral hydrogen number density for the ‘M3’ field region (upper panels) and the ‘S3’ proto-cluster (lower panels), at redshifts, from left to right, $`z=16.5,12`$ and 8.5. The ionized regions (black) clearly grow differently in the two environments. In the field, high density peaks are uncommon and HII regions easily break into the IGM; in the proto-cluster the density is higher, so ionization is more difficult and recombination is faster. Moreover, many photons are initially needed to ionize the high density gas surrounding the sources, before photons can escape into the low density IGM. As a result, more photons are required to ionize the proto-cluster region; although in the proto-cluster the ionizing photon production per unit mass is higher at high redshift (see Fig. 1), filaments of neutral gas are still present after the field region is almost completely ionized.

This is more clearly seen in Fig. 5, where the volume averaged ionization fraction, $`x_v`$, is shown as a function of redshift for ‘M3’ and ‘S3’. This ionization fraction is defined as $`x_v=_ix_iV_i/V`$, where $`x_i`$ and $`V_i`$ are the ionization fraction and volume of the $`i`$-th cell respectively, and $`V`$ is the total volume of the box. The evolution proceeds differently in the two regions: initially the ionization fraction are comparable, but the proto-cluster gas ionizes more slowly, for the reasons discussed above. While the ‘M3’ field region already has $`x_v0.999`$ at $`z=8.3`$, $`x_v0.8`$ for the proto-cluster at the same redshift.

A comparison between volume and mass averaged ionization fractions is useful to understand how the reionization process proceeds. We define the latter quantity as $`x_m=_ix_iM_i/M`$, where $`M_i`$ is the mass in the $`i`$-th cell and $`M`$ is the total mass in the box. In Fig. 6 the redshift evolution of $`x_v`$ (filled circles) and $`x_m`$ (open circles) is compared both for the ‘M3’ field region (upper panel) and for the ‘S3’ proto-cluster (lower panel). The most interesting feature of these plots is that, while for the field region $`x_m`$ is always comparable with $`x_v`$, $`x_m>x_v`$ initially in the proto-cluster, i.e. at first relatively dense cells are ionized. This because the density immediately surrounding the sources must be ionized, before the photons can break into the IGM. This indicates that the assumption adopted by some authors (e.g. Miralda-Escudé, Haehnelt & Rees 2000), that lower density gas always gets ionized before higher density regions, fails in a cluster environment. The trend reverses at later times, when large, low density regions get ionized and the highest density peaks remain neutral or recombine. Although almost 90% of the proto-cluster volume is ionized by $`z8`$, this is true for only $`70\%`$ of the mass.

Some approaches to the study of the reionization process adopt the one-photon-per-baryon approximation. This is valid only in cases when recombination can be safely neglected. An estimate of the photon budget required to ionize the IGM in a more realistic case is shown in Fig. 7, where the ratio between the cumulative number of escaping ionizing photons, $`N_{phot}`$, and the number of baryons, $`N_b`$, is plotted as a function of redshift for the ‘M3’ field region (solid line) and the ‘S3’ proto-cluster (dotted line). High density peaks, where star formation preferentially takes place, are more abundant in the proto-cluster, so photon production is higher there than in the field region, at least up to $`z11`$. At $`z8`$, when $`x_v=0.999`$, about 5 photons per baryon have been produced in the field region. However, the $`15`$ photons per baryon produced in the proto-cluster are not enough to reionize it completely for the reasons discussed above.

## 5 Discussion

In the following we will discuss the convergence of our simulations with respect to mass and grid resolution. To assess the impact of mass resolution on our results, we have run a reionization simulation on the ‘M2’ field region. The redshift evolution of $`x_v`$ is shown in Fig. 8 as open circles. The ‘M2’ and ‘M3’ curves converge only at $`z<7`$. The difference observed at higher redshift is due to the mass resolution, which, for the ‘M2’ simulation, is not enough to resolve the objects that give the main contribution to the ionizing radiation, as is evident from the upper panel of Fig. 1. From Figs. 2-3 we are confident that the ‘M3’ simulation has the appropriate resolution to account for the bulk of ionizing photon production, at least after $`z11`$. It should be noticed that the formation and evolution of objects with $`M<10^9`$ M may be affected by feedback effects, that inhibit their star formation (e.g. Haiman, Rees & Loeb 1997; Omukai & Nishi 1999; Mac Low & Ferrara 1999; Barkana & Loeb 1999; Ciardi et al. 2000; Glover & Brand 2000; Nishi & Tashiro 2000; Kitayama et al. 2001; Glover & Brand 2002; Benson et al. 2002). Thus, these smaller mass structures may well play a minor role in the reionization process.

To assess the impact of grid resolution we have rerun our radiative transfer simulations with $`N_c=64^3`$. The resulting evolution of $`x_v`$ is shown in Fig. 8 by asterisks. From this figure it is evident that the grid resolution does not greatly affect the final results; the value of $`x_v`$ never differs in the two runs by more than 10 or $`20\%`$, although the case with $`N_c=64^3`$ produces slightly higher values. $`N_c>128^3`$ is not allowed by the mass resolution of our N-body simulation, as sampling fluctuations in the density field then become large.

As mentioned above, small mass objects, although negligible for ionizing photon production, are important because small-scale clumping can enhance recombination. Simulations which do not include such objects may underestimate the actual number of ionizing photons needed to reionize the IGM. To mimic clumping on unresolved scales, a so called clumping factor, $`C`$, is often introduced. It is defined as $`C=<n^2>/<n>^2`$, where $`n`$ is the density. The TSC technique used to derive our density distributions will lead us to underestimate the clumping factor. The TSC interpolation has the advantage of producing smoothed density values at the positions of the grid cells, as required for the radiative transfer computation. However, this technique is not adaptive, like an SPH kernel for example, and density variations on sub-grid scales are smoothed out, reducing the contrast of the DM density peaks where our sources are positioned. This results in an underestimate of the recombination rate and thus faster reionization. In order to address the importance of this effect we have measured densities at the particle positions using both the TSC interpolation and a 32-particle SPH smoothing kernel. The result of this comparison for ‘M3’ at redshifts $`z=10`$ and $`z=6`$ is shown in Fig. 9. The overall agreement for low and intermediate densities is reasonable, but the smoothing of the TSC scheme is clearly visible at the high densities. This is especially true at $`z=6`$. An estimate of the error introduced can be quantified by the ratio of the total recombination rate derived with the two different schemes, $`r_{SPH}/r_{TSC}`$. This ratio is equal to 6.24 (20.87) at $`z=10`$ (6), suggesting that we are substantially underestimating the recombination rate in the high density regions surrounding the sources. We are interested primarily in the higher redshifts, where agreement is somewhat better. Since almost all the high density regions where the above corrections are important host sources, the error in the recombination rate can be considered part of the uncertainty on our adopted escape fraction (see below). Moreover, the effect of clumping is alleviated by photoevaporation, which tends to weaken the effect of gas on unresolved scales as ionization precedes (Shapiro et al. 2003).

The number of available ionizing photons depends on the choice of the stellar spectrum, IMF and escape fraction. As the first stars form out of gas of primordial composition, they are thought to be very massive, resulting in a higher ionizing photon emission, when compared with higher metallicity stars (e.g. Larson 1998). For example, Schneider et al. (2002) propose that enrichment is needed to permit further fragmentation into stars with significantly lower masses, inducing a transition from a top-heavy to a more conventional IMF. The debate on the metallicity and the IMF of the first stars is still very lively and many authors argue that metal-free stars are distributed according to a top-heavy IMF. On the other hand, Christlieb et al. (2002) recently discovered a low mass star with metallicity only $`3\times 10^6`$ that of the Sun. To be conservative, in this paper we have assumed metal-free stars with a Salpeter IMF.

Even more uncertain is the value of the escape fraction. Values derived in previous work vary from 0 to 60% (see e.g. Wood & Loeb 1999; Ricotti & Shull 2000; Ciardi, Bianchi & Ferrara 2002 for theoretical work and e.g. Leitherer et al. 1995; Hurwitz et al. 1997; Bland-Hawthorn & Maloney 1999; Steidel, Pettini & Adelberger 2001 for observations), according to the galaxy mass, redshift and density distribution. Given these uncertainties and the impossibility of accurately modeling this parameter for every galaxy (see e.g. Ciardi, Bianchi & Ferrara 2002 for a discussion), we regard $`f_{esc}=5\%`$ as a constant, reference value. It should be noticed that, as the resolution of the radiative transfer calculation is equivalent to the dimension of a single grid cell (which is always larger than the dimension of the sources), $`f_{esc}`$ is actually an “effective” escape fraction, i.e. it includes both the contribution to the absorption of the interstellar medium in the emitting object (“real” escape fraction) and the one of the IGM in the same cell as the object.

Although it is not possible to disentangle different choices for the source emission properties from clumping on small scales, their expected global effect can be estimated introducing the parameter $`𝒜=C_n^1(f_{esc}/0.05)(N_\gamma /1\times 10^{61}\mathrm{M}_{}^1)`$, where $`C_n`$ is the clumping factor normalized to the value of our simulations (e.g. $`C_n=2`$ means that the clumping factor is twice that of the simulations) and $`N_\gamma `$ is the total number of ionizing photons per unit solar mass of formed stars produced by the adopted spectrum and IMF. $`N_\gamma `$ is normalized to a Salpeter IMF and a spectrum typical of metal-free stars (CFMR). From a “back of the envelope” calculation, we conclude that if $`𝒜=10`$, the $`5`$ photons per baryon required to reionize the IGM are already available at $`z13`$, while if $`𝒜=0.1`$ the above condition is never reached. As the most recent observations of high redshift quasars (Becker et al. 2001; Djorgovski et al. 2001; Pentericci et al. 2002) suggest that the universe was already ionized at $`z6.3`$, the lowest value of $`𝒜`$ consistent with this constraint is $`0.7`$. A lower value of $`𝒜`$ would require additional sources of ionizing photons to produce the observed ionization level. In a forthcoming paper we will discuss the possible contribution to ionizing radiation from a primordial population of quasars (see also Wyithe & Loeb 2002). In a recent paper Ricotti (2002) has suggested that a possible contribution to ionizing photons could also come from globular clusters. Given our conservative assumptions on the value of the escape fraction and the IMF, a higher value of the clumping factor would also be compatible with a reionization driven by primordial stellar sources.

Recent results from the WMAP satellite require a mean optical depth to Thomson scattering $`\tau _e0.17`$, suggesting that reionization must have begun at relatively high redshift (e.g. Kogut et al. 2003; Spergel et al. 2003). Reconciling an early reionization with observations of the Gunn-Peterson effect in $`z>6`$ quasars, which imply a volume-averaged neutral fraction above $`10^3`$ at $`z6`$, may be challenging. In Ciardi, Ferrara & White (2003) we run additional simulations and discuss our model in view of these new data.

## 6 Conclusions

We have studied the IGM reionization process in a full cosmological context with structure evolution and a realistic galaxy population. We have used a combination of high-resolution N-body simulations (to describe the dark component of matter), a semi-analytic model of galaxy formation (to track the gas evolution) and the Monte Carlo radiative transfer code CRASH (to follow the propagation of ionizing photons into the IGM; CFMR). The process has been followed in the largest volume ever used for this kind of study, with a comoving length of $`L20h^1`$ Mpc, a mass resolution of $`10^8`$ M and tree from periodic boundary conditions. This allows us to resolve the primary objects responsible for the production of ionizing radiation and, at the same time, to simulate a region of the universe with “mean” properties. To study the effect of environment on the reionization process, we have performed the same radiative transfer simulation on a $`10h^1`$ Mpc comoving box, centered on a proto-cluster region. The main results discussed in this paper can be summarized as follows.

(1) Our field region of the universe gets reionized by the simulated galaxy population at a redshift $`z_{ion}8`$, while the proto-cluster, although it produces a higher number of ionizing photons, gets ionized later. This is due to the fact that high density gas, which is more common in the proto-cluster, is more difficult to ionize and recombines much faster.

(2) While for the field region $`5`$ photons per baryon are enough to produce $`x_v=0.999`$ at $`z8`$, the $`15`$ photons per baryon produced by this time in the proto-cluster are not enough to reionize it.

(3) For the same reason, while for the field region the mass and volume averaged ionization fractions, $`x_m`$ and $`x_v`$ respectively, are always comparable, $`x_m`$ substantially exceeds $`x_v`$ for the proto-cluster, at high redshifts. This is because the high density regions surrounding the sources have to be ionized first, before the photons can break out into the low density IGM.

(4) The primordial stellar sources considered in this study give a value of the reionization epoch consistent with observation, without invoking the presence of additional sources of ionization. Our conservative assumptions on the value of the escape fraction and the IMF can accommodate the additional clumping which remains unresolved in our scheme.

(5) The mass resolution of the simulations can deeply affect the final results. Objects with total masses of $`M10^9`$ M or less must be resolved to account for the bulk of the star formation.

## Acknowledgments

We would like to thank V. Springel for providing us with the cluster simulations and an anonymous referee for his/her useful comments. B.C. is grateful to F.S. for his patience with computer related questions and A. Ferrara for useful comments. This work has been partially supported by the Research and Training Network “The Physics of the Intergalactic Medium” set up by the European Community under the contract HPRN-CT-2000-00126.


format: text

text:
# Magnetic-Field-Induced Localization of Quasiparticles in Underdoped La2-xSrxCuO4 Single Crystals

## Abstract

Magnetic-field-induced ordering of electrons around vortices is a striking phenomenon recently found in high-$`T_c`$ cuprates. To identify its consequence in the quasiparticle dynamics, the magnetic-field ($`H`$) dependence of the low-temperature thermal conductivity $`\kappa `$ of La2-xSrxCuO4 crystals is studied for a wide doping range. It is found that the behavior of $`\kappa (H)`$ in the sub-Kelvin region changes drastically across optimum doping, and the data for underdoped samples are indicative of unusual magnetic-field-induced localization of quasiparticles; this localization phenomenon is probably responsible for the unusual “insulating normal state” under high magnetic fields.

Competing order is an important emerging concept in high-$`T_c`$ superconducting cuprates Sachdev . Recently, Hoffman et al. Hoffman  reported scanning tunneling microscopy (STM) studies of Bi2Sr2CaCu2O8+δ (Bi2212) single crystals in magnetic fields, which revealed a “checkerboard-like” order of quasiparticle (QP) states with four-unit-cell periodicity surrounding vortex cores. This result documented a magnetic-field-induced charge density wave (CDW), which is in good correspondence with the magnetic-field-induced spin density wave (SDW) found by neutron scattering in La2-xSrxCuO4 (LSCO) Lake1 ; Lake2  and La2CuO4+δ Khaykovich , and these phenomena are likely to be results of competing antiferromagnetic and superconducting orders in those materials. In fact, several theoretical models Demler ; Franz ; Zhu ; Chen ; Kivelson  have been proposed to describe the coexistence of $`d`$-wave superconductivity and spin/charge order in the vortex state as a consequence of competing orders. However, it is hardly known how such competition and magnetic-field-induced order affect the QP dynamics, the details of which would allow us to understand the nature of the novel magnetic-field induced states.

In this Letter, we report that the magnetic-field dependence of thermal conductivity, $`\kappa (H)`$, of LSCO at low temperature demonstrates that the magnetic-field-induced order leads to unusual localization of QPs; moreover, we found that there is a distinct change in the behavior of $`\kappa (H)`$ near optimum doping and the QP localization occurs only in the underdoped samples. This magnetic-field-induced localization of QPs in the underdoped regime is clearly in correspondence with the well-known metal-to-insulator crossover, found in the low-temperature normal state under 60-T magnetic field, that occurs at optimum doping Boebinger . Moreover, this result suggests that the unusual “insulating normal state” of the cuprates under 60 T, that is characterized by a log(1/$`T`$) divergence in resistivity Ando1 , is caused by the magnetic-field-induced order.

In high-$`T_c`$ cuprates, it is established that the superconducting gap $`\mathrm{\Delta }`$ has essentially the $`d_{x^2y^2}`$ symmetry, which has four nodes (where the gap magnitude vanishes) along the diagonals of the square Brillouin zone. Because of these nodes, QPs are easily created both by thermal fluctuations and by the impurity scattering, and the thermal conductivity $`\kappa `$ is a useful bulk probe of the QP excitations and their transport behavior Durst ; Kubert . In the mixed state, where the magnetic field enters the superconductor in the form of quantized vortices, it is known that the magnetic field induces extended zero-energy QP states near the gap nodes and their population increases as $`H^{1/2}`$ Durst ; Kubert ; Moler ; this so-called “Volovik effect” is due to the Doppler shift of the QP energy spectrum in the presence of a supercurrent flowing around each vortex Volovik . As a result, the QP heat transport is enhanced in magnetic fields at low temperatures (usually in the sub-Kelvin region) Chiao ; Aubin ; Ando2 ; Proust , while at higher temperatures the QP heat transport is suppressed in magnetic fields because of the dominance of vortex scattering of QPs Krishana ; Ong ; Ando3 . Interestingly, at intermediate temperatures a “plateau” shows up in the $`\kappa (H)`$ profile of Bi2212 Krishana , which, after a long debate (see Ref. Ando3  and references therein), can be understood to be a result of the competition between the increase in the QP population and the decrease in their mean free path with $`H`$ Ando2 . Till now, the magnetic-field-induced enhancement of QP transport was studied in optimally-doped YBa2Cu3O7-δ (YBCO) Chiao  and Bi2212 Aubin ; Ando2  and in overdoped Tl2Ba2CuO6+δ Proust , but its doping dependence has been scarcely known.

Here we choose to study the low-temperature heat transport of LSCO, in which the hole doping can be well controlled over a wide range. High-quality La2-xSrxCuO4 single crystals ($`x`$ = 0.08, 0.10, 0.14, 0.17 and 0.22) are grown by the traveling-solvent floating-zone technique Ando4 . The crystals are cut into rectangular platelets with typical size of $`1.5\times 0.5\times 0.1`$ mm3, where the $`c`$ axis is perpendicular to the platelets within an accuracy of 1, determined by the X-ray Laue analysis. The samples at $`x`$ = 0.08–0.17 are annealed at 800–880 C in air for 60 hours, followed by rapid quenching to room temperature, to remove oxygen defects, while the sample at $`x`$ = 0.22 is annealed in oxygen and quenched to 77 K. Magnetic susceptibility measurements show that the superconducting transition temperature $`T_c`$ is 22, 28.5, 35.5, 36.5 and 28 K for $`x`$ = 0.08, 0.10, 0.14, 0.17 and 0.22, respectively. The zero-field heat transport in these samples at milli-Kelvin temperatures has already been reported in a previous paper Takeya . In this work, both the temperature and the magnetic-field dependences of $`\kappa `$ are measured from 0.3 to 7 K in a 3He refrigerator by using the conventional steady-state “one heater, two thermometer” technique Ando2 . The magnetic field up to 16 T is applied along the $`c`$ axis of the crystals while the heat current flows in the $`ab`$ plane. To avoid complications that are associated with the vortex-pinning-related hysteresis Aubin , the $`\kappa (H)`$ data are taken in the field-cooled procedure Ando2 ; Ando3 .

Figure 1 shows the zero-field $`\kappa (T)`$ of LSCO crystals in a wide temperature range, where the additional data from 5 to 150 K are taken using a Chromel-Constantan thermocouple in a 4He cryostat, and the very low temperature $`\kappa (T)`$ data Takeya  are also shown in the inset. The measured thermal conductivity is a sum of the electron (or QP) term $`\kappa _e`$ and the phonon term $`\kappa _{ph}`$, whose temperature dependences are proportional to $`T`$ and $`T^3`$ in the low-$`T`$ limit, respectively. Thus, in the $`\kappa `$/$`T`$ vs $`T^2`$ plot (Fig. 1 inset), it is clear that the $`\kappa _e`$ component increases with $`x`$. In the main panel of Fig. 1, the superconducting transition can barely be recognized by a weak hump in $`\kappa `$, which is easiest to see in the optimally-doped sample ($`x`$ = 0.17).

Figure 2 shows the magnetic-field dependences of $`\kappa `$ for optimally-doped ($`x`$ = 0.17) and overdoped ($`x`$ = 0.22) LSCO. The optimally-doped sample shows a $`\kappa (H)`$ behavior very similar to that in Bi2212 Aubin ; Ando2  and YBCO Chiao ; Ong ; namely, $`\kappa `$ is enhanced with $`H`$ at very low temperatures, but is suppressed with $`H`$ at higher temperatures. Not surprisingly, a clear “plateau” shows up in the $`\kappa (H)`$ profile at 7 K. It is useful to note that the low-$`T`$ $`\kappa (H)`$ data of the optimally-doped sample can roughly be described by the formula $`\kappa (H)/\kappa (0)=1+aH^{1/2}`$ ($`a`$ is a temperature-dependent parameter), as is shown in Fig. 2(a). This $`H`$-dependence is expected Aubin  if the vortex scattering is negligibly small and the rise in $`\kappa `$ is attributed simply to the increase in QP population due to the Volovik effect. (Note that the phononic contribution to $`\kappa `$ is believed to be independent of $`H`$ in the cuprates Chiao ; Aubin ; Ando2 ; Proust ; Krishana ; Ong ; Ando3 .) The reasonable fits in Fig. 2(a) demonstrate that the vortex scattering becomes relatively insignificant at low temperatures at optimum doping. In overdoped LSCO ($`x`$ = 0.22), the general trend of $`\kappa (H)`$ is similar, although there is some difference at high fields; as is shown in Fig. 2(b), the behavior of $`\kappa (H)`$ follows $`H^{1/2}`$ law at low fields but noticeably deviates from $`H^{1/2}`$ at high fields, presenting a “plateau-like” feature. This flattening of $`\kappa (H)`$ probably indicates an increasing importance of QP scattering at high fields. One possible cause of this feature is the effect of inelastic electron-electron scattering, which may remain important at very low temperatures in overdoped samples Mackenzie . In any case, the magnetic-field dependence of the QP heat transport in optimum and overdoped LSCO is qualitatively similar to that in the optimally-doped Bi2212 and YBCO.

The most surprising result of the present work is the behavior of $`\kappa (H)`$ in underdoped LSCO ($`x`$ = 0.08, 0.10 and 0.14). As shown in Fig. 3, $`\kappa `$ always decreases quickly with $`H`$ at all temperatures down to 0.36 K in all the underdoped crystals studied, which is vastly different from the $`\kappa (H)`$ behavior of optimum and overdoped crystals. To understand the implication of this result, it should first be recognized that the impurity scattering cannot be the source of such a strong magnetic-field suppression of QP heat transport, because it has been demonstrated Ando3  that impurities just strongly diminishes the $`H`$ dependence of $`\kappa `$. Therefore, the effect of magnetic field on the QP heat transport must have been changed fundamentally in underdoped samples. When the data in Fig. 3 are compared to those in Fig. 2, it is obvious that both the magnetic-field enhancement of $`\kappa `$ at low $`T`$ and the “plateau” feature at high $`T`$ are absent in underdoped samples, which points to an increased role of vortex scattering and a diminished role of magnetic-field-induced QPs. Given that the specific heat measurements have already demonstrated SJChen  that the QP density-of-states are always enhanced with $`H`$ (i.e., QPs are certainly created by magnetic fields) in underdoped LSCO, it is most reasonable to conclude the following: (i) the magnetic-field-induced QPs are localized and contribute little to the heat transport in underdoped LSCO, and (ii) vortices strongly scatter QPs in underdoped samples even at low $`T`$, while they do not effectively scatter QPs in optimum and overdoped samples at low $`T`$.

These results are most naturally understood in the light of the recently-found magnetic-field-induced CDW/SDW in the cuprates Hoffman ; Lake1 ; Lake2 ; Khaykovich . For LSCO, neutron scattering experiments have shown Lake2 ; Khaykovich  that the magnetic-field-induced SDW is dynamical at optimum doping, while it becomes a static, zero-energy object in underdoped samples. Thus, the magnetic-field-induced order is expected to be relevant to the low-temperature properties only in underdoped LSCO. Moreover, the STM study of Bi-2212 found checkerboard-like CDW around vortices Hoffman , which is naturally expected to cause enhanced QP scattering. Therefore, it is likely that in underdoped LSCO the magnetic-field-induced order is responsible for both the QP localization and the enhanced QP scattering off the vortices. It is thus reasonable to conclude that the peculiar $`\kappa (H)`$ behavior of underdoped LSCO points to novel magnetic-field-induced localization of QPs that is caused by the magnetic-field-induced order around vortices.

Figure 4 shows the hole-doping dependence of $`\kappa `$ at 0.36 K in 0 and 16 T, which summarizes our main finding. In 0 T, $`\kappa `$ shows a relatively smooth change with $`x`$, but $`\kappa `$ in 16 T shows a pronounced jump at optimum doping. As one can see in Fig. 4, this jump is caused by the fact that $`\kappa `$ is suppressed with $`H`$ in the underdoped regime, while it is enhanced with $`H`$ in the optimum and overdoped regime; we discussed that the former is caused by the magnetic-field-induced localization of QPs that is likely associated with the competing order, while the latter reflects extended QPs created by magnetic fields. Therefore, Fig. 4 depicts the fact that the role of magnetic field changes drastically at optimum doping.

An intriguing implication of the present result is that the magnetic-field-induced QP localization may also be responsible for the “insulating” normal-state resistivity under 60 T Boebinger ; Ando1 . The fact that the “insulating normal state” under 60 T is observed only in underdoped samples Boebinger  is strongly suggestive of the correspondence. In fact, it is easy to imagine that the charge-ordered region around vortices completely overlap in 60-T field and all the electrons are incorporated into the magnetic-field-induced ordered state. Thus, it is tempting to assert that the log(1/$`T`$) behavior of the resistivity under 60 T Ando1  is a property of the state where the magnetic-field-induced order proliferates. In this regard, it is worth noting that a recent Nernst effect measurement found a sizable Nernst signal in the “insulating normal state” in high magnetic fields Capan . In fact, the behavior of the Nernst signal below $`T_c`$ is not much different between underdoped and overdoped LSCO, which probably reflects the fact that the observed Nernst signal primarily comes from the vortex motion and is insensitive to the QP contribution Capan ; Wang . Incidentally, the Nernst effect has been utilized to measure the upper critical field $`H_{c2}`$ of cuprates, where $`H_{c2}`$ was found to show a rather smooth change across optimum doping Wang ; this fact rules out the possibility that the observed change in $`\kappa (H)`$ between $`x`$ = 0.14 and 0.17 might be caused by a sudden change in $`H_{c2}`$, which would cause the reduced field $`H/H_{c2}`$ to vary drastically.

It is useful to note that our previous results in zero magnetic field have shown Takeya  that the residual electronic thermal conductivity divided by $`T`$, $`\kappa _{\mathrm{res}}/T`$, is much larger than what the normal-state resistivity under 60 T would suggest in the underdoped LSCO, which may indicate a strong violation of the Wiedemann-Franz law. The present new data show that in underdoped LSCO the QP transport is strongly suppressed with $`H`$ even at low $`T`$, and thus it is expected that the $`\kappa _{\mathrm{res}}/T`$ value in high magnetic field may smoothly match the value estimated from the normal-state resistivity under 60 T. This is an indication that the magnetic-field-induced order around vortices ultimately changes the global properties of underdoped LSCO in high magnetic field. Therefore, the thermal transport properties give evidence that the magnetic field can cause a drastic change in the low-energy physics of underdoped cuprates because of the magnetic-field-induced localization of QPs, which is likely related to the inherent competition between antiferromagnetic and superconducting orders.

In summary, we measure the magnetic-field dependences of the thermal conductivity of LSCO single crystals for a wide doping range. It is found that the $`\kappa (H)`$ behavior at low-$`T`$ gives evidence for novel magnetic-field-induced localization of QPs in the underdoped regime. We discuss that this phenomenon is associated with the magnetic-field-induced spin/charge order, and is probably responsible for the unusual “insulating normal state” under high magnetic fields.

###### Acknowledgements.

We thank S. A. Kivelson, A. N. Lavrov, and P. A. Lee for helpful discussions.


text:
# Mode Repulsion and Mode Coupling in Random Lasers

## Abstract

We studied experimentally and theoretically the interaction of lasing modes in random media. In a homogeneously broadened gain medium, cross gain saturation leads to spatial repulsion of lasing modes. In an inhomogeneously broadened gain medium, mode repulsion occurs in the spectral domain. Some lasing modes are coupled through photon hopping or electron absorption and reemission. Under pulsed pumping, weak coupling of two modes leads to synchronization of their lasing action. Strong coupling of two lasing modes results in anti-phased oscillations of their intensities.

A random laser is a special kind of laser whose feedback is due to disorder-induced scattering, as opposed to reflection of a conventional laser. Over the past few years, there have been many studies on random lasers with coherent feedback.  One important aspect of a random laser is the interaction of the lasing modes, which has not been well understood.  As previously pointed out , adding optical gain to a random medium creates a new path to study wave transport and localization. Because the modes of coherent random lasers are eigenstates of the random systems, the interaction of the lasing modes reflects the interaction among the eigenstates. So far most studies on eigenstates of random systems are focused on statistical properties such as the density of states and the distribution of spectral spacing between neighboring eigenstates.  However, there has been no detailed study on the interaction between two eigenstates in a disordered system, e.g., the spatial overlap of two eigenstates and the photon hoping between them. Such $`microscopic`$ information is crucial to understanding the behavior of light in the random system. By studying the interaction of random lasing modes, we will gain some knowledge of the interaction among the eigenstates of a disordered system. Furthermore, this study will provide an understanding of the escape channels for photons in an eigenstate, that is important to the localization theory.

In this letter, we present experimental and theoretical studies on the interaction of modes in random lasers with coherent feedback. We find the interaction is rather complicated. On one hand, most lasing modes repel each other either spatially or spectrally. On the other hand, some lasing modes are coupled. Energy exchange between the coupled modes leads to cooperated lasing phenomenon.

The spatial repulsion of lasing modes was first observed in the numerical simulation of random lasers.  It results from gain competition and spatial localization of the modes in an active random medium. The gain spectrum is assumed to be homogeneously broadened. In a finite-sized random medium, the mode repulsion leads to a saturation in the number of lasing modes.  Experimentally we observed the number of lasing modes in micropatterned ZnO polycrystalline films saturates at high pumping intensity.  The saturated number of lasing modes increases linearly with the sample volume.

In random media with inhomogeneously broadened gain spectra, we observed a different kind of mode repulsion. The lasing modes repel each other in the spectral domain instead of the spatial domain. The disordered samples used in our experiment are poly(methyl methacrylate) (PMMA) sheets containing rhodamine 640 perchlorate dye and titanium dioxide (TiO2) microparticles. The concentration of rhodamine 640 is about 10 mM. Each individual dye molecule’s transition frequency is shifted by differing local strains and distortions in the immediate environment of the glassy host. This gives rise to a large inhomogeneous broadening of the dye transition, which is confirmed by good fitting of the absorption band tail with a Gaussian profile (the error of Gaussian curve fitting is ten times smaller than that of Lorentzian curve fitting). The TiO2 particles have an average diameter of 0.4 $`\mu `$m. The density of the microparticles is $``$ 1.4 $`\times `$ 1012 cm-3. The transport mean free path, measured in the coherent backscattering experiment, is 2.1 $`\mu `$m.

The rhodamine 640 dye molecules are optically excited by the second harmonics of a pulsed Nd:YAG laser (10 Hz repetition rate, 25 ps pulse width). The pump spot is about 50 $`\mu `$m in diameter at the sample surface. The spectrum of emission from the sample is measured by a 0.5-meter spectrometer with a liquid nitrogen cooled CCD array detector. The spectral resolution is 0.06 nm.

Below the lasing threshold, the emission spectrum is very broad. When the pumping intensity exceeds the threshold, discrete narrow peaks emerge in the emission spectrum. The number of lasing modes increases with the pumping intensity. Eventually the spectral density of the lasing peaks becomes saturated. Namely, the wavelength spacing between adjacent lasing modes no longer decreases as the pumping intensity increases. New lasing peaks appear not in between the existing peaks but in the spectral region that is away from the gain maximum and has few lasing peaks. The lasing spectrum, shown in Fig. 1 (a), reveals that the lasing peaks are more or less regularly spaced in wavelength. We calculated the normalized spectral correlation function $`C(\mathrm{\Delta }\lambda )I(\lambda )I(\lambda +\mathrm{\Delta }\lambda )𝑑\lambda /[I(\lambda )]^2𝑑\lambda `$, where $`I(\lambda )`$ is the emission intensity at wavelength $`\lambda `$. Fig. 1 (b) shows the normalized spectral correlation function $`C(\mathrm{\Delta }\lambda )`$ obtained from the lasing spectrum in Fig. 1 (a). The periodic oscillation of $`C(\mathrm{\Delta }\lambda )`$ with $`\mathrm{\Delta }\lambda `$ confirms the lasing modes are evenly distributed in wavelength. The oscillation period gives the average separation between adjacent lasing modes, which is 0.8 nm. Note that the minimum of $`C(\mathrm{\Delta }\lambda )`$ does not reach zero because the lasing peaks are on top of a broad band of amplified spontaneous emission.

Recently we developed a spectrally resolved speckle technique which gives the spatial profile of individual lasing modes at the sample surface.  Using this technique, we find the lasing modes in Fig. 1 (a) have different sizes. The radius of the lasing modes varies from 100 $`\mu `$m to 150 $`\mu `$m, which is larger than the radius of the pump spot. To acquire gain, the lasing modes must overlap with the excitation area. This suggests most lasing modes spatially overlap with each other. Owing to a large inhomogeneous broadening of the dye molecules, lasing can occur simultaneously in modes that are spatially overlapped, as long as the frequency separation between the modes is larger than the homogeneous linewidth of the dye molecules. This is because a lasing mode interacts only with the dye molecules whose transition frequency is in resonance with the lasing frequency, and depletes the population inversion of these dye molecules. As a result, a hole is burned in the gain spectrum, and its width is about twice of the homogeneous linewidth of the dye molecules.  Nevertheless, the dye molecules, whose frequency is more than one homogeneous linewidth away from the frequency of the lasing mode, are essentially unaffected by this lasing mode. They provide optical gain for other modes in the same area. Thus the frequency spacing of the lasing modes is limited by the homogeneous linewidth of the dye molecules. At room temperature, the homogeneous linewidth of rhodamine 640 dye molecules in PMMA is $``$ 0.5 nm.  The linewidth of the hole burned by a lasing mode is about 1 nm, close to the wavelength separation between adjacent lasing modes. This confirms the regular spectral spacing of the lasing modes is caused by the hole burning of the inhomogeneously broadened gain spectrum. The spectral repulsion of lasing modes can only be observed when the excitation area is smaller than the spatial extent of lasing modes. When the pump area is much larger than the size of lasing modes, there is no minimum frequency separation between adjacent lasing modes. This is because the lasing modes are not necessarily overlapped in space. If they are spatially separated, the lasing modes can have the same frequency.

In addition to mode repulsion, we also observed mode coupling in random lasers. Despite most lasing modes in Fig. 1 (a) are well separated spectrally, we do see two or three lasing modes partially overlapped in the spectrum, e.g. the two lasing modes centered at 626.9 nm and 627.1 nm. Since their frequency separation is less than the homogeneous linewidth of the dye molecules, these two lasing modes interact with the same group of excited dye molecules. The gain competition would not allow them to lase simultaneously unless they are coupled, i.e., there is an energy exchange between the two modes.

Figure 2 shows the growth of two coupled lasing modes with the pumping intensity. Just above the lasing threshold, a relatively broad lasing peak appears in the emission spectrum. As the pumping intensity increases, this lasing peak is slightly blue-shifted and split into two narrow peaks. It seems initially the linewidth of the two lasing modes is larger than their spectral separation, thus they appear to be a single broad peak. As the optical gain increases, the lasing peaks become narrower. When their linewidth is reduced to below their spectral separation, these two lasing modes can be spectrally resolved, leading to a splitting of the initial peak. We measured the temporal evolution of emission from these two modes with a streak camera. Lasing in these two modes is synchronized, i.e., lasing starts and stops at the same time.

To understand the coupling of the random laser modes, we performed numerical simulations using the time-dependent model for random lasers.  The Maxwell’s equations are coupled to the rate equations of electronic population. Because the coupled modes are within one homogeneous line and interact with the same group of molecules, we can assume the electronic energy levels are homogeneously broadened. The homogeneous linewidth is 5 THz. Temporal evolution of the electromagnetic field is obtained with the finite difference time domain (FDTD) method. Fourier transform is used to obtain the lasing spectra. The one-dimensional random system is made of binary layers with dielectric constants $`ϵ_1`$ = 2.5 and $`ϵ_2`$ = 5.29. The layer thickness $`a=a_0(1+W_a\gamma )`$, $`b=b_0(1+W_b\gamma )`$, where $`a_0`$ = 48 nm, $`W_a`$ = 0.1, $`b_0`$ = 80 nm, $`W_b`$ = 0.6, $`\gamma `$ is a random value in the range \[-1,1\]. In the inset of Fig. 3, we plot the calculated lasing spectra for a system of 100 pairs of binary layers at two pumping rates. The pumping starts at $`t`$ = 0 and stops at $`t`$ = 15 ps. At a pumping rate $`P_r=2\times 10^6`$ s-1, the spectrum has one relatively broad peak. This peak is split into two narrow peaks at a pumping rate $`P_r=4\times 10^6`$ s-1. Figure 3 shows the temporal evolution of the amplitude of the electric field $`A_E`$ for the two lasing peaks. Laser pulses in these two modes rise at the same time. The pulse duration is nearly identical. These results agree with the experimental data.

The synchronization of the two lasing modes results from coupling. The quality factors of these two modes are different. Our previous study of random lasers demonstrates that lasing modes with different quality factors exhibit unsynchronized temporal behaviors.  This conclusion holds only for non-interacting modes. However, lasing modes in a random medium can be coupled through photons and electrons. Namely, photons can hop from one mode to another as long as these two modes have spatial and spectral overlap. In addition, photons in one mode can be absorbed by electrons in the medium, which subsequently remit photons into another mode. To be more specific, let us consider two coupled modes A and B. The quality factor of mode A is higher than that of mode B. After pumping starts, as optical gain increases, the lasing threshold for mode A is reached first. The electromagnetic field of mode A increases quickly. When there are more photons in mode A than in mode B, there is a net energy flow from mode A to mode B through photon hopping or electron absorption and reemission. This seeding increases the photon number of mode B and accelerates lasing process in mode B. Thus mode A provides additional pumping to mode B. If the quality factor of mode B is not much lower than that of mode A, mode B can reach the lasing threshold while mode A still builds up. That is why we observe the two modes lase almost simultaneously. Similarly, when the pumping is stopped, photon number of mode B decreases faster than that of mode A due to higher loss rate of mode B. The difference of photon numbers in these two modes leads to a net energy flow from mode A to mode B. Thus mode B provides an additional leak channel for mode A, and photons in mode A are drained faster. As a result, the two modes stop lasing at almost the same time.

Next we increased the quality factors of the coupled modes by enlarging the system to 120 pairs of binary layers. We also increased the pumping time to 25 ps to make the laser pulses longer. As shown in Fig. 4, the dynamics of the coupled modes is very different from that shown in Fig. 3. Laser light intensities in these two modes oscillate in time. Their oscillations are 180 out of phase. Namely, when one mode reaches the maximum intensity, the other mode is at the minimum intensity. These antiphased oscillations reveal energy exchange between the two coupled modes. From the numerical simulation, we find such oscillations occur only when the coupled modes have high quality factors and the pumping time is not too short. This is because photons must stay in the medium long enough so that they are transferred back and forth between the coupled modes before escaping the system. The oscillation period is inversely proportional to the energy exchange rate or the coupling constant. The coupling of the two modes depends on their spectral and spatial overlap. There are two coupling regimes. In the weak coupling regime, the coupling constant of the two modes is smaller than their decay rates. In the strong coupling regime, the coupling constant is larger than the decay rates of the two modes. The oscillations occur only in the strong coupling regime. In our current experiment with the three-dimensional (3D) random medium, the quality factors of the lasing modes are too low to reach the strong coupling regime. The antiphased oscillations of the coupled modes could not be observed.

To confirm the spectral repulsion of lasing modes, we simulated random laser with inhomogeneously broadened gain spectrum. The inhomogeneous linewidth is 150 THz, 30 times of the homogeneous linewidth. To increase the spatial overlap of the eigenstates, we reduce the randomness to $`W_a`$ = 0 and $`W_b`$ = 0.2 so that the modes are more extended spatially. The size of the system is increased to 700 pairs of binary layers ($``$ 100 $`\mu `$m). However, only the central part of the system with the length 30 $`\mu `$m is pumped. Since the lasing modes must overlap with the pumped region, they overlap with each other. At very high pumping rate, the lasing modes are almost regularly spaced in frequency. The spectral correlation function $`C(\mathrm{\Delta }\nu )`$ exhibits periodic oscillation with $`\mathrm{\Delta }\nu `$. The frequency separation of neighboring lasing modes is about 7 THz. It is much larger than the average frequency spacing of the eigenstates in the random system. This means most eigenstates do not lase even at very high pumping rate due to cross gain saturation. More specifically, there are about three eigenstates within one homogeneous line. Cross gain saturation allows only one of them to lase. Thus the frequency spacing of the lasing modes is determined by the homogeneous linewidth.

Finally we increased the pumping time in our simulation. All the phenomena presented above persist. Despite our experiments were performed with short pump pulses, our numerical simulation results confirm that mode repulsion and mode coupling exist at longer pumping time.

In conclusion, we studied the interaction of lasing modes in random media. In a homogeneously broadened gain medium, cross gain saturation leads to spatial repulsion of lasing modes. In an inhomogeneously broadened gain medium, mode repulsion occurs in the spectral domain. Some lasing modes are coupled through photon hopping or electron absorption and remission. Under pulsed pumping, weak coupling of two modes leads to synchronization of their lasing action. Strong coupling of two lasing modes results in anti-phased oscillations of their intensities.

Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. W07405-Eng-82. This work is supported partially by National Science Foundation under Grant No. DMR-0093949. H.C. acknowledge support from the David and Lucille Packard Foundation and the Alfred P. Sloan Foundation.


text:
# Termination of 4-fold canonical flips

(Date: 2003/1/15)

## Abstract.

We prove the termination of 4-fold canonical flips.

This paper is a supplement of \[KMM, Theorem 5-1-15\] and \[M, Main Theorem 2.1\]. We prove that a sequence of $`4`$-fold log flips for canonical pairs terminates after finitely many steps.

Let us recall the definition of canonical flips, which is slightly different from the usual one (cf. \[S, (2.11) Adjoint Diagram\]).

###### Definition 1 (Canonical flip).

Let $`X`$ be a normal projective variety and $`B`$ an effective $``$-divisor such that the pair $`(X,B)`$ is canonical, that is, $`\mathrm{discrep}(X,B)0`$. Let $`\varphi :(X,B)Z`$ be a small contraction corresponding to a $`(K_X+B)`$-negative extremal face. If there exists a normal projective variety $`X^+`$ and a projective morphism $`\varphi ^+:X^+Z`$ such that

1. $`\varphi ^+`$ is small;
2. $`K_{X^+}+B^+`$ is $`\varphi ^+`$-ample, where $`B^+`$ is the strict transform of $`B`$,

then we call $`\varphi ^+`$ the canonical flip or log flip of $`\varphi `$. We call the following diagram a flipping diagram:

$$\begin{array}{ccc}(X,B)& & (X^+,B^+)\\ \varphi & & \varphi ^+\\ & Z& & .\end{array}$$

The following is the main theorem of this paper.

###### Theorem 2 (Termination of $`4`$-fold canonical flips).

Let $`X`$ be a normal projective $`4`$-fold and $`B`$ an effective $``$-divisor such that $`(X,B)`$ is canonical. Consider a sequence of log flips starting from $`(X,B)=(X_0,B_0)`$:

$$(X_0,B_0)(X_1,B_1)(X_2,B_2)\mathrm{},$$

where $`\varphi _i:X_iZ_i`$ is a contraction and $`\varphi _{i}^{}{}_{}{}^{+}:X_{i}^{}{}_{}{}^{+}=X_{i+1}Z_i`$ is the log flip. Then this sequence terminates after finitely many steps.

###### Corollary 3 (\[M, Main Theorem 2.1\]).

Let $`X`$ be a projective $`4`$-fold with only terminal singularities and $`D`$ an effective $``$-Cartier $``$-divisor. Then any sequence of $`D`$-flops is finite.

We introduce a variant of difficulty. This is slightly different from \[FA, Chapter 4\]. It was inspired by \[FA, 4.14 Remark\]. Note that the notion of difficulty was first introduced by Shokurov in \[S, (2.15) Definition\].

###### Definition 4 (A weighted version of difficulty).

Let $`(X,B)`$ be a pair with only canonical singularities, where $`B=_{j=1}^lb_jB^j`$ with $`0<b_1<\mathrm{}<b_l1`$ and $`B^j`$ is a reduced divisor for every $`j`$. We note that $`B^j`$ is not necessarily irreducible. We put $`b_0=0`$, and $`S:=_{j0}b_j_0`$. Note that $`S=0`$ if $`B=0`$.

We call a divisor $`E`$ over $`X`$ essential if $`E`$ is exceptional over $`X`$ and is not obtained from blowing up the generic point of a subvariety $`WBX`$ such that $`B`$ and $`X`$ are generically smooth along $`W`$ (and thus only one of the irreducible components of $`_{j1}B^j`$ contains $`W`$) and $`dimW=dimX2`$. We set

$$d_{S,b}(X,B):=\underset{\xi S,\xi b}{}\mathrm{}\{E|E\text{is essential and}a(E,X,B)<1\xi \}.$$

Then $`d_{S,b_j}(X,B)`$ is finite by Lemma 5 below. Note that $`\mathrm{discrep}(X,B)0`$ since $`(X,B)`$ is canonical.

###### Lemma 5 (\[FA, (4.12.2.1)\]).

Let $`(X,B)`$ be a klt pair. Then

$$\mathrm{}\{E|E\text{is essential and }a(E,X,B)<\mathrm{min}\{1,1+\mathrm{discrep}(X,B)\}\}$$

is finite.

###### Definition 6.

Let $`\phi :(X,B)(X^+,B^+)`$ be a canonical flip. We say that this flip is of type $`(dimA,dimA^+)`$, where $`A`$ (resp. $`A^+`$) is the exceptional locus of $`\varphi :XZ`$ (resp. $`\varphi ^+:X^+Z`$). We call $`A`$ (resp. $`A^+`$) the flipping (resp. flipped) locus of $`\phi `$. When $`dimX=4`$, the log flip is either of type $`(1,2)`$, $`(2,2)`$ or $`(2,1)`$ by \[KMM, Lemma 5-1-17\].

###### Remark 7.

In \[KMM, §5-1\], the variety is $``$-factorial and every flipping contraction is corresponding to a negative extremal ray. However, the above properties were not used in the proof of \[KMM, Lemma 5-1-17\]. Note that \[KMM, Lemma 5-1-17\] holds in more general setting. For the details, see the original article \[KMM\].

###### Lemma 8 (cf. \[KM, Lemma 6.21\]).

Let $`\phi :(X,\mathrm{\Delta })(X^+,\mathrm{\Delta }^+)`$ be a canonical flips of $`n`$-folds. Let $`\mathrm{\Delta }:=\delta _i\mathrm{\Delta }_i`$ be the irreducible decomposition and $`\mathrm{\Delta }^+`$ (resp. $`\mathrm{\Delta }_i^+`$) the strict transform of $`\mathrm{\Delta }`$ (resp. $`\mathrm{\Delta }_i`$). Let $`F`$ be an $`(n2)`$-dimensional irreducible component of $`A^+`$, and $`E_F`$ the exceptional divisor obtained by blowing up $`F`$ near the generic point of $`F`$. Then $`X^+`$ is generically smooth along $`F`$ and

$$0a(E_F,X,\mathrm{\Delta })<a(E_F,X^+,\mathrm{\Delta }^+)=1\delta _i\mathrm{mult}_F(\mathrm{\Delta }_i^+),$$

where $`\mathrm{mult}_F(\mathrm{\Delta }_i^+)`$ is the multiplicity of $`\mathrm{\Delta }_i^+`$ along $`F`$.

###### Proof.

The pair $`(X^+,\mathrm{\Delta }^+)`$ is terminal near the generic point of $`F`$ by the negativity lemma. Therefore, $`X^+`$ is generically smooth along $`F`$, and the rest is an obvious computation. ∎

###### 9.

Let $`X`$ be a projective variety. Let $`X^{an}`$ be the underlying analytic space of $`X`$. Let $`H_k^{BM}(X^{an})`$ be the Borel-Moore homology. For the details, see \[F, 19.1 Cycle Map\]. Then there exists a cycle map;

$$cl:A_k(X)H_{2k}^{BM}(X^{an}),$$

where $`A_k(X)`$ is the group generated by rational equivalence classes of $`k`$-dimensional cycles on $`X`$. For the details about $`A_k(X)`$, see \[F, Chapter 1\]. We note that the cycle map $`cl`$ commutes with push-forward for proper morphisms, and with restriction to open subschemes. From now on, we omit the superscript $`an`$ for simplicity.

The following is in the proof of \[KMM, Theorem 5-1-15\].

###### Lemma 10 (Log flips of type $`(2,1)`$).

When the $`4`$-dimensional log flip

$$\begin{array}{ccc}(X,B)& & (X^+,B^+)\\ \varphi & & \varphi ^+\\ & Z& & \end{array}$$

is of type $`(2,1)`$, we study the rank of the $``$-module $`cl(A_2(X))`$, where $`cl(A_2(X))`$ is the image of the cycle map $`cl:A_2(X)H_4^{BM}(X)`$. By the following commutative diagrams;

$$\begin{array}{ccccccc}A_2(A)& & A_2(X)& & A_2(XA)& & 0\\ cl& & cl& & cl& & \\ H_4^{BM}(A)& & H_4^{BM}(X)& & H_4^{BM}(XA),\end{array}$$

and

$$\begin{array}{ccc}A_2(Z)& \stackrel{}{}& A_2(Z\varphi (A))\\ cl& & cl& & \\ H_4^{BM}(Z)& \stackrel{}{}& H_4^{BM}(Z\varphi (A)),\end{array}$$

we have the surjective homomorphism;

$$cl(A_2(X))cl(A_2(Z)).$$

We note that $`XAZ\varphi (A)`$ and \[F, p.371 (6) and Lemma 19.1.1\]. For any closed algebraic subvariety $`V`$ on $`X`$ of complex dimension $`2`$, $`V`$ is not numerically trivial since $`X`$ is projective. Therefore, $`cl(V)0`$ in $`cl(A_2(X))`$ (see \[F, Definition 19.1\]). Thus the kernel of the surjection above is not zero. By the similar arguments, we obtain that

$$cl(A_2(X^+))cl(A_2(Z)).$$

We note that $`A^+`$ is one-dimensional and $`X^+A^+Z\varphi ^+(A^+)`$. Therefore, since the rank $`\mathrm{rk}_{}cl(A_2(X))`$ is finite, we finally have the result

$$\mathrm{rk}_{}cl(A_2(X))>\mathrm{rk}_{}cl(A_2(X^+)).$$

Let us start the proof of Theorem 2.

###### Proof of Theorem 2.

Since $`d_{S,b_l}(X,B)`$ does not increase by log flips by the negativity lemma. Note that if $`E`$ is exceptional and $`a(E,X_i,B_i)<1b_l`$, then $`E`$ is essential. If $`B_i^l`$, which is the strict transform of $`B^l`$ on $`X_i`$, contains $`2`$-dimensional flipped locus, then $`d_{S,b_l}(X,B)`$ decreases by easy computations (see Lemma 8). So, after finitely many flips, $`B_i^l`$ does not contain $`2`$-dimensional flipped locus. Thus, we can assume that $`B_i^l`$ does not contain $`2`$-dimensional flipped locus for every $`i>0`$. Let $`\overline{B_i^l}`$ be the normalization of $`B_i^l`$. We consider a sequence of birational maps

$$\overline{B_0^l}\overline{B_1^l}\mathrm{}.$$

By \[M, Lemma 2.11\], after finitely many flips, $`B_i^l`$ does not contain $`2`$-dimensional flipping locus. Thus, we can assume that $`B_i^l`$ does not contain $`2`$-dimensional flipping locus for every $`i0`$. In particular, $`B_i^lB_{i+1}^l`$ is an isomorphism in codimension one for every $`i0`$.

Next, we see $`d_{S,b_{l1}}(X,B)`$. By the same argument as above, $`B_i^{l1}`$ does not contain $`2`$-dimensional flipped locus after finitely many flips. Note that if $`E`$ is essential over $`X_{i+1}`$ and $`a(E,X_{i+1},B_{i+1})<1b_{l1}`$, then $`E`$ is essential over $`X_i`$ and $`a(E,X_i,B_i)a(E,X_{i+1},B_{i+1})`$. It is because $`B_i^lB_{i+1}^l`$ is an isomorphism in codimension one. By using \[M, Lemma 2.11\], we see that $`B_i^{l1}`$ does not contain $`2`$-dimensional flipping locus after finitely many steps. Therefore, we can assume that $`B_i^{l1}B_{i+1}^{l1}`$ is an isomorphism in codimension one for every $`i0`$

By repeating this argument, we can assume that $`B_i^jB_{i+1}^j`$ is an isomorphism in codimenison one for every $`i`$, $`j`$.

If the log flip is of type $`(1,2)`$ or $`(2,2)`$, then $`d_{S,b_0}(X,B)=d_{S,0}(X,B)`$ decreases by Lemma 8. Therefore, we can assume that all the flips are of type $`(2,1)`$ after finitely many steps. This sequence terminates by Lemma 10. ∎

###### Remark 11 ($`3`$-fold case).

By using the weighted version of difficulty, we can easily prove the termination of $`3`$-fold canonical flips without using $``$-factoriality. This is \[FA, 4.14 Remark\].

###### Acknowledgments.

I was partially supported by the Inoue Foundation for Science. I talked about this paper at RIMS in May 2002. I would like to thank Professors Shigefumi Mori, Shigeru Mukai, and Noboru Nakayama, who gave me some comments. I am also grateful to Mr. Kazushi Ueda, who answered my question about Borel-Moore homology.


text:
# Complete identification of alkali sites in ion conducting lithium silicate glasses: a computer study of ion dynamics

## Abstract

The available sites for ions in a typical disordered ionic conductor are determined. For this purpose we devised a straightforward algorithm which via cluster analysis identifies these sites from long time ionic trajectories below the glass transition. This is exemplified for a lithium silicate glass $`(\mathrm{Li}_2\mathrm{O})_x(\mathrm{SiO}_2)_{(1x)}`$ for $`x=0.5`$ and $`x=0.1`$. We find for both concentrations that the number of sites is only slightly bigger than the number of ions. This result suggests a theoretical description of the dynamics in terms of mobile vacancies as most appropriate. Furthermore identification of the ionic sites allows one to obtain detailed characteristics of the ionic motion, e.g. quantification of correlated forward-backward jumps.

Ion conducting glasses have been investigated by various experimental methods, including EXAFS Greaves et al. (1991); Rocca (1992), NMR Maekawa et al. (1991); Yap et al. (1995); Gee and Eckert (1996), and conductivity spectroscopy Funke (1993); Abe et al. (1993). Whereas quite detailed information about the local structure has become available the mechanism of dynamics is still under debate although consensus has been reached that ion dynamics can be described as jumps of the mobile ions Ingram (1987); Baranowskii and Cordes (1999); Bunde et al. (1994). Whereas some authors stress the relevance of the disordered energy landscape Baranowskii and Cordes (1999); Dyre and Schrøder (2000) supplied by the network, others relate the complexity of ion dynamics to the Coulomb interaction among mobile ions Funke (1993). Furthermore it has been argued from structural considerations that the distribution of alkali ions is inhomogeneous Greaves (1985).

For a closer understanding of ion dynamics microscopic information as supplied by molecular dynamics (MD) simulations is highly welcome. Jund et al. Jund et al. (2001) found preferential pathways as a dynamical phenomenon by counting the number of different alkali ions that passed through subvolumes of the simulation box. The resulting subvolumes visited by the largest number of different alkali ions may be interpreted as fast areas. They form a network of conduction paths. Oviedo and Sanz have argued on a qualitative basis that for alkali concentrations lower than 10% the alkali ions are always surrounded by non-bridging oxygens (NBOs) despite the overall small number of NBOs Oviedo and Sanz (1998). This implies that new NBOs are formed via breaking of Si-O bonds along the alkali trajectory. In contrast, for higher concentrations hopping dynamics between so-called micro-channels is proposed such that no formation of new NBOs would be necessary. Horbach Horbach et al. (2002) showed in simulations of $`(\mathrm{Na}_2\mathrm{O})2(\mathrm{SiO}_2)`$, that the coherent intermediate scattering function of sodium relaxes only on the same timescale as those of the network species, proving the existence of stable alkali sites below the glass transition. Cormack Cormack et al. (2002) et al. have recently investigated the mechanism of sodium migration in simulations of $`(\mathrm{Na}_2\mathrm{O})_{0.25}(\mathrm{SiO}_2)_{0.75}`$ glasses, observing a few sequences of jumps between selected sites. He interpreted the resulting dynamics as the motion of vacancies and pointed out that the identification of all sites in the glass, empty or populated at a given time, would be most useful for obtaining a deeper understanding of the kinetic mechanisms.

Determination of all alkali sites in the silicate network is the purpose of this work. These sites can be defined as any regions in the network where alkali ions stay. They must correspond to local minima of the potential energy for the ions. Thus determination of the ionic sites is equivalent to determining the effective energy landscape experienced by the ions. The potential energy involves the network contribution as well as the averaged interaction with other ions.

We analyse lithium silicate glasses $`(\mathrm{Li}_2\mathrm{O})_x(\mathrm{SiO}_2)_{(1x)}`$ with $`x=0.5`$ and $`x=0.1`$ which are common model systems for computer investigations of ion dynamics. The system size has been chosen such that for $`x=0.5`$ we have 384 lithium ions and for $`x=0.1`$ 80 lithium ions. Simulations are performed in the NVT ensemble at temperatures $`640\text{K}`$ and $`920\text{K}`$ and simulation box lengths $`22.96`$ Å and $`25.28`$ Å, respectively. The resulting densities match experimental values Doweidar (1996). The temperatures have been chosen such that the lithium diffusion constants in both systems have similar values. We use a pair potential developed by J. Habasaki Habasaki and Okada (1992). Previous studies with this potential have shown good agreement with experimental results for static and dynamic quantities Habasaki et al. (1995); Banhatti and Heuer (2001); Heuer et al. (2002). The molecular dynamics simulations are performed with a modified version of Moldy Refson (2000). Since both temperatures are below the respective glass transition temperatures we proceed in two steps. We start from a configuration which was first equilibrated just above the glass transition. It is then propagated at the selected temperature for ca. 20 ns in the NVT ensemble until the lithium subsystem is equilibrated. Afterwards we perform a production run for approximately 50 ns. During this run, positions of all particles are recorded every 0.1 ps. More details about the simulations can be found in Heuer et al. (2002). We checked the validity of our results for independent data of up to 10 ns duration. This time is long enough for the lithium dynamics to become diffusive, so that just the dynamics relevant for ionic transport is included.

In Fig. 1 we report the mean square displacement of oxygen and lithium. The dynamics of the lithium ions becomes diffusive for $`t_{diff}1`$ns for $`x=0.5`$ and $`t_{diff}10`$ ns for $`x=0.1`$, as obtained from longer simulations. Evidently, the oxygen ions are localised on the ns-time scale. The network of silicon and oxygen atoms and with this the lithium sites can therefore be considered as static. As a consequence the ionic sites are basically unmodified during our run. Interestingly, as shown recently Sunyer et al. (2002); Angell et al. (1981) the remaining local fluctuations of the network are nevertheless essential in promoting the ionic jumps.

Based on the lithium trajectories we have devised a straightforward algorithm to identify ionic sites. As a basic idea we let the ions decide via their dynamics where good sites are rather than determining the potential energy, as done e.g. in Habasaki et al. (1996). First the simulation box is divided into cubic cells with size approx. (0.3 Å)3. They are small enough to enable resolution of the shape of the ionic sites. During our MD-run we count how many time steps a cell is visited by a lithium ion. The cells with nonzero counts describe the portion of the system that has been visited by lithium ions. These cells include the lithium sites as well as the connecting paths between them. To identify the ionic sites and to eliminate the paths between the sites, cells with less than a minimum value $`M`$ of counts are dismissed. With the remaining cells a cluster analysis is performed. Cells which share a face are grouped into one cluster. Cells that are not directly or indirectly connected over common faces form different clusters. The value of $`M`$ has been determined by the condition that the number of distinct ionic clusters be a maximum. For a smaller threshold some lower populated clusters would not have been detected, for a larger threshold distinct clusters would be connected to single clusters because the transition path would be included. With $`M`$ chosen according to this condition, 2.9% and 1.3% of all cells are included for $`x=0.5`$ and $`x=0.1`$, which contain $`90`$% of all lithium positions in both cases. The precise value of $`M`$ is of no relevance for the latter analysis.

For $`x=0.5`$ and $`x=0.1`$ respectively, 378 and 76 clusters are found. The result is illustrated in Fig. 2, showing a snapshot of a slice through each system. The silicate network is shown as tubes, with silicon bright and oxygens dark. Lithium ions are depicted as large spheres. The white objects are the clusters obtained from the analysis described above, with each small white sphere representing one of the cells. It can be seen that the clusters are in fact discrete objects of compact shape. Most of the lithium ions reside inside of one of these clusters, and most of the clusters are occupied by exactly one lithium ion. As will be quantified further below some clusters can accommodate up to three lithium ions.

For x=0.5 (x=0.1) in total 29 (14) additional small clusters have been excluded. They are set apart by a conspicious drop in the effective radius plot; see Fig. 3. Detailed analysis characterised most of them as simple satellites of bigger clusters. The remaining ones serve as short–lived saddle–like states during transitions between bigger clusters. The small clusters need thus not be considered as sites themselves.

On the basis of the clusters, the dynamics of a lithium ion can be expressed as successive residences in clusters separated by shorter time intervals where the lithium ion is outside of any cluster. When a lithium ion leaves one cluster and moves into a different one, this sequence of events is recorded as a jump. In contrast, if an ion explores part of the volume outside a cluster without entering a second one and comes back to the original cluster, it is regarded as occupying the same cluster during the whole time. Consequently, the residence time of a lithium ion in a cluster is defined as the time between its jump into this cluster and the subsequent jump to another cluster. With this definition on average a lithium ion in the system with $`x=0.5`$ belongs to a cluster at 99.5% of all times. The remaining 0.5% correspond to jumps between different sites. For $`x=0.1`$, the portion of times corresponding to jumps is 1.1%. The typical distance travelled by the lithium ions during a jump coincides well with the first maximum of the radial distribution function of the clusters in both cases. For $`x=0.5`$, 93% of all jumps are inside the first neighbour shell (94% for $`x=0.1`$) as defined by the first minimum of the radial distribution function. Thus the ion dynamics can be described as jumps between adjacent clusters.

Having defined the residence times the mean number of lithium ions occupying each cluster can easily be evaluated. It is plotted in Fig. 3 for $`x=0.5`$, with the clusters ordered from left to right by decreasing volume $`V`$ (the effective radius $`r_{eff}(3V/4\pi )^{1/3}`$ of the individual clusters is plotted in the inset). The occupation data confirms the impression gained from Fig. 2. For most clusters the value is close to unity. These clusters are occupied by a single ion during the largest part of the simulation. Their effective radius $`r_{eff}0.6`$Å is very similar. Only the largest clusters offer space for more than one lithium ion, i.e. they contain more than one site, and often they are indeed multiply occupied.

Generally an empty site must be available for a jump of a lithium ion. The total number of sites can now be estimated from the occupation data in Fig. 3. For the larger clusters with cluster index $`<19`$ we have taken the next higher integer number to estimate the number of available sites. In this way we find 3 clusters with 3 sites and 9 clusters with 2 sites. Starting from cluster index 14 one finds a few clusters with occupation number 1 + $`ϵ`$ with $`ϵ<0.2`$. Here it would be rather unphysical to attribute two sites to this cluster. Rather a more detailed analysis of explicit ion trajectories has revealed the following scenario. At most time instances a single ion stayed in this cluster. During the short times during which this ion explored the immediate neighborhood of this cluster without entering a new cluster (typical excursion length scales are 2 Å from the center of the cluster) a second ion may briefly enter this cluster and immediately leave this cluster again. Thus one may either say that these clusters only contain a single site or, to be more conservative with respect to finding an upper bound of the total number of sites, attribute $`1+ϵ`$ sites to this cluster. With the latter variant we obtain 395 sites for our 384 lithium ions. Analogously, 86 sites were found for the 80 ions for $`x=0.1`$. These results shed new light on the mechanism of ion dynamics. Rather than speaking of individual ions jumping between accessible sites it is more appropriate to speak to first approximation of vacancy dynamics, i.e. the dynamics of the non–interacting unoccupied sites in the energy landscape.

Actually, Dyre Dyre (2002) has recently suggested this type of scenario based on general reasoning. He relates the possible near-equality of sites and ions to the effect that during network formation at the glass transition only a minimum number of sites will be formed due to energetic reasons. This scenario is also compatible with the counter-ion model, proposed by Dieterich and coworkers Dieterich et al. (1994), but is at variance with single-particle approaches.

In the remaining part of this paper we would like to characterize the ion dynamics in terms of the individual clusters. Fig. 4 shows the cumulative distribution function $`S(\tau )`$ of residence times. Thus $`S(\tau )`$ denotes the fraction of clusters with an average residence time less than $`\tau `$. Its median $`\tau _{median}`$, defined by $`S(\tau _{median})=0.5`$, is close to 1 ns for $`x=0.5`$ and 250 ps for $`x=0.1`$. These numbers have to be compared with the median of the transition time for a jump between two clusters, which is 500 fs for $`x=0.5`$ and 200 fs for $`x=0.1`$. This dramatic time scale separation clearly justifies the hopping picture in this kind of materials, as known from previous simulations Habasaki et al. (1997). The lithium dynamics on long time scales can thus be mapped to jumps between the clusters. The generally higher frequency and lower duration of jumps at $`x=0.1`$ can be traced to the higher temperature that was chosen for this system to achieve similar diffusivities.

The effectiveness of jumps for diffusive transport is strongly decreased by correlated forward-backward motion. We determined the probability $`p_{back}`$ that a jump from a cluster A to another cluster B is followed by a direct backjump to the cluster A. The result is also shown in Fig. 4. One can clearly see that clusters with short residence times show a stronger tendency for a correlated forward-backward jump. For the fastest clusters, $`p_{back}`$ reaches more than 0.9, which is several times higher than the statistical value given by the inverse number of nearest neighbors $`(0.2\mathrm{}0.3)`$. The difference in $`p_{back}(\tau _{median})`$, which is $`0.8`$ for $`x=0.1`$ compared to $`0.55`$ for $`x=0.5`$, may explain why the former system becomes diffusive only after most ions have performed several jumps, as shown by the ratio $`\tau _{median}/t_{diff}>>1`$. In contrast one finds $`\tau _{median}/t_{diff}1`$ for $`x=0.5`$, meaning that only half of the ions have jumped at least once at the onset of diffusion. The dispersive regime at shorter times is thus dominated by the fastest lithium ions.

Beyond the results shown in Fig. 4, a detailed analysis of the dynamics, i. e. in terms of cooperativity, becomes feasible. Our approach may easily be generalised to more complex systems and important questions concerning, e.g., the character and the lifetime of the individual sites may be analysed.

We acknowledge interesting discussions with B. Doliwa, J. Dyre, M. Ingram, P. Jund, B. Roling, and M. Vogel and very helpful correspondence with J. Habasaki about this topic.