text:
# Radiative transfer through the Intergalactic Medium

## 1 Introduction

Hydrodynamical simulations of structure formation in the universe have led to fundamental insights into the structure and evolution of the Intergalactic Medium (IGM) (Cen et al. 1994; Zhang, Anninos & Norman 1995; Hernquist et al. 1996; Zhang et al. 1997; Bond & Wadsley 1997; Theuns, Leonard & Efstathiou 1998). Comparisons with the Ly$`\alpha `$forest, as measured in high redshift Quasi-Stellar Object (QSO) spectra, show that the simulations broadly reproduce the statistical properties of the IGM. Precision comparisons with the highest spectral resolution measurements recover the cumulative flux distributions and $`\mathrm{I}`$column density distributions to an accuracy of a few percent (Meiksin, Bryan & Machacek 2001). More problematic are the predicted widths of the absorption features, which appear to require additional sources of broadening, perhaps late He $`\mathrm{II}`$re-ionization, to match the measured widths (Theuns et al. 1999; Bryan & Machacek 2000; Meiksin et al. 2001).

Many of the properties of the IGM may be attributed to the gravitational instability of the dark matter alone, as the baryon density fluctuations closely follow those of the dark matter (Zhang et al. 1998). This has led to the development of pseudo-hydrodynamic simulations based on pure gravity in which the baryon fluctuations are assumed to exactly follow the dark matter and the temperature to be derived from an effective equation of state (Petitjean, Mücket & Kates 1995; Croft et al. 1998; Gnedin & Hui 1998; Meiksin & White 2001). In recent years, such simulations have been increasingly relied on for predicting the flux power spectrum of the Ly$`\alpha `$forest (Meiksin & White 2001; Zaldarriaga, Hui & Tegmark 2001; Croft et al. 2002; Meiksin & White 2003b; Seljak, McDonald & Makarov 2003).

The simulations have generally been done assuming sudden early homogeneous $`\mathrm{I}`$and He $`\mathrm{II}`$re-ionization. While tentative steps have been taken to introduce radiative transfer into the simulations using approximate schemes (Abel, Norman & Madau 1999; Gnedin & Abel 2001; Nakamoto, Umemura & Susa 2001; Ciardi, Stoehr & White 2003), these simulations have emphasized the propagation of $`\mathrm{I}`$ionization fronts and the predicted mean optical depths. Except for the simulations of Gnedin & Abel, they have not solved fully self-consistently for the combined gas dynamics and radiative transfer.

In this Letter, we extend the photon-conserving scheme of Abel et al. to include the treatment of helium and explore the thermal effects of radiative transfer in the inhomogeneous medium predicted in a $`\mathrm{\Lambda }`$CDM cosmology.

## 2 Radiative transfer

### 2.1 Simulation data

We use the data from a previously run $`\mathrm{\Lambda }`$CDM simulation used to investigate effects of radiative transfer on the Ly$`\alpha `$forest (Meiksin & White 2003b). The simulation was run using a pure Particle Mesh (PM) dark matter code, and it was assumed the gas and dark matter have the same spatial distribution. A description of the parallel PM code is given in Meiksin & White (2003a). The parameters used for the simulation are $`\mathrm{\Omega }_\mathrm{M}=0.30`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.70`$, $`\mathrm{\Omega }_b=0.045`$, $`h=H_0/100`$ km s-1 Mpc-1$`=0.70`$, and slope of the primordial density perturbation power spectrum $`n=1.05`$. The model is consistent with existing large-scale structure, Ly$`\alpha `$forest flux distribution, cluster abundance and WMAP constraints (Meiksin & White 2003b). The simulation was run using $`512^3`$ particles and a $`1024^3`$ force mesh, in a cubic box with (comoving) side length $`25h^1`$Mpc, adequate for obtaining converged estimates of the Ly$`\alpha `$pixel flux distribution and flux power spectrum (Meiksin & White 2003a,b). A typical line-of-sight density run at $`z=6`$, used for the radiative transfer calculations here, is shown in Fig. 1.

### 2.2 Equations of radiative transfer

To consider the effect of radiative transfer when modelling the propagation of ionization fronts (I-fronts) into the IGM, we extend the photon-conserving algorithm of Abel et al. to include helium. The advantage of this scheme is that energy is conserved independently of numerical resolution, ensuring that I-fronts propagate at the correct speeds in our simulations. In practice, this means that larger step sizes may be taken on the simulation space grid without the associated loss of accuracy which would occur when solving the radiative transfer equation through direct numerical integration.

The probabilities for the absorption of an ionizing photon by $`\mathrm{I}`$, He $`\mathrm{I}`$and He $`\mathrm{II}`$, respectively, are

$`P_{\mathrm{abs}}^{\mathrm{HI}}`$ $`=`$ $`p_{\mathrm{HI}}q_{\mathrm{HeI}}q_{\mathrm{HeII}}\left[1\mathrm{exp}\left(\tau _\nu ^{\mathrm{tot}}\right)\right]/D,`$ (1)
$`P_{\mathrm{abs}}^{\mathrm{HeI}}`$ $`=`$ $`q_{\mathrm{HI}}p_{\mathrm{HeI}}q_{\mathrm{HeII}}\left[1\mathrm{exp}\left(\tau _\nu ^{\mathrm{tot}}\right)\right]/D,`$ (2)
$`P_{\mathrm{abs}}^{\mathrm{HeII}}`$ $`=`$ $`q_{\mathrm{HI}}q_{\mathrm{HeI}}p_{\mathrm{HeII}}\left[1\mathrm{exp}\left(\tau _\nu ^{\mathrm{tot}}\right)\right]/D.`$ (3)
Here we have defined the auxiliary absorption and transmission probabilities $`p_i=1\mathrm{exp}(\tau _\nu ^i)`$, $`q_i=\mathrm{exp}(\tau _\nu ^i)`$, $`i`$ denoting the species being referred to, $`\tau _\nu ^{\mathrm{tot}}=\tau _\nu ^{\mathrm{HI}}+\tau _\nu ^{\mathrm{HeI}}+\tau _\nu ^{\mathrm{HeII}}`$, where $`\tau _\nu ^i`$ is the optical depth for a given species, and $`D=p_{\mathrm{HI}}q_{\mathrm{HeI}}q_{\mathrm{HeII}}+q_{\mathrm{HI}}p_{\mathrm{HeI}}q_{\mathrm{HeII}}+q_{\mathrm{HI}}q_{\mathrm{HeI}}p_{\mathrm{HeII}}`$. These probabilities can be used to calculate the ionization rate for a given species per unit volume, $`n_i\mathrm{\Gamma }_i`$, as follows. If in a time $`\delta t`$, $`\delta t\dot{N}_\nu ^{l1}`$ photons enter grid zone $`l`$ from zone $`l1`$, then the number of photons that will be absorbed in zone $`l`$ by species $`i`$ is $`\delta t\dot{N}_\nu ^{l1}P_{\mathrm{abs}}^i`$, where $`P_{\mathrm{abs}}^i`$ is the absorption probability for species $`i`$ within zone $`l`$. The ionization rate of species $`i`$ in zone $`l`$ of volume $`V^l`$ is then

$$n_i^l\mathrm{\Gamma }_i^l=\frac{1}{V^l}\underset{g}{}\dot{N}_{\nu _g}P_{\mathrm{abs}}^i(\nu _g),$$
(4)
where we have divided the frequencies into 100 discrete groups $`g`$ evenly spaced between $`\nu _L^{\mathrm{HI}}`$ and $`10\nu _L^{\mathrm{HI}}`$. (While we found this spacing to be adequate, we have made no attempt to optimize it.) These are used in the following set of coupled equations to solve for the positions of the three I-fronts over time:

$$\frac{dn_{\mathrm{HII}}}{dt}=n_{\mathrm{HI}}\mathrm{\Gamma }_{\mathrm{HI}}n_en_{\mathrm{HII}}\alpha _{\mathrm{HI}}(T)3\frac{\dot{a}}{a}n_{\mathrm{HII}},$$
(5)
$`{\displaystyle \frac{dn_{\mathrm{HeII}}}{dt}}`$ $`=`$ $`n_{\mathrm{HeI}}\mathrm{\Gamma }_{\mathrm{HeI}}+n_en_{\mathrm{HeIII}}\alpha _{\mathrm{HeII}}(T)n_{\mathrm{HeII}}\mathrm{\Gamma }_{\mathrm{HeII}}`$ (6)
$``$ $`n_en_{\mathrm{HeII}}\alpha _{\mathrm{HeI}}(T)3{\displaystyle \frac{\dot{a}}{a}}n_{\mathrm{HeII}},`$
$$\frac{dn_{\mathrm{HeIII}}}{dt}=n_{\mathrm{HeII}}\mathrm{\Gamma }_{\mathrm{HeII}}n_en_{\mathrm{HeIII}}\alpha _{\mathrm{HeII}}(T)3\frac{\dot{a}}{a}n_{\mathrm{HeIII}},$$
(7)
where $`n_i`$ denotes number density, $`\mathrm{\Gamma }_i`$ $`(\mathrm{s}^1)`$ is the photoionization rate per atom, $`\alpha _i(T)`$ is the total radiative recombination coefficient, and $`a`$ is the cosmological expansion factor.

The time evolution of the gas temperature $`T`$ is given by:

$$\frac{dT}{dt}=\frac{2(GL)}{3kn}+\frac{T}{n}\frac{dn_e}{dt}3\frac{\dot{a}}{a}\left(\frac{2}{3}T+\frac{n_e}{n}T\right).$$
(8)
where $`n=n_{\mathrm{HI}}+n_{\mathrm{HeI}}+n_{\mathrm{HII}}+n_{\mathrm{HeII}}+n_{\mathrm{HeIII}}+n_e`$, $`n_e=n_{\mathrm{HII}}+n_{\mathrm{HeII}}+2n_{\mathrm{HeIII}}`$, $`k`$ is the Boltzmann constant, and $`G`$ and $`L`$ $`(\mathrm{J}\mathrm{m}^3\mathrm{s}^1)`$ are the atomic heating and cooling rates, respectively. The last term is the adiabatic cooling term resulting from cosmological expansion, which will dominate the thermal effects of gas motions at the densities we consider. While computing the photoionization, we assume the gas over-density (not the density) stays frozen, which is a good approximation on the scales relevant to the Ly$`\alpha `$forest (Zhang et al. 1998). As in similar previous studies, we do not include the effects of adiabatic compression heating as this will only affect the temperature for the relatively rare virialised halos with circular velocities exceeding $`20`$ km s-1 (Meiksin 1994; Meiksin & White 2003b).

The heating rate $`G^l`$ in cell $`l`$ is due to the photoionization of $`\mathrm{I}`$, He $`\mathrm{I}`$and He $`\mathrm{II}`$according to $`G^l=G_{\mathrm{HI}}^l+G_{\mathrm{HeI}}^l+G_{\mathrm{HeII}}^l`$, where for each species $`i`$, $`G_i^l`$ is evaluated in a similar manner to the ionization rate per unit volume as given in Eq. (4),

$$n_i^lG_i^l=\frac{1}{V^l}\underset{g}{}(h\nu _g\chi _i)\dot{N}_{\nu _g}P_{\mathrm{abs}}^i(\nu _g),$$
(9)
where $`\chi _i`$ is the ionization potential of species i. The atomic cooling rate $`L`$ includes radiative recombination cooling to $`\mathrm{I}`$, He $`\mathrm{I}`$and He $`\mathrm{II}`$, electron excitation of $`\mathrm{I}`$, and Compton cooling off the Cosmic Microwave Background photons. The radiative recombination and cooling rates are taken from Meiksin (1994).

Eqs. (5)–(8) are solved using explicit forward Euler integration. Although an implicit numerical scheme such as backward Euler integration results in a more stable solution at low numerical resolution (Anninos et al. , 1997), we find that to obtain the required numerical accuracy as well as to maintain stability, both methods require similar numerical resolution. For the sake of speed and simplicity, the forward method is preferred. The timestep is restricted to being no more than several times greater that the hydrogen ionization timescale, $`\mathrm{\Gamma }_{\mathrm{HI}}^1`$, for numerical accuracy.

We take the frequency specific luminosity of the ionizing source in our simulations to be described by a power law spectrum with index $`\alpha =1.5`$, and a luminosity of $`L_\nu ^Q=10^{23}\mathrm{W}\mathrm{Hz}^1`$ at the Lyman edge. This spectrum is typical of that of QSOs, which are probable candidate sources for contributing to the re-ionization of the IGM. The assumed mass fraction of helium in the IGM is $`Y0.235`$.

We tested the photon-conserving algorithm on the photoionization of gas with uniform density around a point source. The resulting solutions for the ionization and temperature profiles were accurate to better than 10 per cent for optical depths at the Lyman edge of up to $`\mathrm{\Delta }\tau _\nu 20`$ per cell on the space grid. This results in a significant reduction in the computational resources compared with a radiative transfer scheme in which the ionization rates are computed by direct numerical integration over the intensity ($`L_\nu ^Qe^{\tau _\nu }/[4\pi r]^2`$) and photoionization cross-sections (eg, Madau, Meiksin & Rees 1997). Typically, we find that to obtain an accuracy comparable to the photon-conserving scheme at $`\mathrm{\Delta }\tau _\nu ^{\mathrm{HI}}20`$, a direct integration algorithm requires $`\mathrm{\Delta }\tau _\nu ^{\mathrm{HI}}1/4`$ in each cell. For the simulation results presented in this Letter, the original PM grid of 1024 separate cells was used over a line of sight 25$`h^1`$ comoving Mpc in length. This provided accurate spatial resolution, as at most $`\mathrm{\Delta }\tau _\nu ^{\mathrm{HI}}6.8`$ in each cell at the average baryonic density assumed for the neutral IGM at $`z=6`$.

## 3 Results

We set up a problem of a QSO source of luminosity $`L_\nu ^Q`$ at $`z=6`$ placed at a comoving distance of 10$`h^1`$ Mpc from the left edge of the density run shown in Fig. 1, and assume the gas surrounding the QSO up to that point has already been ionized. Placing the source at this distance avoids having to impose the restriction that the I-fronts propagate no faster than the speed of light. We indicate the displaced position of the source in the figures by starting the (comoving) spatial axes at $`R=10h^1`$ Mpc.

Two cases were compared: computing the ionization and temperature profiles by incorporating radiative transfer as discussed in section 2.2 (hereafter the RT simulation), and a second case neglecting radiative transfer, assuming instead an instantaneous uniform ionization rate across the whole line of sight (hereafter the nRT simulation), as is usually done in simulations of the Ly$`\alpha `$forest. Our intention is to compare the final gas temperatures to determine how large an effect including radiative transfer may have.

Fig. 2 shows a comparison between the IGM temperatures at $`z=5`$ computed with and without radiative transfer. The inclusion of radiative transfer results in a significant boost to the temperature of the ionized IGM, a point illustrated by Abel & Haehnelt (1999) in the context of a smooth IGM. The positions of the $`\mathrm{II}`$and He $`\mathrm{II}`$I-fronts at $`z=5`$ in the RT simulation are 24$`h^1`$ comoving Mpc from the source, while the He $`\mathrm{III}`$front lags behind at about 19$`h^1`$ comoving Mpc; their speeds of propagation are limited by the reduction in the ionization rate due to the increasing optical depth through the box. Consequently, the heating of the IGM, which we find to be dominated by the photoionization of neutral hydrogen at the $`\mathrm{II}`$I-front, is also constrained to lie behind the $`\mathrm{II}`$I-front. All gas beyond 24$`h^1`$ Mpc is still much cooler. In contrast, the temperature computed by the nRT simulation is much lower and the heating extends across the whole line of sight. The temperature is also much less sensitive to variations in the gas density for the nRT simulation. Fig. 3 compares the results obtained at $`z=3`$. By this time, both the $`\mathrm{II}`$and He $`\mathrm{II}`$I-fronts have reached 35$`h^1`$ Mpc from the source. The material which has been ionized early on in the RT simulation is now beginning to cool to the same temperature as that computed by the nRT simulation, although the temperature differences remain large. In particular, the gas that was photoionized at $`z=5`$ ($`R<24h^1\mathrm{Mpc}`$), remains nearly twice as hot as the uniformly ionized gas by $`z=3`$. By contrast, the regions which were ionized more recently are much hotter than for the nRT simulation. We find that only by $`z=1`$, once the gas has been in an ionized state for a sufficiently long period of time, do the temperatures predicted in the two scenarios converge.

The reason for the temperature boost in the RT simulation is made clear by considering the energy per ionization ($`E=G/\mathrm{\Gamma }`$) at the I-front. In Fig. 4, we show that the average energy of a photon absorbed by either $`\mathrm{I}`$, He $`\mathrm{I}`$or He $`\mathrm{II}`$increases at the relevant I-front, as proportionally more of the lower energy photons have been absorbed already – a consequence of including radiative transfer. This increase to the heating energy per ionization is responsible for the temperature boost. The assumption of an instantaneous ionization field across the line of sight results in an underestimate of the IGM temperature at high redshift.

A comparison of the IGM temperature-density relation computed by the two different simulations is made in Figs. 57. Only data points where He $`\mathrm{II}`$has been 90 per cent ionized are plotted, so that gas which is neutral or in the process of being ionized is not considered. In Fig. 5, the underdense regions in the line of sight (Fig. 1), are initially heated to a higher temperature in the RT simulation. A non-monotonic temperature-density relation results. The relation actually splits into two trends. To understand the origin of the split, we divide the data points into two distinct sets corresponding to two regions in which the average density fluctuations are slightly higher or lower than the overall average in the line of sight (see caption to Fig. 5). The gas in the denser region, which is closer to the source and thus is ionized sooner, shows the same overall trend as the gas in the less dense region but is slightly cooler at a fixed over-density. In contrast, in the nRT simulation the temperature-density relation is both monotonic and essentially single-valued. By $`z=3`$ (Fig. 6), the ionized gas has cooled adiabatically, with the overdense regions cooling at a slower rate than the underdense regions, as the denser regions are better able to maintain thermal balance because of the shorter cooling (and photoionization heating) times. The IGM temperatures computed by the two simulations are now beginning to converge at high density, but at low density there is still a significant difference between the results. A third set of data points is also evident, corresponding to the latest locally underdense region to be ionized.

Finally, by $`z=1`$, the gas temperatures computed by the RT and nRT simulations have nearly converged. The gradient of the line in the temperature-density plane indicates that adiabatic cooling due to expansion is largely responsible for cooling the gas. We find the entropy is now nearly constant throughout the line of sight. The striations in the RT data points have also disappeared. By $`z=1`$, the results from the two simulations are in good agreement.

## 4 Conclusions

We extend the probabilistic method of radiative transfer described for hydrogen by Abel et al. (1999) to include helium and tested the method on the propagation of an ionization front through a uniform medium. We compare the method against direct numerical integration of the exact time-steady radiative transfer equation. We confirm their conclusion that the probabilistic method accurately reproduces the position of the ionization fronts even for large optical depths at the photoelectric edges. Specifically, we recover the numerically integrated converged solutions for the ionization fractions and temperature to an accuracy of 10 per cent for an incremental optical depth per grid zone at the $`\mathrm{I}`$photoelectric edge as high as 20, while direct numerical integration of the radiative transfer equations requires an incremental optical depth less than 1/4 to obtain comparable accuracy.

We use the probabilistic method to compute the propagation of an I-front generated by a QSO source through the IGM, using the density distribution drawn from a numerical simulation in a $`\mathrm{\Lambda }`$CDM universe. Including the effects of radiative transfer results in a substantial boost in the post-ionized gas temperature compared with the case of sudden uniform ionization with no radiative transfer. The boost is due chiefly to the increased energy per photoionization within the ionization fronts, a consequence of the higher probability for high energy photons to pass through the gas compared with lower energy photons.

A consequence of the radiative transfer is to alter the dependence of temperature on density from the approximately polytropic relation found for moderate to low densities assuming uniform photoionization (Meiksin 1994; Gnedin & Hui 1998) to a generally non-monotonic trend. While the differences between the temperatures for over-densities above $`5`$ are small (less than 25 per cent), at lower densities the two temperatures converge only slowly with time, nearly recovering the polytropic relation only by $`z1`$.

At $`z>1`$, the relation between temperature and density is not single-valued: multiple temperatures may result for different gas parcels of the same over-density. Although this is in part a consequence of the delay in time for gas at different distances from the source to be photoionized as the ionization front propagates, an additional contributing factor is the environment of the gas parcel: if the ionization front has passed through a region of dense gas, sufficient to drive the optical depth at the $`\mathrm{I}`$photoelectric edge to near unity, before reaching the parcel, the radiation field will deliver a higher energy per photoionization to the parcel. Thus the determination of the gas temperature is strongly non-local due to the dependence of the heating rate on the optical depth to photoionizing photons between the gas parcel and the source of ionizing radiation.

Although the inclusion of radiative transfer results in a substantial and persistent boost in gas temperature, it is insufficient in itself to account for the broader measured Ly$`\alpha `$absorption systems than predicted by numerical simulations assuming sudden uniform photoionization. Meiksin, Bryan & Machacek (2001) find that a boost of $`\mathrm{\Delta }T1.7\times 10^4\mathrm{K}`$ is required at $`z\stackrel{<}{}3.5`$, larger than the temperatures we find here. Although it may be possible to obtain a larger boost for an alternative re-ionization scenario, such a large boost may still require invoking late He $`\mathrm{II}`$re-ionization. Our principal conclusion here is that for numerical simulations to accurately predict the temperature of the IGM, they must include the effects of radiative transfer. We are currently extending our methods to incorporate radiative transfer with ray tracing (eg, Abel & Wandelt 2002), into fully self-consistent simulations to explore this and related issues, including alternative sources of re-ionization like low luminosity AGNs and stars.

A.M. thanks the University of Edinburgh Development Trust for its financial support.


format: text

text:
# Comment on ”Dielectric behavior of paraelectric KTaO3, CaTiO3, and (Ln1/2Na1/2)TiO3 under a dc electric field”

## Abstract

Chen Ang, Bhalla and Cross \[Phys. Rev. B 64, 184104 (2001)\] have studied the low-frequency (20 Hz – 100 kHz) dielectric dispersion of the KTaO3 crystal under a dc electric field. Performing fits of the electric field dependence they came to conclusion that an appreciable contribution to the dielectric permittivity originates from polar clusters. In this Comment we show that the dielectric permittivity at low frequencies (100 Hz – 1 MHz) equals to that in the THz region, close below the polar phonon response. This excludes the possibility of any appreciable dielectric dispersion due to polar clusters. In addition, we demonstrate that correct treatment using Landau-Ginzburg-Devonshire theory allows to fit the electric field dependence of the dielectric constant without assuming any polarization mechanism besides the polar phonon modes.

The family of incipient ferroelectric materials with perovskite structure has attracted a great attention due to its high value of dielectric constant, low dielectric losses and high degree of tunability by the external electric field and temperature. The best known members of this family are SrTiO3 (STO), KTaO3 (KTO)Samara (2001) and CaTiO3 (CTO)Lemanov et al. (1999) which are also called ”quantum paraelectrics”. Microwave measurements of these materials by Rupprecht and BellRupprecht and Bell (1964) demonstrated that there is no sign of appreciable dielectric dispersion up to the GHz range. Therefore it is considered that dielectric properties of pure incipient ferroelectrics are determined by the contribution of the lowest polar-phonon mode (the soft-mode) whose frequency is decreasing on cooling roughly according to the Cochran law. Recently, combining the far infrared (FIR) spectroscopy with low-frequency (LF) dielectric measurements some of present authors confirmed the absence of appreciable dielectric dispersion below the soft-mode response in STO ceramicsPetzelt et al. (2001) and CTO single crystalsŽelezny et al. (2002).

Chen Ang, Bhalla and Cross in their paperChen Ang et al. (2001) investigated the dielectric response of KTO, CTO and (L$`n_{1/2}`$Na1/2)TiO3 (Ln = La and Nd) (LNTO) under a strong dc electric field (up to 60 kV/cm). To describe the field dependence of the permittivity, they attempted to use an equation which is a simple polynomial expansion by even powers of electric field up to the fourth order. They were able to describe the field dependence in CTO and LNTO, but failed to perform a satisfactory fit for KTO. This fact brought the authors to conclusion that there should be another Langevin-type term describing the relaxational contribution of polar clusters to the total polarization. The final equation was of the formChen Ang et al. (2001)

$$ϵ(E)=ϵ_1ϵ_2E^2+ϵ_3E^4+[P_rx/ϵ(0)][\mathrm{cosh}(Ex)]^2$$
(1)
where $`x=P_rL^3E/2k_BT`$, $`P_r`$ is the effective polarization of the polar clusters with the size L and $`k_B`$ is the Boltzmann’s constant. Using Eq. 1 to fit the measured field dependence of the permittivity, they obtained the following parameters of polar clusters: $`P_r=0.4`$ $`\mu `$C/cm2, L = 4.7 – 6 nm and $`\mathrm{\Delta }ϵ_{cluster}`$ = 656 (at 14 K). This would mean that 23% of the total dielectric response at low frequencies are due to some relaxational dispersion of polar clusters below the soft-mode response.

In the following we will show that our recent FIR and LF dielectric measurements of KTO exclude the possibility of any appreciable contribution to the LF permittivity besides the contribution of the soft phonon mode. Moreover we will demonstrate that a correct treatment using the Landau-Ginzburg-Devonshire thermodynamic theory allows to describe the field dependence of the dielectric constant in KTO without using additional Langevin-type terms.

Our measurements have been performed on nominally pure KTO single crystalsTrepakov et al. (2000). The samples were prepared in the form of thin $`<001>`$-oriented plane-parallel plates with dimensions of $`7\times 7\times 0.2`$ mm3. LF dielectric measurements were carried out using a HP 4192A impedance analyzer in the frequency range 100 Hz – 1 MHz with gold electrodes deposited on the main surfaces of the sample. The FIR reflectivity was measured with Bruker IFS 113v Fourier Transform spectrometer in the range 30 – 650 cm-1. The complex dielectric function in the frequency range 5 – 30 cm-1 was obtained by means of time-domain THz transmission spectroscopy on a 65 µm thick plane-parallel sample. The frequencies of the three polar phonon modes and their contribution to the dielectric function were obtained from the fit to the FIR reflectivity using the factorized form of the model dielectric function combined with the directly measured values of the THz complex permittivityKamba et al. (2003) . The temperature dependence of the soft-mode frequency was found to be in perfect agreement with the hyper-Raman dataVogt (1995).

In Fig. 1 the measured permittivity at 43.5 kHz is compared with the low-frequency limit of our fitted FIR dielectric function. The permittivity data from Ref. Chen Ang et al., 2001 are also added for comparison. One can observe good agreement between LF permittivities of KTO crystals and between LF and FIR permittivities in the whole measured temperature range. Somewhat higher value of FIR permittivity compared to the LF values at 20 K lies within the instrumental errors. Anyway, Fig. 1 clearly shows that at all measured temperatures there is no space for any appreciable dispersion below the soft-mode response down to kHz range.

To demonstrate the fitting of the field dependence of permittivity, following Fleury and WorlockFleury and Worlock (1968) let us expand the crystal free energy per unit volume in a power series of lattice polarization:

$`G=G_0+{\displaystyle \frac{1}{2}}\chi (T)(P_x^2+P_y^2+P_z^2)+{\displaystyle \frac{1}{4}}\xi (P_x^4+P_y^4+P_z^4)`$
$`+{\displaystyle \frac{1}{2}}\xi ^{^{}}(P_x^2P_y^2+P_y^2P_z^2+P_x^2P_z^2)+{\displaystyle \frac{1}{6}}\zeta (P_x^6+P_y^6+P_z^6)`$
$`+\mathrm{other}\mathrm{sixth}\mathrm{}\mathrm{order}\mathrm{terms}+\mathrm{}`$ (2)
In our case polarization has only one component directed along the z-axis ($`P_x=P_y=0`$). Thus we can obtain relation between the electric field and polarization:

$$E_z=(G/P_z)_T=\chi (T)P_z+\xi P_z^3+\zeta P_z^5+\mathrm{}$$
(3)
which is a special case of the general relation for arbitrary polarizationFleury and Worlock (1968); Lawless (1978). Then one can derive the expression for the inverse susceptibility:

$$\left(\frac{1}{ϵ_0[ϵ(T,E)1]}\right)_{zz}=E_z/P_z=\chi (T)+3\xi P_z^2+5\zeta P_z^4+\mathrm{}$$
(4)
Here $`\chi (T)=1/ϵ_0[ϵ(T,0)1]`$ is the inverse susceptibility in the zero bias field and $`ϵ_0=8.8510^{12}`$ F/m is the permittivity of vacuum. To calculate the permittivity under dc bias field, one has to use Eq. 3 to calculate the polarization and then insert it into Eq. 4. The calculations should be performed numerically because Eq. 3 is a nonlinear polynomial equation of the 5-th order. Therefore the analytical expression given in Ref. Chen Ang et al., 2001 cannot be valid for highly anharmonic crystals and can be applied only in the weak field approximation. This approximation apparently holds in the case of CTO and LNTO because of their weak nonlinearity for bias fields $`E<60`$ kV/cm. However, in the case of KTO the lattice anharmonicity is essential. The fit using Eq. 3 and 4 of KTO permittivity under the dc electric field at 14 K taken from Ref. Chen Ang et al., 2001 (the first run of measurements) is shown in Fig. 2. The fit uses two independent parameters $`\xi =10^{10}`$ and $`\zeta =10^{11}`$ in SI units and the value $`ϵ`$(14 K, 0) = 2821 taken from Ref. Chen Ang et al., 2001. The fitted nonlinear coefficients are in reasonable agreement with the previously published data (see Ref. Fleury and Worlock, 1968 and references therein). Fig. 2 demonstrates that the expansion up to the 6-th order in Eq. Comment on ”Dielectric behavior of paraelectric KTaO3, CaTiO3, and (Ln1/2Na1/2)TiO3 under a dc electric field” is sufficient to ensure a good agreement between the theory and experiment.

The critical limit for the weak field approximation can be found from the condition that 4-th and 6-th order terms in Eq. Comment on ”Dielectric behavior of paraelectric KTaO3, CaTiO3, and (Ln1/2Na1/2)TiO3 under a dc electric field” are small compared to the quadratic term

$$\xi P_{cr}^2+\zeta P_{cr}^4=\lambda \chi ,$$
(5)
where $`\lambda 1`$ is a dimensionless constant. For the typical values of nonlinear coefficientsFleury and Worlock (1968) we can neglect the 4-th order term in Eq. 5 and using Eq. 3 we obtain

$$E_{cr}=\chi \sqrt{\frac{\lambda \chi }{\xi }}.$$
(6)
Using our fit parameters for KTO and $`\lambda =0.1`$ we calculate $`E_{cr}1.3`$ kV/cm. One can see that this value is much smaller than the highest dc electric field in Ref. Chen Ang et al., 2001 (60 kV/cm).

It is believed that polar nano-clusters are induced by unavoidable dipolar impurities in nominally pure KTOVogt (1991). The presence of such clusters appears in the LF dielectric spectra of KTO as a universal low-temperature dielectric dispersionTrepakov et al. (2000) with thermally activated behavior. Below 30 K this dielectric dispersion shifts below the frequency range of conventional LF dielectric spectroscopy. Dielectric spectroscopy can only establish that the contribution to permittivity of polar clusters is very small comparing to the phonon contribution above 30 K and cannot be observed in LF dielectric spectra below 30 K. On the other hand, the polar clusters can be detected by light scattering as it has been done in RamanUwe et al. (1986) and hyper-RamanVogt (1991) scattering experiments. These experiments, however, bring no information about the dynamics of the clusters because their relaxation is quasi-static comparing to the frequencies resolved in the light scattering experiments. The local lattice symmetry breaking due to the polar clusters relaxes the selection rules and induces activation of new lines. These changes are in agreement with existence of the polar clusters and do not contradict with our dielectric spectroscopy results.

Taking into account all the mentioned arguments, we can conclude that the estimated size of polar clusters found in Ref. Chen Ang et al., 2001 does not have any reliable physical background and its possible correspondence with the Raman resultsUwe et al. (1986) is accidental.

Analysis of the field dependence of permittivity using equation similar to Eq. 1 was previously applied also to $`\mathrm{Sr}_{1\mathrm{x}}\mathrm{Ca}_\mathrm{x}\mathrm{TiO}_3`$ systemDec et al. (1995); Bianchi et al. (1995) where, due to the Ca-doping, appreciable dielectric dispersion indeed exists. However, the applied electric fields were much smaller ($`E<1`$ kV/cm) than those in Ref. Chen Ang et al., 2001. Estimation of the critical electric field from Eq. 6 gives for STO $`E_{cr}0.2`$ kV/cm (nonlinear parameters are taken form Ref. Fleury and Worlock, 1968). Thus even 1 kV/cm field can introduce some error if one uses Eq. 1. In our opinion, the correct thermodynamic approach described here is always preferable for the analysis of nonlinear behavior under strong electric fields.

This work was supported by the Grant Agency of the Czech Republic (projects No. 202/01/0612 and 202/02/0238) and Czech Academy of Sciences (project AVOZ-010-914).


text:
# The pants complex has only one end

## 1. Definitions and statement of the main theorem

The purpose of this note is to prove the following theorem:

Theorem 4.1. Let $`S`$ be a closed orientable surface with genus $`g(S)3`$. Then the pants complex of $`S`$ has only one end. In fact, there are constants $`K=K(S)`$ and $`M_3=M_3(S)`$ so that, if $`R>M_3`$, any pair of pants decompositions $`P`$ and $`Q`$, at distance greater than $`KR`$ from a basepoint, can be connected by a path which remains at least distance $`R`$ from the basepoint.

Recall that a pants decomposition of $`S`$ consists of $`3g(S)3`$ disjoint essential non-parallel simple closed curves on $`S`$. Each component of the complement of the curves is a sphere with $`3`$ holes or pair of pants. Then the pants complex $`C_P(S)`$ is the metric graph whose vertices are pants decompositions of $`S`$, up to isotopy. Two vertices $`P,P^{}`$ are connected by an edge if $`P,P^{}`$ differ by an elementary move. In an elementary move all curves in the pants are fixed except for one curve $`\alpha `$. Removing $`\alpha `$, the component of the complement of the remaining curves that contains $`\alpha `$ is either a punctured torus or a sphere with $`4`$ holes. Then $`\alpha `$ is replaced by a curve $`\beta `$ contained in this domain that intersects $`\alpha `$ minimally; in the punctured torus case once, and in the sphere case twice. All edges of $`C_P(S)`$ are assigned length $`1`$. We let $`d(,)`$ be the distance function in $`C_P(S)`$. The pants complex $`C_P(S)`$ is known to be connected \[HT\].

Recall that a metric space $`(X,d)`$ has one end if for any basepoint $`OX`$ and any radius $`R`$ the complement of $`B_R=B_R(O)`$, the ball of radius $`R`$ centered at $`O`$, has only one unbounded component. It is easy to see that the definition does not depend on the choice of the point $`O`$. Clearly having one end is a quasi-isometry invariant. So, following Brock \[B\], our theorem implies:

###### Corollary 1.1.

Teichmüller space, equipped with the Weil-Peterson metric, has only one end.

Finally, recall that the curve complex $`𝒞(S)`$ is the complex whose $`k`$-simplices consist of $`k+1`$ distinct isotopy classes of essential simple closed curves on $`S`$ that have disjoint representatives on $`S`$. Or, in the case of a once-punctured torus and four-punctured sphere, $`𝒞(S)`$ is the Farey graph. From the metric point of view we will only interested in the $`1`$-skeleton of $`𝒞(S)`$. Each edge is assigned length $`1`$. We let $`d_S(,)`$ denote the distance function in $`𝒞(S)`$.

## 2. The set of handle curves is connected

In this section we prove two combinatorial facts. First, the set of handle curves in the curve complex is connected and second, any pants decomposition is a bounded distance (in the pants complex) from a decomposition containing a handle curve.

Again assume $`S`$ is a closed orientable surface with genus three or greater. We will call $`\alpha `$ a handle curve in $`S`$ if $`\alpha `$ separates $`S`$ into two surfaces: the handle $`S(\alpha )𝕋^2\{\mathrm{pt}\}`$ and the rest of the surface.

We will need the following result. It was first proved by Farb and Ivanov \[FI\] by different methods. Another proof was given by McCarthy and Vautaw \[MV\] by methods similar to ours. We include a proof for completeness.

###### Proposition 2.1.

If $`g(s)3`$, the subcomplex $`h(S)𝒞(S)`$ of handle curves is connected.

###### Remark 2.2.

Note that the hypothesis $`g(S)3`$ cannot be removed; it is easy to check that the handle curves for a closed surface of genus $`2`$ do not form a connected set.

###### Remark 2.3.

Note that Proposition 2.1 immediately implies that the set of separating curves in $`𝒞(S)`$ is also connected.

###### Remark 2.4.

We note that Proposition 2.1 may be generalized to the case $`S\mathrm{}`$. A still open question is the higher connectivity of $`h(S)`$.

Before we begin the proof we will require a bit of terminology. Let $`i(,)`$ denote the geometric intersection number of two essential simple closed curves in $`S`$. Also, if $`\delta `$ is a dividing curve in $`S`$ we say that an arc $`\beta ^{}`$ is a wave for $`\delta `$ if $`\beta ^{}\delta =\beta ^{}`$ and $`\beta ^{}`$ is essential as a properly embedded arc in $`S\delta `$. We say that two waves $`\beta ^{}`$ and $`\beta ^{\prime \prime }`$ for $`\delta `$ link if $`\beta ^{}\beta ^{\prime \prime }=\mathrm{}`$, both $`\beta ^{}`$ and $`\beta ^{\prime \prime }`$ meet the same side of $`\delta `$, and $`\beta ^{}`$ separates $`\beta ^{\prime \prime }`$ inside $`\delta `$.

Finally we define double surgery as follows. Suppose we are given a linking pair of waves $`\beta ^{}`$ and $`\beta ^{\prime \prime }`$ for an essential dividing curve $`\delta _0`$. Form the closed regular neighborhood $`U=\mathrm{neigh}(\delta _0\beta ^{}\beta ^{\prime \prime })`$. Let $`\delta _1`$ be the component of $`U`$ which is not homotopic in $`S`$ to $`\delta _0`$. We say that $`\delta _1`$ is obtained from $`\delta _0`$ via double surgery along $`\beta ^{}`$ and $`\beta ^{\prime \prime }`$. Note that $`\delta _1`$ is necessarily a dividing curve and is disjoint from $`\delta _0`$. If the component of $`S\delta _0`$ containing $`\beta ^{}\beta ^{\prime \prime }`$ is not a handle then $`\delta _1`$ is also essential.

We are now equipped to prove the proposition:

###### Proof of Proposition 2.1.

Let $`\alpha ,\beta h(S)`$ be handle curves and $`S`$ a closed orientable surface of genus at least three. Suppose that $`\alpha `$ and $`\beta `$ are tight: $`\alpha `$ has been isotoped to make $`|\alpha \beta |=i(\alpha ,\beta )`$. If $`i(\alpha ,\beta )=0`$ then there is nothing to prove. If $`i(\alpha ,\beta )>0`$ we will find a curve $`\gamma h(S)`$ with $`i(\gamma ,\alpha )=0`$ and $`i(\gamma ,\beta )<i(\alpha ,\beta )`$. By induction, $`\gamma `$ will be connected to $`\beta `$ in $`h(S)`$, proving the proposition.

We find $`\gamma `$ via the following inductive procedure. Recall that $`S(\alpha )`$ is the handle which $`\alpha `$ bounds. To begin, we define $`\delta _0SS(\alpha )`$ to be a parallel copy of $`\alpha `$, still intersecting $`\beta `$ tightly. At stage $`k`$ by induction we will be given an essential dividing curve $`\delta _k`$ where

* $`i(\alpha ,\delta _k)=0`$,
* $`\delta _k`$ is tight with respect to $`\beta `$, and
* $`i(\delta _k,\beta )<i(\delta _{k1},\beta )`$, if $`k>0`$.

Let $`T_k`$ be the component of $`S\delta _k`$ which does not contain $`\alpha `$. If $`T_k`$ is a handle, then we take $`\gamma =\delta _k`$ and we are done with the inductive procedure. If $`i(\delta _k,\beta )=0`$ then we may take $`\gamma `$ to be any handle curve inside $`T_k`$. As this $`\gamma `$ satisfies $`i(\alpha ,\gamma )=i(\beta ,\gamma )=0`$ this would finish the proposition. From now on we assume that $`T_k`$ is not a handle and that $`i(\delta _k,\beta )>0`$.

We now attempt to do a double surgery of $`\delta _k`$ into $`T_k`$. Either we will find $`\gamma `$ directly or the curve resulting from double surgery, $`\delta _{k+1}`$, will satisfy the induction hypothesis. As the geometric intersection with $`\beta `$ is always decreasing, this procedure will stop after finitely many steps yielding the desired handle curve.

So all that remains is to do the double surgery. Recall that we are given $`\alpha ,\beta `$ tightly intersecting handle curves and we are also given $`\delta _k`$ satisfying the induction hypotheses. Recall also that $`T_k`$ is the component of $`S\delta _k`$ which does not contain $`\alpha `$. Recall $`T_k`$ is not a handle and that $`i(\delta _k,\beta )>0`$.

Suppose further that $`\beta ^{},\beta ^{\prime \prime }\beta T_k`$ are linking waves for $`\delta _k`$. As described above we may form $`\delta _{k+1}`$ via a double surgery along $`\beta ^{}`$ and $`\beta ^{\prime \prime }`$. Isotope $`\delta _{k+1}`$, in the complement of $`\delta _k`$, to be tight with respect to $`\beta `$. As noted in the definition of double surgery, $`\delta _{k+1}`$ is an essential dividing curve which is disjoint from $`\alpha `$. Finally note that $`i(\delta _{k+1},\beta )i(\delta _k,\beta )4`$. Thus all of the induction hypotheses are satisfied.

Suppose now that we cannot find linking waves among the arcs $`\beta T_k`$. Choose instead an outermost wave $`\beta ^{}\beta T_k`$: that is, there exists an arc $`\delta _k^{}\delta _k`$ such that $`\delta _k^{}\beta =\delta _k^{}=\beta ^{}`$.

Here there are two remaining cases. If $`\delta _k^{}\beta ^{}`$ is a separating curve take $`\delta _{k+1}=\delta _k^{}\beta ^{}`$ and note that the induction hypotheses are easily verified. The final possibility is that $`\delta _k^{}\beta ^{}`$ is not separating. In this case choose a properly embedded essential arc $`\beta ^{\prime \prime }T_k`$ such that $`\beta ^{\prime \prime }\beta =\mathrm{}`$ and $`|\beta ^{\prime \prime }\delta _k^{}|=1`$. Then $`\beta ^{}`$ and $`\beta ^{\prime \prime }`$ link. Do a double surgery along these waves to obtain $`\delta _{k+1}`$. Isotope $`\delta _{k+1}`$, in the complement of $`\delta _k`$, to be tight with respect to $`\beta `$. Again, all of the induction hypotheses are easily verified, as we have $`i(\delta _{k+1},\beta )i(\delta _k,\beta )2`$. This completes the second induction step and hence completes the proof of Proposition 2.1. ∎

We also require

###### Lemma 2.5.

There is a constant $`M_3=M_3(S)`$ such that the pants decompositions containing a handle curve are $`M_3`$-dense in the space of all pants decompositions

###### Proof.

The mapping class group acts co-compactly on the space of pants decompositions. ∎

## 3. Subsurface projections and distances

Here we give two lemmas studying the pants complex. The first gives a condition for a pants decomposition to lie outside of a large ball about the origin in $`C_P(S)`$ while the second provides us with useful paths laying outside of such a ball. This uses techniques developed by Masur and Minsky \[MM\].

Given a subsurface $`WS`$ and a curve $`\gamma `$ that intersects $`W`$, we may define a projection of $`\gamma `$ to $`𝒞(W)`$ which associates to $`\gamma `$ a collection of curves in $`𝒞(W)`$. Namely, the intersections of $`\gamma `$ with $`W`$ fall into finitely many homotopy classes of disjoint arcs (and curves) relative to $`W`$. Let $`\alpha `$ be any such arc or curve. Let $`U`$ be a regular neighborhood of $`\alpha W`$. Then $`UW`$ is a curve in $`𝒞(W)`$ if $`\alpha `$ is a curve or is an arc connecting distinct components of $`W`$. Otherwise $`UW`$ is a pair of curves in $`𝒞(W)`$. Clearly this curve or pair of curves in $`𝒞(W)`$ depend only on the homotopy class of $`\alpha `$. If $`\alpha ,\beta `$ are disjoint homotopy classes of arcs, then any pair of curves $`\alpha ^{},\beta ^{}`$ built out of this surgery satisfy $`d_W(\alpha ^{},\beta ^{})2`$ (\[MM\] Lemma 2.3). Thus we can define $`\pi _W(\gamma )`$ to be the corresponding subset of diameter at most $`2`$ in $`𝒞(W)`$.

Similarly, given a pants decomposition $`P`$ we may project each curve of $`P`$ that intersects $`W`$ into $`W`$. We denote the resulting image set in $`𝒞(W)`$, which has diameter at most $`4`$, by $`\pi _W(P)`$. By $`d_W(P,P^{})`$ we mean the distance in the curve complex of $`W`$ between the sets $`\pi _W(P)`$ and $`\pi _W(P^{})`$.

Let $`[x]_C`$ be the function on $``$ giving zero if $`x<C`$ and giving $`x`$ if $`xC`$. We will need the following result from \[MM\] (see Theorem 6.12 and Section 8 of that paper):

There is a constant $`C_0=C_0(S)1`$ such that for any $`CC_0`$ there are constants $`M_1=M_1(C)C`$ and $`M_2=M_2(C)0`$ with the following property: for any pants decompositions $`P,P^{}`$ we have

$$\frac{1}{M_1}\underset{V}{}[d_V(P,P^{})]_CM_2d(P,P^{})M_1\underset{V}{}[d_V(P,P^{})]_C+M_2,$$
(1)
where the sums range over subsurfaces $`VS`$ with essential boundary and where $`V`$ is not an annulus nor is it a thrice-punctured sphere.

Fix now such a $`C>2`$. It follows from equation (1) that the projections that are at the critical values $`C,C+1`$ cannot account for the entire distance. Namely there are constants $`c=c(C)>1`$ and $`M_4=M_4(C)>0`$ such that

$$\underset{V}{}[d_V(P,P^{})]_Cc\underset{V}{}[d_V(P,P^{})]_{C+2}+M_4.$$
(2)
Choose $`K=K(C)>0`$ so that for all $`R1`$,

$$\frac{1}{2cM_1}\left((K1)RM_2M_1M_4cM_1^2(R+M_2)\right)>M_1(R+M_2)$$
(3)
Also, choose a basepoint $`OC_P(S)`$ and let $`B_R=B_R(O)`$ be the ball of radius $`R`$ centered at $`O`$.

###### Lemma 3.1.

Fix a handle curve $`\alpha `$ and some curve $`\alpha ^{\prime \prime }S(\alpha )`$ satisfying $`d_{S(\alpha )}(O,\alpha ^{\prime \prime })>M_1(R+M_2)`$. For any pants decomposition $`P`$ containing $`\alpha ^{\prime \prime }`$ we have $`PB_R`$.

###### Proof.

Note that $`d_{S(\alpha )}(O,\alpha ^{\prime \prime })C`$. So, by the left inequality of (1) any pants $`P`$ containing $`\alpha ^{\prime \prime }`$ has

$$d(P,O)\frac{1}{M_1}[d_{S(\alpha )}(P,O)]_CM_2\frac{1}{M_1}d_{S(\alpha )}(\alpha ^{\prime \prime },O)M_2>R.$$

As the Farey graph for $`S(\alpha )`$ has infinite diameter, and as the diameter of $`\pi _{S(\alpha )}(O)`$ is bounded, such curves $`\alpha ^{\prime \prime }`$ exist in abundance. We now turn to the existence of paths lying outside of the $`R`$-ball about the basepoint.

###### Lemma 3.2.

Suppose $`P_0`$ is a pants decomposition of $`S`$ such that $`P_0B_{(K1)R}`$ and $`P_0`$ contains a curve $`\alpha `$ which bounds a handle $`S(\alpha )`$. Then there is a path $`P_t`$ starting at $`P_0`$ such that

* for all $`t`$, $`P_t|(SS(\alpha ))=P_0|(SS(\alpha ))`$.
* for all $`t`$, $`P_tB_R`$
* The endpoint of the path, $`P_1`$, contains a curve $`\alpha ^{\prime \prime }S(\alpha )`$ which does not appear in any pants decomposition in $`B_R`$.

###### Proof.

Let $`\alpha ^{}P_0`$ be the curve strictly contained in $`S(\alpha )`$. Consider a geodesic segment in the Farey graph connecting $`\alpha ^{}`$ to $`\beta \pi _{S(\alpha )}(O)`$, where $`\beta `$ is chosen as close as possible to $`\alpha ^{}`$. Extend this segment through $`\alpha ^{}`$ to a geodesic ray $`L`$ in the direction opposite $`\beta `$. The ray $`L`$ meets the segment only at $`\alpha ^{}`$. Move along $`L`$ distance more than $`M_1(R+M_2)`$ from $`\alpha ^{}`$ to a point $`\alpha ^{\prime \prime }`$. Let $`P_t`$ be the path obtained by making elementary moves along the curves in $`L`$ and fixing the pants in $`SS(\alpha )`$.

Suppose first that $`d_{S(\alpha )}(\beta ,\alpha ^{})>M_1(R+M_2)`$. Then by Lemma 3.1 any pants $`\widehat{P}`$ containing any $`\alpha _tL`$ has $`d(\widehat{P},O)>R`$. So Lemma 3.2 holds in this case.

Next suppose that $`d_{S(\alpha )}(\beta ,\alpha ^{})M_1(R+M_2)`$. Then by (1) and (2)

$$(K1)Rd(P_0,O)M_1\underset{V}{}[d_V(P_0,O)]_C+M_2$$

$$cM_1\underset{V}{}[d_V(P_0,0)]_{C+2}+M_1M_4+M_2$$

$$cM_1\underset{VS(\alpha )}{}[d_V(P_0,O)]_{C+2}+cM_1^2(R+M_2)+M_1M_4+M_2.$$

Let $`V`$ be any subsurface disjoint from $`S(\alpha )`$. Since $`P_t`$ is constant in $`V`$, the projection $`\pi _V(P_t)`$ is constant. Now let $`V`$ be a subsurface that intersects $`S(\alpha )`$ or strictly contains $`S(\alpha )`$. Since $`\alpha P_t`$, it follows that $`\pi _V(P_t)`$ contains $`\pi _V(\alpha )`$. Since each $`\pi _V(P_t)`$ has diameter at most $`2`$, $`d_V(P_t,O)d_V(P_0,O)2`$. Thus for any subsurface $`V`$ not isotopic to $`S(\alpha )`$, as $`C>2`$, we have $`[d_V(P_t,O)]_C\frac{1}{2}[d_V(P_0,O)]_{C+2}`$. Thus,

$$\underset{VS(\alpha )}{}[d_V(P_t,O)]_C\frac{1}{2}\underset{VS(\alpha )}{}[d_V(P_0,0)]_{C+2}$$

$$\frac{1}{2cM_1}\left((K1)RM_2M_1M_4cM_1^2(R+M_2)\right)>M_1(R+M_2),$$

the last inequality following from (3). So, by (1),

$$d(P_t,O)>R.$$

Finally, as the pants decomposition $`P_1`$ contains $`\alpha ^{\prime \prime }`$ and

$$d_{S(\alpha )}(O,\alpha ^{\prime \prime })d_{S(\alpha )}(\alpha ^{},\alpha ^{\prime \prime })>M_1(R+M_2)$$

we have $`P_1B_R`$, by Lemma 3.1. ∎

## 4. Proof of the theorem

Recall the statement:

###### Theorem 4.1.

Let $`S`$ be a closed orientable surface with genus $`g(S)3`$. Then the pants complex of $`S`$ has only one end. In fact, there are constants $`K=K(S)`$ and $`M_3=M_3(S)`$ so that, if $`R>M_3`$, any pair of pants decompositions $`P`$ and $`Q`$, at distance greater than $`KR`$ from a basepoint, can be connected by a path which remains at least distance $`R`$ from the basepoint.

###### Proof.

We take $`M_3`$ as defined in Section 2 and $`K`$ as defined in Section 3.

First move $`P`$ and $`Q`$ a distance at most $`M_3<R`$ to obtain pants decompositions $`P_0`$ and $`Q_0`$ which contain handle curves $`\alpha _PP_0`$, $`\alpha _QQ_0`$ and such that $`P_0,Q_0B_{(K1)R}`$.

Apply Lemma 3.2 twice in order to connect $`P_0`$ and $`Q_0`$ to pants decompositions $`P_1`$ and $`Q_1`$ satisfying all of the conclusions of the lemma. Let $`\alpha _P^{\prime \prime }`$ and $`\alpha _Q^{\prime \prime }`$ be the curves lying in the handles $`S(\alpha _P)P_1`$ and $`S(\alpha _Q)Q_1`$ respectively.

We must now construct a path from $`P_1`$ to $`Q_1`$. Consider first the case where $`\alpha _P\alpha _Q`$.

Applying Proposition 2.1 we connect $`\alpha _PP_1`$ and $`\alpha _QQ_1`$ by a path $`\{\alpha _i\}_{i=1}^n`$ of handle curves in $`𝒞(S)`$, the curve complex of $`S`$. Here we have $`\alpha _1=\alpha _P`$, $`\alpha _n=\alpha _Q`$, and $`n>1`$. Note that in this step the hypothesis $`g(S)>2`$ is used. Choose, for $`i\{2,3,\mathrm{}n1\}`$, any curve $`\alpha _i^{\prime \prime }S(\alpha _i)`$ such that $`d_S(\alpha _i)(O,\alpha _i^{\prime \prime })>M_1(R+M_2)`$. Set $`\alpha _1^{\prime \prime }=\alpha _P^{\prime \prime }`$ and $`\alpha _n^{\prime \prime }=\alpha _Q^{\prime \prime }`$. Let $`P_n=Q_1`$.

Inductively, we connect $`P_i`$ by a path to $`P_{i+1}`$ where, first, every pants decomposition in the path contains $`\alpha _i`$ and $`\alpha _i^{\prime \prime }`$ and, second, $`P_{i+1}`$ also contains $`\alpha _{i+1}`$ and $`\alpha _{i+1}^{\prime \prime }`$. (This is possible because $`C_P(SS(\alpha _i))`$ is connected.) By Lemma 3.1 this path lies outside of the ball of radius $`R`$ and we are done.

In the case which remains we have $`\alpha _P=\alpha _Q`$. So there is no need for Proposition 2.1. Instead we choose any handle curve $`\beta `$ which is disjoint from $`\alpha _P`$. Note that $`\beta `$ exists as $`g(S)>2`$. Choose also any $`\beta ^{\prime \prime }`$ satisfying the hypothesis of Lemma 3.1. We now consider the sequence $`\alpha _P,\beta ,\alpha _Q`$ as a path of length two in the curve complex and connect $`P_1`$ to $`Q_1`$ as in the previous paragraph. This completes the proof. ∎


text:
# Ubiquity of metastable-to-stable crossover in weakly chaotic dynamical systems

## 1 Introduction

The statistical-mechanical comprehension of the occurrence and eventual disappearance of out-of-equilibrium quasistationary states (QSS), referred here also as metastable states, is an important undertaking that nowadays attracts the attention of researchers in various fields, ranging from weakly chaotic nonlinear dynamics to slow glassy dynamics in condensed matter. Under specific classes of initial conditions, the dynamics in these systems displays two-step relaxation processes, entering first one or more periods of little variation before a crossover to the final equilibrium state. The QSS are conjectured to emerge as the result of long-ranged correlations or interactions amongst the elements of the system or long-ranged time correlations in the orbits of iterated maps. An interesting statistical-mechanical model that displays QSS is the Hamiltonian Mean Field (HMF) model , that in its simplest version consists of a set of $`N`$ inertial XY classical spins or rotors all of which interact equally with each other. But also, interestingly, the most relevant features of QSS dynamics are displayed by considerably simpler nonlinear dynamical maps, as it is the case of conservative coupled maps  and also by single, dissipative, one-dimensional maps, like the prototypical logistic map close to the edge of chaos and in the presence of external noise .

In this paper we present a comparative study of the aforementioned models, starting our discussion with low-dimensional nonlinear maps and ending it with the HMF model. We emphasize the analogies found between their dynamical properties and indicate the mechanisms that give raise to the QSS as well as those responsible for the crossover to the equilibrium state. In all of these cases the quantities that describe dynamical evolution seem to be associated to the predictions of the nonextensive generalization  of the Boltzmann-Gibbs (BG) statistical mechanics. Also, in all of these examples there is clear evidence of the non-commutability of the infinite time limit with other important parameter limits, such as the thermodynamic limit ($`N\mathrm{}`$) or the limit of vanishing noise amplitude. Moreover, there is a common indication that, in the characterization of the phase space available for dynamical evolution, there appears evidence of fractalization and of ergodicity failure.

The paper is organized as follows. In the next section we describe how the dynamics at the edge of chaos of the logistic map leads, via the addition of external noise, to the crossover from QSSs to a strongly chaotic state that is the analog of the BG equilibrium state . We present here new numerical evidence by means of an ensemble implementation of the dynamics in this system. In section 3 we note that the emergence of QSSs and eventual departure from them appear associated to the complex structure of phase space accessibility in symplectic maps governed by the Kolmogorov-Arnold-Moser (KAM) theorem . We also show new results that indicate that the QSS, present at the lowest possible dimensions, persists as hundreds of symplectic maps are (globally) coupled. Section 4 is devoted to the QSS of the HMF model. Particular attention is paid to the time evolution of the anomalous probability distribution functions (PDFs) for the velocities, and also to basic questions referring to the measurability of temperature and to the zeroth principle of thermodynamics as it relates to the QSS . Finally, in section 5 we draw our conclusions, we make reference to the existing analogies amongst the different models described and discuss the underlying mechanism that gives raise to the QSS.

## 2 QSS in the logistic map with external noise

Let us start our analysis by recalling  that phenomena like two-step relaxation and aging can be displayed at the lowest possible dimension, namely, by means of a one-dimensional map. For this purpose we consider the paradigmatic logistic map with additive white noise,

$$x_{\tau +1}=1\mu x_\tau ^2+\xi _\tau \sigma ,\tau =0,1,\mathrm{},x[1,1],\mu [0,2],$$
(1)
where $`\xi _\tau `$ is Gaussian-distributed and delta-correlated, $`\xi _\tau _\xi =0`$, $`\xi _\tau \xi _\tau ^{}_\xi =\delta _{\tau \tau ^{}}`$, and $`\sigma `$ measures the noise intensity. As it is well known, for $`\sigma =0`$ the Feigenbaum attractor at the edge of chaos $`\mu =\mu _c1.401155189\mathrm{}`$ is the accumulation point of both the period-doubling and band-splitting cascades. The Lyapunov coefficient is zero  and the attractor is a cantor set of fractal dimension $`d_f=0.5338\mathrm{}`$. It has been analytically proved by means of the conventional renormalization-group approach for this type of map that the dynamics at this critical point is unmistakably associated to the nonextensive statistical mechanics . Specifically, the sensitivity to initial conditions is actually given in a closed form by a $`q`$-exponential function and it is related to the rate of entropy production via a $`q`$-generalized Pesin identity, linking the sensitivity to initial conditions to the nonextensive expression for the entropy, $`S_q=\mathrm{ln}_qW`$ (where $`W`$ is the number of equally-probable positions of an ensemble of iterates at a given time). For this stationary state the index $`q`$ takes the specific value $`q=0.2445\mathrm{}`$. Let us remind the reader that the $`q`$-exponential and its inverse the $`q`$-logarithm are defined, respectively, as $`\mathrm{exp}_q(x)[1+(1q)x]^{1/1q}`$ and $`\mathrm{ln}_q(x)(x^{1q}1)/(1q)`$, and that these expressions reduce to the usual exponential and logarithmic functions in the limit $`q1`$.

More precisely, the Feigenbaum attractor can be described by means of a set of position subsequences generated when the iterative map is given the initial position $`x_0=0`$ . See the empty circles in Fig. 1(a) where we have plotted the absolute value of these positions $`\left|x_\tau \right|`$ and time $`\tau `$ at which they are reached in logarithmic scales. Each of the subsequences is labelled by an integer $`k=0,1,\mathrm{}`$ and each is generated by the time subsequence $`\tau _k=(2k+1)2^{nk}`$, with $`n=0,1,\mathrm{}`$ and $`nk`$. This can be appreciated in Fig. 1(a). In the following we will focus on the principal subsequence $`\tau _0=2^n`$. Another important feature of Fig. 1(a) is that the logarithm of the attractor positions appear clearly divided into an infinite number of equally-spaced horizontal bands. In the empty space between any two bands there are points called repellers, the unstable solutions of the iterated equation $`x=1\mu _cx^2`$. (See the small dots in Fig. 1(a)). At the chaos threshold of the logistic map with $`\sigma 0`$ the dynamics is governed by two competing phenomena: On the one hand the action of the noise is to drive the iterates out of the bands, and on the other hand the instability in the neighborhood of the repellers sends the positions back into the attractor bands. When the noise intensity is large enough to fill a gap between bands, the result is that the iterates are able to jump from one band to the next. (See full and dashed lines in Fig. 1(a)). Since the band gaps in $`\left|x\right|`$-space are of variable size, increasing in size as $`\left|x\right|1`$, the motion becomes confined among some few bands for a certain amount of time but at longer times the iterates are able to spread to other more distant bands.

The analysis of the crossover time $`t_c`$ for the broadening of the attractor bands produced by the noise can be made quantitative. At the edge of chaos $`\mu =\mu _c`$, for $`\sigma 0`$, the iterated position of $`x_0=0`$ at shifted time $`t=\tau _01`$ is given by (see  for details):

$$x_t=\mathrm{exp}_{2q}(\lambda _qt)[1+\xi _t\sigma \mathrm{exp}_r(\lambda _rt)],$$
(2)
where $`q=1\mathrm{ln}2/\mathrm{ln}\alpha =0.2445\mathrm{}`$, $`\lambda _q=\mathrm{ln}\alpha /\mathrm{ln}2=1.3236\mathrm{}`$ is the generalized Lyapunov coefficient, $`r=1\mathrm{ln}2/\mathrm{ln}\kappa 0.633`$, $`\lambda _r=\mathrm{ln}\kappa /\mathrm{ln}22.727`$ ($`\alpha =2.50290\mathrm{}`$ being one of Feigenbaum’s universal constants and $`\kappa 6.619`$ a constant associated to the noise term). Let us consider now a different trajectory $`x_t^{}`$ also starting at the origin at $`t=0`$, $`x_0^{}=0`$. Since the noise acts independently on each trajectory, after time $`t`$ the orbit’s separation is given by

$$x_tx_t^{}=\mathrm{exp}_{2q}(\lambda _qt)\mathrm{exp}_r(\lambda _rt)\sigma (\xi _t\xi _t^{}).$$
(3)
At time $`t=2^n1`$ the separation between bands at which the iterate is located is of order $`\mathrm{exp}_{2q}(\lambda _qt)`$. Using $`\xi _t\xi _t^{}1`$ we obtain an estimate of the crossover time $`t_c`$ at which the noise-assisted jumps into different bands starts to take place. This is

$$\mathrm{exp}_r(\lambda _rt_c)\sigma 1t_c\sigma ^{r1}.$$
(4)
In order to display the many-plateaux structure of time evolution in this dynamical model we consider next an ensemble of $`M=10^3`$ independent copies of the map (1) at the edge of chaos $`\mu =\mu _c`$. The initial positions are $`x_0=0M`$ and we measure the variance $`x_t^2x_t^2`$ (which plays a role analogous to the temperature in the statistical-mechanical description of Hamiltonian systems) of the distribution of positions at time $`t=2^n1`$. Fig. 1(b) shows a numerical corroboration of previous results , for different values of the noise intensity $`\sigma `$. While in the first plateau the variance is zero, as the dynamics of all $`M`$ trajectories is confined to the attractor, at the crossover time $`t_c\sigma ^{r1}`$ the ensemble of orbits starts spreading and the variance increases until a second plateau is reached where it saturates as a consequence of the confinement induced by the unstable repeller regions. As it can be well appreciated in Fig. 2, for long enough times it is possible to see another increment of the variance as the iterates are able to reach other bands of the attractor. The whole process ends when all bands become occupied and the variance is consequently of order $`1`$.

A two-step process of relaxation is known to be one of the main dynamical properties displayed by supercooled liquids close to glass formation . Two other characteristic features of glassy dynamics are the so-called Adam-Gibbs formula and the scaling property known as aging. The first is a relation between the relaxation time $`t_c`$ (e.g. the viscosity or the inverse of the difussivity) and the entropy $`S_{conf}`$ associated to the number of possible molecular minimum energy configurations in the fluid . The second is the loss of time translation invariance, named aging , that is due to the fact that properties of glasses depend on the procedure by which they are obtained. The time decrease of relaxation functions and correlations exhibit a scaling dependence on the ratio $`t/t_w`$ where $`t_w`$ is a waiting time. As the counterpart to the Adam-Gibbs formula, it has recently been shown  that the logistic map at $`\mu _c`$ for $`\sigma 0`$ presents a relationship between the plateau duration $`t_c`$, and the entropy $`S_{conf}`$ for the state that comprises the largest number of (iterate positions) bands allowed by noise intensity. This entropy is obtained from the probability of chaotic band occupancy at position $`x`$ . Also, the trajectories at $`\mu _c`$ in the limit $`\sigma 0`$ have been shown  to obey a scaling property characteristic of aging, of the form $`x_{t+t_w}=g^{(t_w)}(0)\mathrm{exp}_q(\lambda _qt/t_w)`$ where $`g^{(t_w)}(0)`$ is the $`t_w`$-times composition of the Feigenbaum fixed-point map function $`g(0)`$.

## 3 QSS in coupled symplectic maps

Symplectic maps are important tools in the study of Hamiltonian systems. If we consider a time-independent Hamiltonian system with $`n`$ degrees of freedom ($`2n`$-dimensional Gibbs $`\mathrm{\Gamma }`$-space), a $`(2n2)`$-dimensional symplectic map is the result of taking a Poincaré section over the constant energy hypersurface . The recurrence time is discrete and the map is useful in displaying the stationary, recurrent properties of the original Hamiltonian system. Interestingly enough, a $`2n`$-dimensional symplectic map is also the result of a Poincaré section in the phase space of a time-dependent Hamiltonian system with $`n`$ degrees of freedom ($`2n+1`$ dimensions), with a periodic dependence on time . As a consequence of the symplectic structure, the Lyapunov spectra in the $`d`$-dimensional phase space of the map ($`d`$ being an even number) is characterized by $`d/2`$ pairs of Lyapunov coefficients, where each element of the pair is the negative of the other.

A prototypical $`2`$-dimensional symplectic map, that displays the KAM mechanism for the transition from regularity to chaos, is the Chirikov-Taylor standard map (see, e.g.,  and references therein). Since fundamental dynamical processes, like for example Arnold diffusion, occur only for $`d>2`$, it is interesting to consider more generally the case of $`N`$ symplectic and globally coupled standard maps ($`d=2N`$), described, for example, by the equations:

$`\begin{array}{cccc}\hfill \theta _i(t+1)& =& p_i(t+1)+\theta _i(t)+\frac{b}{N1}\underset{\begin{array}{c}m=1\\ mi\end{array}}{\overset{N}{}}p_m(t+1)\hfill & \hfill (\mathrm{mod}\mathrm{\hspace{0.33em}1}),\\ \hfill p_i(t+1)& =& p_i(t)+\frac{a}{2\pi }\mathrm{sin}[2\pi \theta _i(t)]\hfill & \hfill (\mathrm{mod}\mathrm{\hspace{0.33em}1}),\end{array}`$ (7)
where $`t=1,2,\mathrm{}`$, $`i=1,2,\mathrm{},N`$, $`b`$ is the coupling constant (in the following we will fix $`b=2`$), and $`a`$ is a parameter that we will use to control chaoticity; $`\theta _i`$ may be regarded as an angular variable and $`p_i`$ as an angular momentum (defined on a compact set though). Notice that in order to describe a system with a phase space of finite size we are considering, as usual, the map on the torus ($`\mathrm{mod}\mathrm{\hspace{0.33em}1}`$). This map is typically chaotic for large values of $`|a|`$, while for small $`|a|0`$ a complex phase-space structure develops of regular regions surrounded by a stochastic web. The stochastic web is connected for $`d>2`$ (i.e., $`N>1`$) and disconnected for the marginal case $`d=2`$ (i.e., $`N=1`$) (see, e.g., ). Intentionally we introduced a global coupling among the $`N`$ standard maps: to obtain time evolution properties similar to those of the long-ranged interaction many-body Hamiltonian that we shall discuss below.

Our coupled maps display a two-plateau structure with striking analogies with both the one-dimensional map discussed in Section 2 and with the HMF model that we analyze in section 4. As in the previous section, we start by considering an ensemble of $`M`$ independent copies of map (7). In trying to deduce statistical properties directly from the dynamical behavior of a Hamiltonian system – with diagonal kinetic matrix and zero average momentum – it is common practice to identify the temperature with the average squared momentum per particle (see, e.g., ). Because we study here situations with nonzero ”bulk” motion, a natural measure of temperature is given by the variance of the total angular momentum,

$$T(t)\frac{1}{N}\underset{i=1}{\overset{N}{}}\left(p_i^2(t)p_i(t)^2\right),$$
(8)
where $``$ means ensemble average. As the classical microcanonical ensemble is based on the equal a priori probability postulate, we will call BG temperature the temperature associated to the uniform distribution in phase space:

$$T_{BG}\frac{1}{N}\underset{i=1}{\overset{N}{}}\left[_0^1𝑑p_ip_i^2\left(_0^1𝑑p_ip_i\right)^2\right]=1/120.083(N).$$
(9)
Of course, if the system is sufficiently chaotic or if the initial ensemble is sufficiently close to equilibrium (i.e., to the uniform distribution) the temperature (8) rapidly relaxes to $`T_{BG}`$ as a function of the iteration time. Nevertheless, if the phase space presents complex structures and the initial ensemble is very far-from-equilibrium, partial barriers (see, e.g.,) may confine the ensemble in some limited region for quite long times, before allowing for a larger occupation of phase space. Consequently, it appears a plateau with $`T=T_{QSS}T_{BG}`$.

In close analogy with the case of section 4, for our out-of-equilibrium initial conditions we set all momenta $`p_i`$ randomly distributed inside a properly chosen small interval of phase space (the so-called “water bag” initial conditions), while the angles $`\theta _i`$ are distributed at random inside the whole interval $`[0,1]`$. Fig. 3(a) displays the case of a single standard map ($`N=1`$). For this marginal case the formation of total barriers for $`aa_c=0.971635406\mathrm{}`$ prevent the formation of a second plateau below this critical point. We also notice that unless $`a`$ is large enough, the temperature associated to the second plateau is smaller than $`T_{BG}`$. However, the $`2`$-dimensional case is very useful for visualizing the mechanism for plateau formation (see Fig. 3(b)). Fig. 4(a) indicates that in $`d=4`$ the topology of phase space differs from $`d=2`$. Now total barriers for diffusive processes do not exist, and it is possible to have out-of-equilibrium initial data that leads to a dynamics with a $`T_{QQS}`$ temperature two-plateau of long duration before a crossover to $`T_{BG}`$. Interestingly enough, during the first plateau the geometry of the ensemble is characterized by a nontrivial fractal dimension ($`d_f2.7<4=d`$) that at equilibrium crosses over to the dimension of the whole phase space (see Fig. 4(b)). As shown by the new results presented in Fig. 5(a) for the $`200`$-dimensional phase space obtained by coupling $`N=100`$ standard maps, the familiar two-plateau structure persists for large values of $`N`$. Preliminary results suggest that even for large values of $`N`$ a first plateau develops with $`d_f/d<1`$. This fact strongly reinforces the picture that there are nonzero-measure sets of initial conditions such that the system is maintained dynamically in a subset of the whole phase space, and only later occupies the entire allowed space, and adopts a BG-like statistics. This behavior is the more pronounced the closer to zero is the largest Lyapunov exponent.

It is important to stress that in all these cases the crossover time to BG statistics diverges as $`a0`$ (or as $`aa_c`$ for $`d=2`$), as it is shown in Fig. 5(b) and in the inset of Fig. 4(a) and 3(a). The consequence of this is that if the $`t\mathrm{}`$ limit is taken after the $`a0`$ (or $`aa_c`$ for $`d=2`$) limit, the temperature associated to these ensembles is forever $`T_{QSS}T_{BG}`$. The same feature is revealed if, fixing $`a`$, we analyze the scaling of the crossover time for a large $`N`$ (see Fig. 6(a) where we have taken $`a=0.4`$). Moreover, the exponent in the scaling relation $`t_cN^{1.2}`$ is very close to what have been found for the QSS of the HMF model: $`t_cN`$ . It is also interesting to notice that, for the case we are discussing, this scaling is effective only for large value of $`N`$ (say, $`N>20`$). For small $`N`$ a more intricate, richer behavior appear (see Fig. 6(b)). In particular, for $`N=5`$ we observe for the first time in the case of globally coupled symplectic maps evidence of many-plateau structures (like those in section 2).

## 4 QSS in the Hamiltonian Mean Field model

Along the previous lines of thought we can draw our attention now to a more complex dynamical scenario, that is, Hamiltonian many-body dynamics. In analogy with the globally-coupled symplectic maps of the previous section, we will here focus on long-range Hamiltonians. Long-ranged interacting systems constitute nowadays an important subject in many areas of physics, e.g. astrophysics, hydrodynamics, nuclear physics, Bose-Einstein condensates, atomic clusters, plasma physics, and several others . If interaction terms decay as $`r^\alpha `$ and $`d`$ is the dimension of the system, we say that a classical system has long-range interactions when $`\alpha /d>1`$. Long range interactions induce long range correlations, and, as a consequence of this, it has been found that these systems show, in a variety of situations, sensible deviations from the BG equilibrium behavior.

In this section we will overview some known aspects and present as well new results for the HMF model . This model has been largely considered in the literature, as it has been found to display a surprisingly rich behavior (see, e.g.,  for a review). Here we are particularly interested in the appearance, under specific initial conditions (see below), of QSSs that after a certain amount of time cross over to the BG equilibrium. During the QSSs a variety of anomalous behaviors has been detected. Among them, vanishing Lyapunov spectrum, Lévy walks and anomalous diffusion, anomalous distribution of the momenta, non-commutability of the infinite time limit with the thermodynamical limit (see, e.g.,  and references therein), and, more recently, aging  and glassy behavior . As some connections have already been established with nonextensive statistical mechanics , these QSSs are candidates for a situation that is correctly described by this formalism.

The HMF Hamiltonian has the following form:

$$H=K+V=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2}+\frac{1}{2N}\underset{i,j=1}{\overset{N}{}}\left[1cos(\theta _i\theta _j)\right],$$
(10)
where $`\theta _i[0,2\pi )`$ is the $`i`$th angle and $`p_i`$ is the conjugate variable representing the angular momentum (with unit inertial momenta). It is an inertial version of the XY ferromagnetic spin model where the interaction terms couple globally all the spins. Because of the presence of the kinetic part in the Hamiltonian, this model (in opposition with the usual Ising-like spin models) is naturally provided of a dynamics. Also note that it is common use (but not strictly necessary ) to divide the potential term by $`N`$ in order to make the Hamiltonian (formally) extensive .

This system can be analytically solved within the BG canonical formalism, which predicts a second-order phase transition from a low-energy ferromagnetic phase with magnetization $`m|𝐦|0`$, where $`𝐦_{i=1}^N𝐦_i/N`$ (with $`𝐦_i=(\mathrm{cos}\theta _i,\mathrm{sin}\theta _i)`$), to a high-energy phase with all spins uniformly distributed on the interval $`[0,2\pi )`$ (and consequently $`m=0`$). The critical point is at specific energy $`u=u_c0.75`$, and the complete equilibrium caloric curve $`T_{BG}(u)`$ can be exactly obtained . On the other hand, it is possible to numerically integrate the Hamilton equations at fixed total energy, therefore allowing an analysis in the microcanonical setup. Results show nonequivalence between the two approaches for a range of energy densities below the critical point and for a certain class of initial conditions (which includes the so-called “water bag” initial conditions, characterized by a random distribution of the angular momenta around zero). Numerical simulations show in fact a first regime where the system stays in a non-equilibrium metastable state. It is possible to distinguish this QSS by observing the evolution of the dynamical temperature, here defined in the standard way $`T(t)2K(t)/N`$ since the total angular momentum is zero. After a short transient the system evolves into a first stage with a temperature $`T_{QSS}<T_{BG}`$ and it is only at longer times that it relaxes to the predicted equilibrium value. Very similarly to the case discussed in the previous section, the duration of the QSS diverges linearly as $`N`$ goes to infinity .

Inside this framework, we present numerical results using a textbook construction (due to Gibbs in fact) of the canonical ensemble as a subsystem of the microcanonical one. It is a purely dynamical approach, in contrast to the analytical derivation mentioned above. We consider in fact a subset of $`M`$ spins of the total $`N`$ spins in the isolated (microcanonical) system, and study the temperature evolution of this subset, particularly during the metastable state. Temperature is defined as $`T_M(t)2K_M(t)/M`$, where $`K_M`$ is the kinetic energy of the $`M`$ spins. Different criteria may be used in order to choose the $`M`$ subsystem. We show two examples, namely selecting the highest and the lowest temperature subsystem. Result are illustrated in Fig. 7 for $`u=0.69`$ : Two curves coming from above and below, respectively. In both cases $`T_M`$ is seen to relax after some time to $`T_N`$, *within the metastable state.* Afterwards, both temperatures evolve sharing the same value and perform together the crossover to the predicted $`T_{BG}`$. This result is repeated for two system sizes, $`M=100`$ (solid line) and $`M=500`$ (dashed line), yielding the same behavior. It is important to stress that in the case of the hot subsystems, during the relaxation to $`T_{QSS}`$ the temperature value $`T_{BG}`$ is crossed with no signs of relaxation to it. Simulations are made using the $`4`$th order symplectic Neri-Yoshida integrator  with energy conservation $`\mathrm{\Delta }U/U10^4`$, and “water bag” initial condition. This result gives a first step into the study of thermal meta-equilibrium of two systems, i.e., the zeroth principle of thermodynamics for the QSS and, furthermore, allows for the possibility of a generalized canonical treatment of the metastable state, which constitutes in fact an ongoing study .

We also present a new analysis of how the PDF of the momenta of the $`M`$-subset evolves with time. Results are given in Fig.8. We observe that the PDF for the $`M`$-subset essentially coincides with that of the total $`N`$-system. After a transient of order $`t10^3`$, the PDF stabilizes into an non-Maxwellian distribution that lasts up to the crossover time $`t10^6`$. During this whole metastable stage the PDF is well fitted by a $`q`$-Maxwellian with a cut-off, with temperature $`T=T_{QSS}`$ and $`q3.7`$. The value of $`q`$ that we find is different from that calculated in . This may be due to the fact that we are using slightly different (more ordered) initial conditions or to finite-size effects. Our present conjecture is that not only the value of $`q`$ can change when the thermodynamic limit is adequately approached (first $`N\mathrm{}`$, then $`M\mathrm{}`$, and finally $`t\mathrm{}`$), but also the location of the cut-off (in momenta) may go further away, towards infinity. After the crossover time the PDF leaves the QSS anomalous distribution and relaxes to the equilibrium quasi-Maxwellian with $`T=T_{BG}`$ (strictly Maxwellian, with no cutoff at all, only for $`N\mathrm{}`$).

## 5 Summary and discussion

Let us summarize our results. We have reviewed and analyzed three quite different dynamical systems: First, the logistic map at the edge of chaos with additive noise, secondly a system composed of globally coupled standard maps, and finally the HMF model. In all three cases a common dynamical feature is found: The occurrence, for certain classes of initial conditions, of metastable states that after a certain time evolve into the stable equilibrium state. These metastable or quasistationary states appear to be intimately associated to a partial, non-trivial, occupation of phase space, typically (multi)fractal (scale-free and perhaps network-like). This feature retains the system in a non-ergodic state (while held in reserve from entering into all regions of the phase space available to the system). A detailed description for this situation can be given for the simpler case of the logistic map with additive noise at the onset of chaos, i.e., the interplay between a multifractal attractor and a multifractal repeller. This circumstance can also be recognized in the more complex system of two coupled standard maps, since our calculations clearly show the occurrence of quasistationary states and their crossover to final standard chaotic behavior. As we have seen too, the characteristic two-plateau feature is maintained for the case of hundreds of coupled standard maps, and so, we may propose the hypothesis that a fractal occupation holds even for a many-dimensional phase space. The similarities observed between globally coupled symplectic maps and the HMF model indicate that the same state of affairs possibly holds in the case of the QSS detected in the HMF model, a system for which a direct analysis of the $`\mathrm{\Gamma }`$-space is much harder.

The importance of these results is underscored by the fact that the QSSs do become permanent when conventional limits are taken in a specific order, e.g. the thermodynamic limit before the infinite time limit in the HMF model. There is a nonuniform convergence feature in the dynamical evolution that leads to an atypical stationary state. Non-commutability of limits occurs too in the case of the maps. For the coupled standard maps this is appreciated when the thermodynamic limit is replaced (for fixed $`a`$) by the limit of infinite number of coupled maps, or when the chaoticity parameter $`a`$ tends (for fixed $`N`$) towards a specific value ($`a=0`$ for $`d>2`$, and $`a_c`$ if $`d=2`$); while in the case of the logistic map it is the limit of vanishing noise that plays an analogous role. Several other common features are found for the QSSs occurring in these different systems. For instance, in all of them the whole Lyapunov spectrum vanishes, which seems to pinpoint the cause of the observed anomalies. Also, aging and other glassy properties are both present in the logistic map with additive noise and in the HMF model. These common features suggest a deeper, basic, connection among these systems. As we have seen, the usual BG statistical theory appears to be inadequate in explaining the metastable states, and that a generalization of the usual BG statistics is required. In this respect, the nonextensive theory  comes out as the strongest candidate to meet this test.

## Acknowledgments

We acknowledge C. Anteneodo and A. Rapisarda for useful discussions and comments. FB, LGM and CT have benefitted from partial support by FAPERJ, CNPq and PRONEX (Brazilian agencies). APM was partially supported by SECYT-UNC (Argentinean agency). AR was partially supported by CONACyT and DGAPA-UNAM (Mexican agencies).