text:
# Fidelity Decay as an Efficient Indicator of Quantum Chaos

## Abstract

Recent work has connected the type of fidelity decay in perturbed quantum models to the presence of chaos in the associated classical models. We demonstrate that a system’s rate of fidelity decay under repeated perturbations may be measured efficiently on a quantum information processor, and analyze the conditions under which this indicator is a reliable probe of quantum chaos and related statistical properties of the unperturbed system. The type and rate of the decay are not dependent on the eigenvalue statistics of the unperturbed system, but depend on the system’s eigenvector statistics in the eigenbasis of the perturbation operator. For random eigenvector statistics the decay is exponential with a rate fixed precisely by the variance of the perturbation’s energy spectrum. Hence, even classically regular models can exhibit an exponential fidelity decay under generic quantum perturbations. These results clarify which perturbations can distinguish classically regular and chaotic quantum systems.

Over the last two decades a great deal of insight has been achieved regarding the manifestations of chaos and complexity in quantum systems. We are interested in the problem of identifying such signatures in the context of quantum simulation on a quantum information processor (QIP). It is known that a QIP enables efficient simulation of the dynamics of a wide class of quantum systems Lloyd ; ZalkaWeisner . In the case of quantum chaos models, quantum simulation provides an exponential speedup over direct classical simulation Schack ; Shep ; Casati . Recently, the quantum baker’s map has been implemented using a nuclear magnetic resonance QIP Yaakov . These developments highlight the importance of devising efficient QIP methods for the measurement of quantum chaos signatures and related properties of complex quantum systems.

Perhaps the most established signature of quantum chaos is given by the (nearly) universal correspondence between the eigenvalue BT77BGS84  and eigenvector IZR87KMH88HZ90  statistics of quantized classically chaotic systems and those of the canonical ensembles of random matrix theory (RMT). Unfortunately direct detection of these spectral signatures is algorithmically inefficient by any known technique. However, following the original observation of Peres Peres , some recent work has demonstrated that, under sufficiently strong perturbation, the fidelity exhibits a characteristic exponential inthe case of classical chaotic systems JP01Cucchietti01 ; Jacquod01 ; Casati01 ; Prosen02 . Below, we demonstrate that this charecteristic decay may be measured by an efficient algorithim, and analyze in detail how the observation of this decay may be applied as an indicator of canonical RMT statistics (quantum chaos) in the unperturbed system.

By first considering RMT models we show that Wigner-Dyson fluctuations in the system eigenvalue spectrum are not necessary to produce this characteristic decay. More importantly, we find that the canonical RMT statistics for the system eigenvector components, expressed in the eigenbasis of the perturbation operator, are sufficient to produce the characteristic exponential decay. These observations are checked in the case of a dynamical model, where we demonstrate that an exponential fidelity decay can arise not only for a classically chaotic system, but also for classically regular system under a generic choice of perturbation. Our analysis implies specific restrictions that must be imposed on the choice of applied perturbation operator in order to extract useful information about the statistical signatures of quantum chaos in the unperturbed system.

The aim of this Letter is to characterize certain static properties of a unitary map, $`U`$, by observing the rate of divergence between a fiducial Hilbert space vector evolved under this map, $`|\psi _u(n)=U^n|\psi (0)`$, and the same initial vector evolved under this map but subject also to a sequence of small perturbations $`|\psi _p(n)=(U_pU)^n|\psi (0)`$, where $`U_p=\mathrm{exp}(iV\delta )`$ is some unspecified perturbation operator and $`n`$ denotes the number of iterations. The fidelity,

$$O(n)=|\psi _u(n)|\psi _p(n)|^2,$$
(1)
provides a natural indicator of this divergence. The value of $`O(n)`$ may be determined by an efficient algorithm on a QIP as follows. We start by preparing the fiducial state $`|\psi (0)=U_o|0^{n_q}`$, where $`n_q`$ is the number of qubits required to span the system’s Hilbert space $`N=2^{n_q}`$. As shown below, the choice of fiducial state is not critical and the computational basis states, which are simplest to implement, are a convenient set. After applying the sequence $`(U^{})^n(U_pU)^n`$, the system register contains a final state $`|\psi (n)`$. The circuit implementation of $`U`$ requires only Poly($`n_q`$) operations for the simulation of a wide class of quantum systems to arbitrary accuracy Lloyd ; ZalkaWeisner ; Shep ; Schack ; Casati . Here we observe that $`|\psi _p(n)|\psi _u(n)|=|\psi (n)|\psi (0)|=|\psi _f|0|`$, where the state $`|\psi _f`$ is obtained by time-reversing the initial state preparation $`|\psi _f=U_o^{}|\psi (n)`$. The magnitude of $`O(n)`$ is then determined from sampling the population of the state $`|0`$. The entire algorithm therefore scales as Poly($`n_q`$).

Recently, Jacquod and coworkers Jacquod01  observed that for $`|\psi (0)=|v_o`$, an initial eigenstate of $`U`$ with eigenphase $`\varphi _o`$, the fidelity relates to the local density of states (LDOS), $`\eta (\varphi _o\varphi _m^{})=|v_o|v_m^{}|^2`$, via Fourier transform,

$$O(n)=|\underset{m}{}\eta (\varphi _o\varphi _m^{})\mathrm{exp}(i(\varphi _o\varphi _m^{})n)|^2.$$
(2)
In the above $`|v_m^{}`$ are eigenstates of the perturbed map $`U_pU|v_m^{}=\mathrm{exp}(i\varphi _m^{})|v_m^{}`$. These observations locate the decay of $`O(n)`$ in an existing theoretical framework. In the non-perturbative regime $`\sigma /\mathrm{\Delta }1`$, where $`\sigma ^2=\delta ^2\overline{V_{mn}^2}`$ denotes a typical off-diagonal matrix element and $`\mathrm{\Delta }`$ is the average level spacing, previous work suggests that when the perturbed system is complex the LDOS is typically Lorentzian Wigner ; FCIC96 ,

$$\eta (\varphi _o\varphi _m^{})\frac{\mathrm{\Gamma }}{(\varphi _o\varphi _m^{})^2+(\mathrm{\Gamma }/2)^2}$$
(3)
with width $`\mathrm{\Gamma }=2\pi \sigma ^2/\mathrm{\Delta }`$ determined by the Fermi golden rule (FGR). From (2) and (3) one expects the exponential decay,

$$O(n)\mathrm{exp}(\mathrm{\Gamma }n).$$
(4)
The onset of the exponential decay (4) has been confirmed recently in a few classically chaotic models Jacquod01 ; Casati01 , though the rate is not always given by the golden rule Wisniacki02 ; JP01Cucchietti01 . However, in the case of integrable $`U`$ the situation is less clear since under some perturbations the LDOS is known to take on a Lorentzian shape FCIC96 ; BGI98 .

Below, we examine which statistical properties of the unperturbed system lead to the FGR decay, and how this characteristic decay depends on the properties of the perturbation operator. This approach is motivated by the context of QIP simulation, in which the eigenbasis of this perturbation may be mapped onto an arbitrary basis of the simulated system $`U`$ Cory99 . Our choice of perturbation eigenvalue structure is motivated from the perspective of quantum control studies. Specifically, we consider

$$U_p=\mathrm{\Pi }_{j=1}^{n_q}\mathrm{exp}(i\delta \sigma _z^j/2),$$
(5)
where $`\sigma _z`$ is the usual Pauli matrix and (5) therefore corresponds to a collective rotation of all the qubits by an angle $`\delta `$. Eq. (5) is a model of coherent far-field errors Fortunato , and for this type of error model a better understanding of the fidelity decay is a subject of intrinsic interest.

As a first test of the LDOS/FGR framework in the case of the qubit perturbation (5) we evaluate the fidelity decay for a map $`U=U_{\text{CUE}}`$ drawn from the circular unitary ensemble (CUE). These matrices form an established model for classically fully chaotic time-periodic systems (without additional symmetries) since these systems (almost always) exhibit the same characteristic (Wigner-Dyson) spectral fluctuations and eigenvector statistics as the CUE Haake01 . The randomness of the system eigenvectors enables a system-independent estimate of the rate $`\mathrm{\Gamma }`$ of the FGR decay. Since the components of random eigenstates are distributed uniformly over the basis states and uncorrelated with the distribution of eigenvalues, the second moment of the matrix elements $`V_{mn}`$ may be directly evaluated,

$$\overline{V_{mn}^2}=\overline{\lambda ^2}/N$$
(6)
assuming $`\overline{\lambda }=0`$ and where $`\overline{\lambda ^2}=N^1_{i=1}^N\lambda _i^2`$ denotes the variance of the eigenvalues of $`V`$. As a result, the rate of the FGR decay is determined by the eigenvalues of the perturbation,

$$\mathrm{\Gamma }=\delta ^2\overline{\lambda ^2},$$
(7)
where we have used $`\mathrm{\Delta }=2\pi /N`$. For the qubit perturbation (5) the variance of the eigenvalues has a simple form,

$$\overline{\lambda ^2}=\frac{1}{N}\underset{k=0}{\overset{n_q}{}}\left(\frac{2kn_q}{2}\right)^2C_k^{n_q},$$
(8)
where the $`C_k^{n_q}`$ are binomial coefficients. Using our RMT estimate (7), for $`n_q=10`$, the rate is,

$$\mathrm{\Gamma }=2.50\delta ^2.$$
(9)
While a CUE map may be generated on a QIP using the gate decomposition devised in Ref. Zyc , for our numerical study we construct $`U=U_{\text{CUE}}`$ directly from the eigenvectors of a random Hermitian matrix. Since computational basis states are easiest to prepare in the QIP setting, we consider the fidelity decay for both single computational basis states and averages over 50 such states. The behaviour of the fidelity decay for a matrix typical of CUE is displayed in Fig. 1. The three perturbation values displayed in the figure are chosen near the onset of the non-perturbative regime ($`\delta >0.1`$) and it is evident that the fidelity decay even for individual computational basis states exhibits FGR decay (4) at the expected rate. The FGR decay persists for a time-scale $`\mathrm{\Gamma }^1\mathrm{log}(N)`$ until saturation at a time-average that decreases as $`1/N`$.

We next consider $`O(n)`$ for the Gaussian unitary ensemble (GUE) in order to clarify the relationship between the FGR decay and the distinct statistical features of RMT that represent signatures of quantum chaos. The GUE consists of Hermitian matrices with independent elements drawn randomly with respect to the unique unitarily-invariant measure Haake01 . GUE forms the relevant RMT model for the important class of chaotic or complex autonomous Hamiltonian systems that are unrestricted by any additional symmetries. We may examine the sensitivity to perturbations for the GUE by constructing the unitary operator $`U_{\text{GUE}}=\mathrm{exp}(iH_{\text{GUE}}\tau )`$, where $`\tau `$ is a time-delay between perturbations. We consider the same perturbation as for the CUE case. For sufficiently small $`\tau `$ the propagator approaches identity and the overlap decay is dominated by the perturbation operator, $`O(n)=|\mathrm{exp}(i\delta V)|^2+O(\delta ^2\tau ^2n^2)`$. This behaviour is demonstrated in Fig. 2 for $`\tau =0.001`$ and $`\tau =0.01`$ and with $`\delta =0.3`$. For larger values of $`\tau `$ the fidelity decay under $`U_{\text{GUE}}`$ obeys the FGR with the RMT rate (9).

The important point is that for sufficiently large $`\tau `$ the eigenphases of $`U_{\text{GUE}}`$ become spread pseudo-randomly in the interval $`[0,2\pi )`$. Under these conditions the eigenphases of the map $`U_{\text{GUE}}`$ exhibit the Poissonian spectral fluctuations that are characteristic of classically integrable (time-periodic) systems. We checked the nearest-neighbor spacing distribution of $`U_{\text{GUE}}`$ and found that for $`\tau =100`$ the statistics are in excellent agreement with the Poissonian distribution $`P(s)\mathrm{exp}(s)`$ (see inset to Fig. 2). However, the eigenvectors of $`U_{\text{GUE}}`$ are random (by construction) and independent of $`\tau `$ (for finite $`\tau `$). From these observations it is clear that the presence of Wigner-Dyson spectral fluctuations in the implemented $`U`$, which comprises the only basis-independent criterion of quantum chaos, is not actually necessary for the onset of exponential (FGR) decay at the rate (7). This suggests that it is the RMT statistics of the eigenvectors of $`U`$ that lead to the FGR decay with the RMT rate.

We next consider the fidelity decay for the quantum kicked top, which is an exemplary dynamical model of quantum chaos HKS87 ; Haake01 . The kicked top is a unitary map $`U_{\text{QKT}}=\mathrm{exp}(i\pi J_y/2)\mathrm{exp}(ikJ_z^2/j)`$ acting on the Hilbert space of dimension $`N=2j+1`$ associated with an irreducible representation of the angular momentum operator $`\stackrel{}{J}`$. In previous fidelity decay and LDOS studies the choice of perturbation has usually been tied to a physical coordinate of the system $`U`$. We first follow this convention and identify the eigenbasis of the perturbation (5) with the eigenbasis $`|m_j`$ of the system coordinate $`J_z`$ (where $`m_j=\{j,\mathrm{},j\}`$). In Fig. 3 we compare the fidelity decay for the chaotic and regular regimes of the kicked top for averages over 50 initial computational basis states. The fidelity decay for the chaotic top ($`k=12`$) is well described by the FGR prediction (4) and the RMT rate (9), whereas the regular top ($`k=1`$) shows a slower non-exponential decay rate. Similarly, if we associate the perturbation eigenbasis with the basis of the $`J_y`$ coordinate, the fidelity decay for the chaotic top remains in agreement with the RMT rate and the regular top again exhibits non-exponential decay, though in this case with a faster decay than the RMT rate (9). However, we now demonstrate that an exponential decay at the RMT rate arises even for the regular kicked top when the qubit perturbation (5) is diagonal in a generic basis relative to the eigenbasis of $`U_{\text{QKT}}`$. Specifically, we leave the perturbation eigenvalue spectrum unchanged but set

$$U_p=T\left[\mathrm{\Pi }_j\mathrm{exp}(i\delta \sigma _z^j/2)\right]T^1,$$
(10)
where $`T`$ is drawn from CUE. As demonstrated in the inset to Fig. 3, under this type of perturbation the fidelity decay for the regular top is indistinguishable from that of the chaotic top and is very accurately described by the FGR at the RMT rate (9).

The sensitive dependence of the type of fidelity decay on the eigenbasis of the applied perturbation suggests a close connection with the basis-dependence of the eigenvector statistics of classically regular quantum models. Expressed in a generic quantum basis, the eigenvectors of any quantized classical system $`U`$ will have randomly (Gaussian) distributed components. In contrast, in the eigenbases of the system coordinates the components of classically chaotic and integrable systems are known to be different IZR87KMH88HZ90 , with the former Gaussian distributed and the latter exhibiting substantial deviation from the canonical Gaussian distribution. In light of this connection, in the case of quantized classical models we infer that exponential (non-exponential) fidelity decay can be correlated with the presence (absence) of characteristic RMT spectral fluctuations in the unperturbed system provided that the applied perturbation commutes with a system coordinate.

In summary, we have shown that the fidelity decay may be measured efficiently on a QIP. We then examined which statistical properties of the unperturbed system determine the type and rate of the decay. In the case of random unitary and Hermitian matrices, as well as a classically chaotic dynamical model, we have shown that the fidelity decays exponentially with a characteristic rate given precisely by the variance of the perturbation’s eigenspectrum. The occurrence of the exponential decay is not directly dependent on the Wigner-Dyson fluctuations of the unperturbed spectrum, but does depend sensitively on the RMT statistics of the system eigenvectors in the eigenbasis of the applied perturbation. Hence, the fidelity decay for both classically regular and chaotic dynamical systems is given by the FGR under all but a small subset of unitary perturbation operators. In the case of classical models, we conclude that the fidelity decay provides a reliable indicator of RMT statistics (quantum chaos) in the unperturbed system only when the applied perturbation is restricted to the subset of perturbations that commute with a classical coordinate.

We are grateful to K. Zyczkowski for suggesting the relevance of the eigenvector statistics, Ph. Jacquod, E.M. Fortunato, C. Ramanathan, and T.F. Havel for helpful discussions, and DARPA and the NSF for financial support.


format: text

text:
# Ordered phase and phase transitions in the three-dimensional generalized six-state clock model

## I INTRODUCTION

The systems of Z6 symmetry in the three dimensions (3D) have attracted interests because of their peculiar nature of orderings. In particular, there have been many works on the nature of intermediate temperature phases.

In two dimensions, the six state clock model shows three phases. At intermediate temperature, the discreteness of the six states is irrelevant and the system shows the same properties as that of XY model, which is the Kosterlitz-Thouless phase where the long range order does not exist but the correlation length of the order diverges .

In three dimensions, the six state clock model also shown natures which are similar to the ordered phase of the XY model. However it found that there exists only a phase where one of the six state of the model is chosen as the long range order . The apparent behavior is due to the large fluctuation at the higher temperature region of the ordered phase. However, if we generalize the energy structure of the model, a new type of intermediate phases appears.

Ueno et al have introduced three-dimensional generalized six-state clock model (3D-6GCL model). The 3D-6GCL model can be regarded as prototypes of wide range of Z6 models. That is, it represents various categories of the Z6 models according to the values of the energy parameters of the model. The existence of intermediate phase have attracted interest for the following models: three-dimensional three-state antiferromagnetic Potts model (3D-3AFP model), and stacked triangular antiferromagnetic Ising model (STAFI model). Although, the existence of various type of intermediate phases for these models have been proposed, there is no general understanding on the nature of the intermediate phase(s) . If we consider only one universality class for models with Z6 symmetry, we conclude that the intermediate phase is an apparent phase as we found in the regular six state model. However, it has been clarified that the situation is not so simple and interesting new type of phase exists in the generalized case. In fact, the 3D-6GCL model was introduced in order to understand general properties of other Z6 models . In the present paper, we study the 3D-6GCL model with a set of the parameters which is considered to correspond to the STAFI model and the 3D-3AFP model.

The Hamiltonian of the 3D-6GCL model is given by

$`={\displaystyle \underset{i,j}{}}\epsilon _0\delta _{n_in_j}+\epsilon _1\delta _{n_in_j\pm 1}+\epsilon _2\delta _{n_in_j\pm 2}+\epsilon _3\delta _{n_in_j\pm 3},`$ (1)
where $`i,j`$ runs over nearest neighbor pairs, and $`n_i`$ is the spin variable which takes one of $`1,2,\mathrm{}`$ and $`6`$. We set $`\epsilon _0=0\epsilon _1\epsilon _2\epsilon _3`$ as in Fig.1. Namely, we only study ferromagnet cases where all the spins occupy the same state in the ground state. We expect that this ferromagnetic ordered phase appears at the low temperature. This ordered phase is called completely ordered phase (COP). Since the COP has six degenerate states, the order parameter space of this model is illustrated by a hexagon as in Fig.2. The bold solid line represents the one of the six COP states. The dotted line represents a mixture of neighboring two COP states which is a candidate of an intermediate phase, which we call incompletely ordered phase (IOP1). As to the distribution of the spin state in the incompletely ordered phase, there have been proposed the two kinds depicted in Fig.3. The IOP1 discussed above corresponds to Fig.3(a). The type of Fig.3(b) was called IOP2 which was one of the possible mixed state. For $`\epsilon _1=0`$, the COP does not appear because the ground state is macroscopically degenerate. This degenerate state is IOP1 and corresponds to the ground state of the 3AFP model. Taking the approach of physical percolation to phase transitions, Ueno rigorously proved that for $`\epsilon _1=0`$, at least one kind of IOPs can exist at $`T<\epsilon _2=\epsilon _3`$ . For $`\epsilon _1>0`$ , Ueno and Kasono proposed that this model has two incompletely ordered phases (IOP1, IOP2) besides the COP by using the Monte Carlo twist method (MCTM) . In this paper, we reexamine the nature of the intermediate phase(s) in detail.

On the 3D-3AFP model, so far various types of phases at intermediate temperatures have been proposed in the previous studies, such as the permutationally symmetric sublattice phase (PSS phase) which corresponds the IOP2  or the gapless phase which is similar to the ordered phase of the three dimensional XY model which we call XY phase . Recent theoretical and numerical studies on the 3D-3AFP model have revealed that the existence of these phases are fake owing to the finite size effect . There is only a low-temperature phase that corresponds to the IOP1.

We will show that the correlation length of the fluctuation of the order parameter in IOP1 is finite. However the correlation length just below the critical point is so long that the phase looks like the ordered phase of the 3D XY model where the correlation length is infinite although the long range order exists. In order to distinguish the difference between the XY phase and IOP1 we have to be careful in the examination of properties of the phase.

In this paper, we study a 3D-6GCL model by using the Monte Carlo (MC) simulation and clarify properties of the phase at intermediate temperature and natures of the phase transitions. In this study, we choose the energy parameters of the 3D-6GCL model as $`\epsilon _0=0`$, $`\epsilon _1=0.1`$ and $`\epsilon _2=\epsilon _3=1.0`$ which are the same as those in .

In order to examine the possibility of many ordered phases by the MCTM, one needs a full set of appropriate boundary conditions which cause different types of domain-walls, and calculate size dependence of the excess free energy. In the present study, we deny the existence of IOP2 as intermediate phase from the dependence of the excess free energy on the boundary conditions. The intermediate phase shows apparent property of the XY phase where Z6 discreteness is irrelevant. However, we also deny the existence of XY phase in the thermodynamic limit. This apparent property is attributed to the long correlation length at the high-temperature region of the intermediate phase. Finally, we conclude that only IOP1 exists as the intermediate temperature. Appearance of the fake XY phase near the critical point has been also found in the six-state clock model, where the Z6 discreteness is very weak in the ordered phase near the critical point. However the Z6 discreteness is always relevant in the ordered phase .

Moreover, we examine the properties of the phase transitions. We consider that the high temperature phase transition belongs to 3D-XY universality class. In order to confirm whether this transition belongs to 3D-XY universality class, we study this transition by using non-equilibrium relaxation method (NERM) and estimate the critical temperature and critical exponents. Ueno and Kasono pointed out that the low-temperature transition between the IOP1 and the COP is of the first order because the two phases have no symmetry relation of group and subgroup. We investigate this phase transition by using MCTM and the usual finite size scaling analysis, and confirmed that it is of the first order.

The outline of this paper is as follows. In Sec.II, we study properties of the intermediate phase. In Sec.III, we study that the high temperature transition and in Sec.IV, we investigate the low temperature phase transition, and the Sec.V is for summary and discussion.

## II INTERMEDIATE PHASE

The properties of the phases of the 3D-6GCL model have been studied by the MCTM , which is a simulation method under special boundary conditions to detect the domain-wall excess free energy. We prepare systems with fixed boundaries in one direction and impose a periodic boundary condition in the other directions. As a reference system, we set the same state $`\alpha `$ at the boundaries ($`\alpha `$-$`\alpha `$) in which we expect that the ordered phase $`\alpha `$ appears. We also prepare a system where we set the states $`\alpha `$ and $`\beta `$ at the boundaries ($`\alpha `$-$`\beta `$) in which we expect an interface. The interfacial free energy is defined as the excess free energy between the two systems:

$`\mathrm{\Delta }F^{\alpha \beta }(T,L)=F^{\alpha \beta }(T,L)F^{\alpha \alpha }(T,L),`$ (2)
where $`L`$ is a liner size of the system. Here we consider an $`L\times L\times L`$ system. The size dependence of $`\mathrm{\Delta }F^{\alpha \beta }(T,L)`$ is given by the following asymptotic behavior for large $`L`$

$`\mathrm{\Delta }F(T,L)A(\alpha ,\beta )L^{\psi _{\alpha \beta }(T)},`$ (3)
where $`\psi _{\alpha \beta }(T)`$ is called stiffness exponent. The ordered phase is classified by the value of stiffness exponent as

$`\psi _{\alpha \beta }(T)=\{\begin{array}{c}D1\text{for a domain-wall type interface (}\xi <\mathrm{}\text{)}\hfill \\ D2\text{for a gepless interface (spin-wave type }\xi =\mathrm{}\text{)}\hfill \\ \text{non-integer for a new type on interface}\hfill \end{array}.`$ (7)
When the interface does not appeare, the stiffness exponent becomes negative because $`\mathrm{\Delta }F(T,L)0`$ for $`L\mathrm{}`$.

In the present work, we use three sets of boundary conditions which are shown in Fig.4. There are two kinds of boundary conditions; in $`\varphi _1`$ we fix the boundaries to be two of the low temperature phases. The relative angle (twist angle) between the two phases can be $`\frac{\pi }{3}`$, $`\frac{2\pi }{3}`$ and $`\pi `$. Here we take the cases $`\frac{\pi }{3}`$ and $`\pi `$ which are shown in Fig.4 (a) and 4 (b), respectively. In $`\varphi _2`$ we fix the boundaries to be two of the IOP1s in Fig.4 (c) (twist angle $`\pi /3`$). As shown in Table I, we distinguish the five phases from the signs and the values of the stiffness exponent. We perform MC simulation for the lattices with liner size $`L=8,10`$, and 12 with boundary conditions $`\varphi _1(\pi )`$, $`\varphi _1(\pi /3)`$ and $`\varphi _2(\pi /3)`$. At each temperature, after discarding first $`5000L`$, $`5000L`$ and $`10000L`$ steps, we calculate quantities of interest using data of next $`25000L`$, $`25000L`$ and $`50000L`$ steps for the three sets of boundary conditions, respectively.

In Fig.5, we show the temperature dependence of $`\psi `$ under the boundary conditions. From this figure, we may conclude the existence of the disorder phase, the XY phase, the IOP1 and the COP, as the temperature decreases. The high temperature transition point is estimated to be $`T_\text{C}1.35(5)`$, and the low temperature transition point is to be $`T_\text{F}0.30(5)`$. We see a change around $`T0.7`$ which can be a phase transition between the XY phase and IOP1.

In order to confirm the region of phases, we investigate the size dependence of the stiffness exponents for $`\varphi _1(\pi /3)`$ for larger sizes. The region of the XY phase decreases as $`L`$ becomes large (Fig.6). The size-dependence strongly suggests that the region of the XY phase disappears in the limit $`L\mathrm{}`$ and only the IOP1 exists in the intermediate temperature region.

If we consider the domain wall between the ordered states of the type of IOP1, there is no domain wall in the $`\varphi _1(\pi /3)`$. On the other hand, one domain wall appears in the case of $`\varphi _2(\pi /3)`$, and two domain walls appear in the $`\varphi _1(\pi )`$. Therefore, stiffness exponent for the cases $`\varphi _1(\pi )`$ and $`\varphi _2(\pi /3)`$ should be the same. However, we notice that the difference of the values of stiffness exponents of $`\varphi _1(\pi )`$ and $`\varphi _2(\pi /3)`$ is significant at the IOP1 as we see at Fig.5. The reason for this difference is that the effect of twisting in the system by $`\varphi _1(\pi )`$ is stronger than one of the twisting by $`\varphi _2(\pi /3)`$. Actually, when we examine $`\psi _1(\pi )`$ and $`\psi _2(\pi /3)`$ in much larger lattice sizes and we confirm that the values of the stiffness exponents $`\psi _1(\pi )`$ and $`\psi _2(\pi /3)`$ approach the same value $`\psi =2`$ which indicates the domain-wall type order. In Fig.8, we show a snapshot of the domain-wall of the boundary in $`\varphi _2(\pi /3)`$ of $`L=160`$ at $`T=0.4`$, where we find a localized domain-wall clearly.

Next, following the idea of Oshikawa  for the 3D-3AFP model, we make the finite size scaling plots of $`\mathrm{cos}6\theta `$ by ordinary Monte Carlo method on lattices $`L=12,18,24`$ and $`30`$. Here, $`\mathrm{cos}6\theta `$ is the order parameter for the Z6 symmetry breaking and $`\theta `$ is the angle of magnetization from one of the COP states. After discarding first $`10000L`$ steps, we calculate quantities of interest for next $`40000L`$ steps at each temperature. In Fig.7, we show the scaling plots in the form of the scaling function given by Oshikawa

$`\mathrm{cos}6\theta f(cL^2(T_CT)^{\nu |y_6|}).`$ (8)
Here, we put the value of the critical temperature $`T_\text{C}=1.384`$ which was obtained by using the NERM which will be explain a later. We obtain a good scaling plot when we use $`\nu |y_6|=3.2`$. It is important to notice here that behavior of the order in the region of the IOP1 is expressed by a scaling function using the high temperature critical point $`T_\text{C}`$. This fact indicates that other phases do not exist between the IOP1 and the disordered phase because the fluctuation in the IOP1 phase is controlled by the high temperature critical point. However, it is very difficult to estimate the accurate value of this exponent because the strong finite size effect appears in the intermediate phase even at the largest system which we can calculate. After all, we conclude that the 3D-6GCL model with $`\epsilon _0=0`$, $`\epsilon _1=0.1`$ and $`\epsilon _2=\epsilon _3=1.0`$ has two phases, i.e. IOP1 as the intermediate phase and COP as the low temperature phase.

In order to investigate the nature of IOP1, we examine the following quantity;

$`\chi ={\displaystyle \frac{1}{L^3}}((p_1p_2)^2(p_1p_2)^2),`$ (9)
where $`p_i`$ denotes the number of the spin variables which take the value of the $`i`$-th state. We adopt the boundary conditions $`\varphi _2(0)`$. The asymptotical form of $`\chi `$ is written $`\chi AL^x`$. The exponent $`x`$ depends on the type of the order. When the fluctuation is of the XY phase, we expect $`x=2`$. On the other hand, when the fluctuation has a finite correlation length, $`x=0`$.

We perform MC simulations for the lattices with the liner sizes $`L=8,10,\mathrm{}`$ and $`24`$. For each temperature, after discarding first $`10000L`$ steps, we calculate $`\chi `$ using data of next $`50000L`$ steps. In Fig.9, we show the size and temperature dependence of $`x`$. We find that the exponent $`x`$ decreases in the IOP1 region. Moreover, we perform MC simulation for the lattice up to the liner size $`L=80`$ at $`T=0.4`$ (in the IOP1 region). In Fig.10, we show the temperature dependence of $`\chi `$ at $`T=0.4`$. The value of $`\chi `$ approaches the constant as $`L`$ increases. Thus we consider that the value of $`x`$ decreases and approaches 0 as $`L\mathrm{}`$, and we conclude that the fluctuation in IOP1 is the type of the domain-wall and the correlation length of the fluctuation of the order parameter is finite.

## III HIGH TEMPERATURE PHASE TRANSITION

The high temperature phase transition in Z6 models belongs to the 3D-XY universality class . Indeed, the transitions in the 3D-3AFP model  and the six-state clock model  are found to belong to the 3D-XY universality class. We consider that the high temperature phase transition in the 3D-6GCL model also belongs to 3D-XY universality class.

We investigate this phase transition by the NERM . The NERM is an efficient numerical technique to estimate the critical point and critical exponents from a dynamical process toward the equilibrium state from an ordered state. The decay of the order parameter $`m(t)`$ shows a power-low only at the critical point. Detecting such point, a precise determination of the critical temperature can be done. The asymptotical form of $`m(t)`$ at the critical temperature is written as

$`m(t)t^{\lambda _m}.`$ (10)
To examine asymptotic behavior of $`m(t)`$ clearly, we introduce a local exponent $`\lambda _m(t)`$

$`\lambda _m(t){\displaystyle \frac{d\mathrm{log}m(t)}{d\mathrm{log}t}}.`$ (11)
This exponent $`\lambda _m(t)`$ corresponds to $`\beta /z\nu `$. Further, we consider the following function and the local exponent.

$`f_{mm}(t)\left[{\displaystyle \frac{m^2}{m^2}}1\right],\lambda _{mm}{\displaystyle \frac{d\mathrm{log}f_{mm}(t)}{d\mathrm{log}t}},`$ (12)
where $`\lambda _{mm}(t)`$ corresponds to $`d/z`$. Therefore, we obtain the exponents $`\beta /\nu `$ and $`z`$ independently from these quantities. We simulate the relaxation process starting from an initial state which is set to be the IOP1 state, and measure the magnetization $`m(t)`$. In Fig.11, we show $`\lambda _m(t)`$ for the temperatures near the $`T_\text{C}`$. Here, MC simulation is performed in a lattice $`L=60`$. We show the MC steps and the number of samples in Table II. The curve for $`T=1.385`$ turns up while the curve for $`T=1.383`$ turns down. Therefore we locate the transition temperature in $`1.383<T_\text{C}<1.385`$, denoting $`T_\text{C}=1.384(1)`$.

At this accurate value of $`T_\text{C}`$, we calculate the relaxation of quantities. We perform about 160000 independent runs up to 200MCS and average the process to obtain $`m(t)`$. We obtain the critical exponents $`\beta /\nu `$ and $`z`$ by making use of (11) and (12). In Figs.12 and 13, we show $`z(t)`$ and $`\beta /\nu (t)`$, respectively. From the extrapolated values of $`z(t)`$ and $`\beta /\nu (t)`$, to $`1/t=0`$, we estimate $`z=2.05(5)`$ and $`\beta /\nu =0.515(10)`$.

Moreover, assuming $`T_\text{C}=1.384(1)`$ and $`\beta /\nu =0.515(10)`$, we estimate $`\beta `$ from the scaling plots of data of the magnetization which are obtained from the ordinary equilibrium MC simulation. Seeking the value of $`\beta `$ by which the date collapse into a scaling function, we estimate $`\beta =0.34(1)`$ for the best fit (Fig.14). Since the values, $`\beta =0.34(1)`$ and $`\nu =0.66(4)`$, are close to the 3D-XY universality class ($`\beta =0.345`$, $`\nu =0.669`$), we conclude that this phase transition belongs to the 3D-XY universality class.

## IV LOW TEMPERATURE PHASE TRANSITION

While many researchers have studied the high temperature transition in the Z6 model, but few study has been done for the low temperature one. Ueno and Kasono pointed out that the low temperature transition between the IOP1 and the COP is of first-order because they have no symmetry relation between group and subgroup which is usually seen in the second-order transition.

In this paper we study this problem by a finite-size-scaling analysis of $`\mathrm{cos}6\theta `$. If the transition is of first order, the order parameter is scaled by $`\mathrm{cos}6\theta f(ϵL^d)`$ where $`ϵ=(TT_\text{C})/T_\text{C}`$, $`d`$ is the dimension of the system and $`f(x)`$ is a finite size scaling function. We obtain a good scaling plot at $`T_\text{F}=0.295(5)`$ depicted Fig.15. Thus we confirm this phase transition is of first-order.

We also confirm this result by an analysis of extension of the MCTM . Now, we study the interfacial energy $`\mathrm{\Delta }E`$ is defined in the same way for the excess energy due to the boundary condition causing the interface, eq.(2). In a completely ordered phase, the size dependence of the excess energy becomes $`\mathrm{\Delta }EL^{d1}`$. At the critical point, it becomes $`\mathrm{\Delta }EAL^{1/\nu }`$. We define the stiffness exponent of the energy $`\psi _E`$ as $`\mathrm{\Delta }EAL^{\psi _E}`$. For the 3D Ising model at the critical point $`\psi _E=1.59`$ which is considered to be a large value among the 3D models exhibiting the second order phase transition. We study $`\psi _E`$ for the present phase transition. We use the boundary condition $`\varphi _1(\pi /3)`$. In Fig.16, we show the temperature and the size dependence of $`\psi _\text{E}`$. We confirm that $`\psi _\text{E}`$ approach to 3 which is much larger than the values of $`\psi _E`$ of the second order phase transition in 3D models. For the first order phase transition, it is known that $`\psi _E=D`$ from the argument of wetting. This result also suggests that the present phase transition is of first order.

## V SUMMARY and DISCUSSION

We study successive phase transitions in the 3D-6GCL model and found an intermediate phase and two phase transitions. First, it has been revealed that the intermediate phase is single phase of the IOP1. In the high temperature region of the intermediate phase, the correlation length is very large and the system shows an apparent XY behavior. In this region the stiffness constant is very sensitive to the boundary condition. Therefore, the wrong results were obtained in the previous numerical studies . The present results agree with the phase diagram of the 3D-3AFP model obtained by Kishi and Ueno, and our results of the STAFI model . Thus, we have obtained the same intermediate phase in the three 3D Z6 models, i.e. the 3D-6GCL model, the 3D-3AFP model and STAFI model. Moreover, we found that the mixing rate of the two states in the intermediate phase is steady and does not show anomaly fluctuation of the rate proposed preciously . This IOP1 with the steady mixing of two COP states can be regarded as partially disordered phase proposed by Mekata in the original work on the triangular antiferromagnets. Here we understand that the entropy effect due to the flustration allows different type of intermediate phases in two and three dimensional model.

Secondly, we examine properties of the high temperature phase transition, and the transition temperature is estimated to be $`T_\text{C}=1.384(1)`$ with the critical exponents $`\beta =0.34(1)`$, $`\nu =0.66(4)`$ and $`z=2.05(5)`$. These values of the exponents are close to those of 3D-XY universality class. We conclude that the universality class of this transition belongs to the 3D-XY universality class.

Thirdly, we investigate the low temperature phase transition in the 3D-6GCL model by using the finite-size-scaling analysis and the MCTM. All of the results support that this transition is of first-order. This first-order phase transition is a transition at which the $`\mathrm{cos}6\theta `$ changes from $`1`$ to $`1`$ discontinuity.

## ACKNOWLEDGMENTS

We would like to thank F. Matsubara, M. Oshikawa, N.Ito and Y. Ozeki for useful discussions. The numerical calculations were performed on the supercomputers in the Computer Center of the University of Tokyo.


text:
# Models of Galaxy Clusters with Thermal Conduction

## 1 Introduction

For many years it was thought that radiative losses via X-ray emission in galaxy clusters leads to a substantial inflow of gas, and mass dropout, in the form of a “cooling flow” (see Fabian 1994 for a review). Mass deposition rates were estimated to be as much as several hundred $`M_{}\mathrm{yr}^1`$ in some clusters (e.g., Peres et al. 1998). However, there was little direct evidence for the mass dropout in any band other than X-rays (Fabian, 1994).

Recent high resolution X-ray data from XMM-Newton and Chandra indicate that there is no evidence for the multi-temperature gas that one expects if there is substantial mass dropout (Peterson et al. 2001; Tamura et al. 2001; Böhringer et al 2001; Fabian et al. 2001; Molendi & Pizzolato 2001; Matsushita et al. 2002; Johnstone et al. 2002). The new data strongly suggest that mass dropout must be prevented by some additional source of heat that balances radiative losses. Two possibilities are currently being investigated for the heat source: (i) energy in jets, outflows, and radiation from a central active galactic nucleus (Pedlar et al. 1990; Tabor & Binney 1993; Churazov et al. 2000, 2002; Ciotti & Ostriker 2001; Brüggen & Kaiser 2002), and (ii) thermal energy from outer regions of the cluster transported to the central cooling gas by conduction (Narayan & Medvedev 2001, hereafter NM01).

Conduction in clusters has been discussed by several authors over the years in various contexts (Binney & Cowie, 1981; Tucker & Rosner, 1983; Bertschinger & Meiksin, 1986; Bregman & David, 1988; Gaetz, 1989; Rosner & Tucker, 1989; David et al., 1992; Pistinner & Shaviv, 1996), but its importance was always considered doubtful. For conduction to have any significant effect on the cooling gas in a cluster, the effective isotropic coefficient of conduction $`\kappa `$ has to be a reasonable fraction of the classical Spitzer (1962) conductivity $`\kappa _{Sp}`$. However, conventional wisdom is that magnetic fields strongly suppress conduction perpendicular to the field. Therefore, while conduction may be very efficient parallel to the field, the overall isotropic $`\kappa `$ is expected to be $`\kappa _{Sp}`$.

This picture has changed with the recent work of NM01 who, following earlier work by Rechester & Rosenbluth (1978), Chandran & Cowley (1998), Chandran et al. (1999) and Malyshkin & Kulsrud (2001), showed that a turbulent MHD medium in which the fluctuation spectrum extends over two or more decades of spatial scales has an effective $`\kappa `$ that is a substantial fraction of $`\kappa _{Sp}`$. This is because of chaotic transverse wandering of field lines as a result of strong MHD turbulence (Goldreich & Sridhar, 1995, 1997), which leads to a large enhancement of cross-field diffusion. NM01 estimated that the ratio $`f=\kappa /\kappa _{Sp}`$ in a turbulent MHD medium is approximately $`0.2`$, though, given the uncertainties in their model, the value could probably lie anywhere in the range $`0.10.4`$.

NM01 showed via a simple order-of-magnitude estimate that thermal conduction with $`f0.2`$ is sufficient to balance the radiative losses of the cooling gas in X-ray clusters. Based on this, they suggested that conduction might be a large part of the explanation for the lack of mass dropout in X-ray clusters. Their suggestion finds support in the work of Voigt et al. (2002) and Fabian et al. (2002) who show that conduction with $`f`$ in the range $`0.11`$ can indeed be a dominant heat source in X-ray clusters except perhaps in the very central regions.

The aim of the present paper is to study conduction in galaxy clusters in more detail. Our goals are two-fold. First, we consider a very simple model in which the gas is in hydrostatic equilibrium and in thermal balance, with cooling exactly compensated by heat conduction. We investigate whether this basic model can fit the observed shapes of the electron density profile $`n_e(r)`$ and temperature profile $`T(r)`$ of clusters for which high-resolution observations are available. We believe this approach, which is complementary to that of Voigt et al. (2002), is a good test of the conduction hypothesis, and a reliable method of estimating what value of $`f`$ is needed in various clusters to explain the observations. Second, having found a steady state model for a given cluster, we check whether the hot gas is thermally stable. We briefly discuss the importance and plausibility of other energy sources.

The paper is organized as follows. In §2 we write down the relevant differential equations and boundary conditions for a steady equilibrium model of hot gas in a cluster. In §3 we compare model predictions with data on 10 clusters, and in §4 we discuss the thermal stability of the gas in these clusters. We conclude with a discussion in §5. Throughout the paper, we take $`H_0=70`$ km sec-1 Mpc-1, $`\mathrm{\Omega }_M=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, rescaling observational results whenever the original papers have used a different cosmology for their analysis.

## 2 The Model

We consider a spherically symmetric cluster containing hot gas in equilibrium with no inflow. The condition of hydrostatic equilibrium gives

$$\frac{1}{\rho }\frac{dp}{dr}=\frac{d\varphi }{dr},$$
(1)
where $`\rho (r)`$ is the gas density, $`p(r)`$ is the pressure, and $`\varphi (r)`$ is the gravitational potential. We assume that gas pressure dominates and that magnetic fields have a negligible dynamical effect, so that

$$p=\frac{\rho k_BT}{\mu m_u}=\frac{\mu _e}{\mu }n_ek_BT\rho c_s^2,$$
(2)
where $`T(r)`$ is the temperature of the gas, $`n_e(r)`$ is the electron number density, $`\mu `$ is the mean molecular weight, $`\mu _e`$ is the molecular weight per electron, $`m_u`$ is the atomic mass unit, and $`c_s(r)`$ is the isothermal sound speed of the gas. We assume that the gas has hydrogen fraction $`X=0.7`$ and helium fraction $`Y=0.28`$, so that $`\mu =0.62`$, $`\mu _e=1.18`$ and $`\mu _e/\mu =1.91`$ (taking hydrogen and helium to be fully ionized).

The energy equation of the gas is a simple balance between radiative cooling and heating by thermal conduction. Thus, in steady state, we have

$$\frac{1}{r^2}\frac{d}{dr}(r^2F)=j,$$
(3)
where $`F(r)`$ is the heat flux due to electron conduction and $`j(r)`$ is the rate of radiative cooling per unit volume of the gas. We assume that the conductivity of the gas $`\kappa `$ is a fixed fraction $`f`$ of the classical Spitzer (1962) conductivity. Thus,

$$F=\kappa \frac{dT}{dr},$$
(4)
$$\kappa =f\kappa _{Sp}=f\frac{1.84\times 10^5T^{5/2}}{\mathrm{ln}\mathrm{\Lambda }}\mathrm{erg}\mathrm{s}^1\mathrm{K}^1\mathrm{cm}^1,$$
(5)
where $`\mathrm{ln}\mathrm{\Lambda }37`$ is the usual Coulomb logarithm. In the relevant range of temperatures the gas cools primarily by thermal bremsstrahlung. We therefore write (Rybicki & Lightman, 1979)

$$j=2.1\times 10^{27}n_e^2T^{1/2}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1,$$
(6)
where the numerical coefficient corresponds to the particular gas composition that we have assumed; the coefficient also includes an average Gaunt factor.

The gravitational acceleration $`d\varphi (r)/dr`$ due to the mass in the cluster is obtained by solving Poisson’s equation,

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d\varphi }{dr}\right)=4\pi G(\rho _{DM}+\rho ),$$
(7)
where $`\rho _{DM}(r)`$ is the density distribution of the dark matter, which we assume has the NFW form (Navarro et al., 1997) with a softened core:

$$\rho _{DM}(r)=\frac{M_0/2\pi }{(r+r_c)(r+r_s)^2}.$$
(8)
Here, $`r_s`$ is the standard scale radius of the NFW profile, $`r_c`$ is a core radius inside which the density becomes constant, and $`M_0`$ is a characteristic mass (e.g., if $`r_c=0`$, the mass inside a sphere of radius $`r_s`$ is $`(2\mathrm{ln}21)M_0`$).

We relate the parameters of the NFW profile to the temperature $`T_{out}`$ of the hydrostatic gas at $`rr_s`$, well outside the cooling region of the cluster. We use the mass-temperature relation from Afshordi & Cen (2002) and the mass-scale relation from Maoz et al. (1997), so that, once we know $`T_{out}`$, we can estimate the virial mass $`M_{vir}`$ and $`r_s`$. The former is defined as the mass inside the radius within which the mean density of the dark matter is $`200\rho _{crit}(z)`$, where $`\rho _{crit}(z)`$ is the critical density of the universe at the redshift of the cluster: $`\rho _{crit}(z)=3H(z)^2/8\pi G`$. The normalization of the NFW profile is then obtained via a simple calculation:

$$M_0=2\pi \rho _{crit}(z)r_s^3\left(\frac{200}{3}\right)\frac{c^3}{\mathrm{ln}(1+c)c/(1+c)},$$
(9)
where $`c`$ is the concentration parameter:

$$c=\frac{1}{r_s}\left(\frac{3M_{vir}}{4\pi 200\rho _{crit}(z)}\right)^{1/3}.$$
(10)
It is interesting that simulations of cluster formation allow us almost completely to determine the cluster potential from a single observed number, $`T_{out}`$. The only remaining uncertainty is the shape of the density distribution in the very inner regions of the cluster, which we parametrize through the core radius $`r_c`$.

The value of the core radius $`r_c`$ is somewhat uncertain. It is zero in the original NFW model, but cluster lensing studies (Tyson et al., 1998; Shapiro & Iliev, 2000) suggest that it may be tens of kpc in some clusters. As explained in §3 below, we try two values of $`r_c`$: $`r_c=0`$ and $`r_c=r_s/20`$.

Equations (1), (3), (4) and (7) are four ordinary differential equations for the four variables $`n_e(r)`$, $`T(r)`$, $`F(r)`$ and $`d\varphi (r)/dr`$. We solve these equations numerically between $`r=0`$ and $`r=2r_s`$ with the following four boundary conditions. Two boundary conditions are taken directly from the observations of a given cluster, namely the temperature $`T_{in}=T(0)`$ at the cluster center and $`T_{out}=T(2r_s)`$. For $`T_{out}`$ we use the observed temperature at around $`rr_s`$, which is well outside the cooling region in the center. The other two boundary conditions are $`F(0)=0`$ and $`M(0)=0`$, where $`M(r)`$ is the the mass (dark matter plus gas) inside radius $`r`$.

We note that we make a slight approximation by assuming a perfect steady state. In real clusters with conductive heating, the energy for the cooling ultimately comes from gravity and so the gas in the cluster must slowly settle in the gravitational potential. This effect would be accentuated by conduction to the outside, as discussed by Loeb (2002). We do not include gravitational settling in our model. Instead, we assume that the outer regions of the cluster behave like an infinite isothermal heat source at a temperature of $`T_{out}`$. We have experimented with different radii for the outer boundary of our models, from $`r=r_s`$ to $`r=5r_s`$, and find that the results are not sensitive to this choice. The models discussed below correspond to $`r=2r_s`$.

## 3 Results

We have analysed the following 10 clusters for which high angular resolution data are available from Chandra or XMM-Newton or both: A1795, A1835, A2052, A2199, A2390, A2597, Hydra A, RXJ 1347.5-1145, Ser 159-03, and 3C295. Table 1 lists the relevant parameters for each cluster. In the case of several of these clusters, deprojected density and temperature profiles are available, and we have used them in our modeling. In a few cases, however (A2597, A2052, Hydra A, Ser 159-03), the original authors have not published deprojected profiles, and we had to use the un-deprojected data. Given the observational uncertainties, deprojection only affects the profiles in the innermost regions of the cluster, within about 30 kpc. On the other hand, our model solutions only depend slightly on the inner temperature. Therefore, we believe that the shapes of the profiles are essentially the same outside about 30 kpc, whether we use deprojected or un-deprojected temperatures.

Table 1 lists for each cluster the inner and outer temperatures, $`T_{in}`$ and $`T_{out}`$, which were used as boundary conditions in the modeling. It also lists the parameters of the NFW profile of the dark matter, which were estimated from $`T_{out}`$ as described in §2. The last column gives references to the original publications for the X-ray data.

We carried out the modeling as follows. We first assumed that the core radius $`r_c=0`$ and selected a value of $`f`$ (say 0.2). We adjusted the inner density $`n_{in}=n_e(0)`$ of the gas such that on integrating the four differential equations (1), (3), (4) and (7) from $`r=0`$ to $`r=2r_s`$, the model temperature agreed with the measured outer temperature $`T_{out}`$. We then compared the model profiles of $`n_e(r)`$ and $`T(r)`$ with the data; we restricted ourselves to a comparison by eye since neither the deprojected data nor the models are reliable enough for a detailed $`\chi ^2`$ fit. If the fit was not satisfactory, we tried other values of $`f`$. If no value of $`f`$ gave a satisfactory fit, we repeated the process by setting $`r_c=r_s/20.`$ (The ratio 20 was selected to be consistent with estimates from lensing; the solutions are insensitive to the exact value.) For most clusters, we obtained a good fit for some choice of $`f`$ and one of the two values of $`r_c`$ we tried. The values of the best-fit parameters are listed in Table 2.

We consider a model to be acceptable if the model profiles of $`n_e(r)`$ and $`T(r)`$ agree well with the data and if the required value of $`f`$ lies within the range $`f0.10.4`$ predicted by theory (NM01). By this criterion, 5 of the 10 clusters we studied are consistent with a pure thermal conduction model: A1795, A1835, A2199, A2390, and RXJ 1347.5-1145. Figures 1, 2 and 3 show the fits for three of these clusters (A1795, A2199 and A2390).

We consider the good agreement between the model and the data encouraging since the model has very few adjustable parameters. Note in particular that our primary adjustable parameter $`f`$ affects only the overall normalization of $`n_e(r)`$, but has no effect on the shapes of either $`n_e(r)`$ or $`T(r)`$. The fact that the model curves agree well with the data is thus a significant result. To our knowledge, this kind of modeling of the shapes of profiles has not been reported before. There is a slight inconsistency in the density profile below 10 kpc in A2199 (and also RXJ 1347.5-1145), but this is hardly surprising since the central galaxy is likely to introduce large perturbations to our assumed NFW profile in this region.

Of the 5 clusters for which we obtain acceptable models, we note that A1835 is a complicated case because the Chandra and XMM-Newton data do not agree with each other; indeed, Markevitch (2002) has shown that there are problems with the analysis of both datasets. We obtain a reasonable model for this cluster if we use the density profile from XMM-Newton data (Majerowicz et al., 2002) corrected for flares. For A1795, although both Chandra and XMM-Newton data sets are available (Ettori et al., 2002; Tamura et al., 2001) we do not use the XMM-Newton temperature profile because it is un-deprojected.

The remaining 5 clusters, A2052, A2597, Hydra A, Ser 159-03 and 3C295, all require large values of $`f`$, generally greater than unity. Since such values are unphysical, we conclude that these clusters must have some other source of energy in addition to thermal conduction. Heating from a central active galactic nucleus is a promising possibility (Churazov et al., 2002; Brüggen & Kaiser, 2002). The idea is supported by the fact that all these clusters possess strong active galactic nuclei (AGN) in their centers which are also all powerful extended radio sources (3C317 in A2052, PKS2322-122 in A2597, 3C218 in Hydra A, PMNJ2313-4243 in Sersic 159-03, and the cluster 3C295 is named after the radio source in its center). Strong extended radio emission is indicative of powerful outflows interacting with the intracluster medium (e.g. Bridle & Perley 1984). In contrast, of the 5 well-fit clusters, only two (A1795 and A2199) are known to have strong extended radio sources in their centers (4C26.42 and 3C338, respectively).

Mass inflow and dropout is, of course, another possible source of energy, as has been discussed for many years in the context of cooling flows. However, as mentioned in §1, spectroscopic analyses of Chandra and XMM-Newton data of several clusters have failed to find any evidence for the multi-phase gas expected if there is significant mass inflow and dropout. The observations place an upper limit on the mass accretion rate in a few clusters of about 10-20% of the nominal value needed to explain the observed X-ray flux (e.g. David et al. 2001; Schmidt et al. 2001; Johnstone et al. 2002). Accretion at this level cannot substantially modify the energy balance, so one needs either conduction or mechanical energy from an AGN to fit the observations. We have computed models with larger levels of mass accretion, at the level of hundreds of $`M_{\mathrm{}}`$ yr-1, and we do find that accretion can then explain the X-ray data. Although the details of accretion and energy release during the mass dropout are uncertain (see different versions of energy balance in Bertschinger & Meiksin 1986 and Sarazin 1988), we find that we can reproduce the observed density and temperature profiles with reasonable assumptions.

## 4 Thermal Stability

To analyze the stability of the steady-state solutions described in §3, we work with the full time-dependent equations. Equation (1) generalizes to the Euler equation,

$$\frac{𝐯}{t}+(𝐯)𝐯=\frac{1}{\rho }p\varphi ,$$
(11)
where $`𝐯`$ is the gas velocity. We also have the continuity equation,

$$\frac{\rho }{t}+(\rho 𝐯)=0.$$
(12)
The energy equation (3) now becomes the time-dependent entropy equation,

$$\rho T\left[\frac{s}{t}+(𝐯)s\right]=j𝐅+q^+,$$
(13)
where the specific entropy $`s`$ of the gas is defined in terms of the specific internal energy $`ϵ`$ and specific enthalpy $`w`$ by

$$Tds=dϵ\frac{p}{\rho ^2}d\rho =dw\frac{dp}{\rho }.$$
(14)
For an ideal gas,

$$ϵ=\frac{1}{(\gamma 1)}c_s^2,w=\frac{\gamma }{(\gamma 1)}c_s^2,$$
(15)
where $`\gamma `$ is the adiabatic index and $`c_s`$ is the isothermal sound speed (Eq. 2). In our calculations, we set $`\gamma =5/3`$.

For completeness, we have included in equation (13) a heating term $`q^+`$ in addition to the term due to thermal conduction. This term may refer, for instance, to heating from a central AGN. We assume that the heating rate per unit volume has a power-law dependence on the density,

$$q^+=q_0^+\rho ^\alpha .$$
(16)
The heat flux equation (4) and the Poisson equation (7) are still valid, except that we generalize them to include non-spherical perturbations:

$$𝐅=\kappa T,$$
(17)
$$^2\varphi =4\pi G\left(\rho _{DM}+\rho \right).$$
(18)
The above five equations, (11), (12), (13), (17) and (18), describe the general time-dependent evolution of the gas. Note that, while we include the effect of a turbulent magnetic field on conduction through the parameter $`f<1`$ in equation (5), we do not include magnetic pressure or tension in the Euler equation (see Loewenstein 1990 for a discussion of these effects).

For the stability analysis, we take a steady state spherical solution of the equations, e.g., one of the cluster models described in §3, and consider local linear perturbations of the solution within the WKB approximation; that is, we assume that the perturbations are small and are proportional to $`\mathrm{exp}(ikxi\omega t)`$. Balbus (1988) has shown that radial perturbations, with wavevector $`k`$ parallel to the local direction of gravity, are not well-described by a local analysis; these perturbations also do not show thermal instability. In the following we therefore consider only tangential perturbations, with $`k`$ and the corresponding $`x`$ taken to be in a tangential direction with respect to the radius vector.

We carry out the linear stability analysis in the standard way; see Field (1965) for a classic treatment of the thermal stability problem. In brief, we linearize the equations, substitute perturbations of the form $`\mathrm{exp}(ikxi\omega t)`$ and derive the dispersion relation for the mode frequency $`\omega `$. If any of the solutions for $`\omega `$ has a positive imaginary part, the corresponding mode grows with time and the system is unstable. If all modes for all allowed values of $`k`$ have $`\mathrm{Im}(\omega )0`$, we say that the system is stable.

Before presenting the results, we note that there are several physical effects present in this problem, each with a characteristic frequency or time scale. In addition to the acoustic frequency $`c_sk`$, we have frequencies associated with gravity (Jeans frequency $`\omega _J`$), cooling ($`\omega _{cool}`$) and heating ($`\omega _{heat}`$), and a characteristic time scale associated with conduction ($`\tau _{cond}`$):

$$\omega _J=(4\pi G\rho )^{1/2},\omega _{cool}=\frac{j}{\rho c_s^2},\omega _{heat}=\frac{q^+}{\rho c_s^2},\tau _{cond}=\frac{\kappa T}{\rho c_s^4}.$$
(19)
Each of the parameters, $`\omega _J`$, $`\omega _{cool}`$, $`\omega _{heat}`$ and $`\tau _{cond}`$, has been defined in such a way that the parameter goes to 0 as the corresponding physical effect becomes negligible.

We discuss first a simple case in which we assume that conduction is negligible ($`\tau _{cond}0`$) and, for simplicity, we ignore self-gravity ($`\omega _J0`$). Energy balance of the equilibrium system requires $`q^+=j`$, and thus $`\omega _{heat}=\omega _{cool}`$. Assuming that $`q^+`$ and $`j`$ scale with density and temperature as in equations (16) and (6), and taking $`\gamma =5/3`$, we then obtain the following dispersion relation for $`\omega `$:

$$(i\omega )^3+\frac{1}{3}\omega _{cool}(i\omega )^2+\frac{5}{3}c_s^2k^2(i\omega )=\left(1\frac{2}{3}\alpha \right)\omega _{cool}c_s^2k^2.$$
(20)
By inspection we see that one of the solutions for $`\omega `$ is always positive imaginary if $`\alpha <1.5`$. Thus, if $`\alpha <1.5`$, the system is always thermally unstable. Physically, we expect many heating mechanisms to be described by $`\alpha `$ in the range between 0 (constant heating rate per unit volume) and 1 (constant heating rate per unit mass). According to equation (20), such systems will be thermally unstable. This suggests that models which try to solve the cooling flow problem by invoking local heating (Churazov et al., 2002; Brüggen & Kaiser, 2002) may still have some difficulties. While they allow equilibrium configurations, the equilibria they obtain are likely to be unstable. Unless the heating occurs preferentially in high density regions, with an index $`\alpha 1.5`$ (which seems implausible), these models require thermal conduction or some other agency for stability.

We move next to the problem of interest to this paper, namely models in which cooling is balanced entirely by heat conduction ($`q^+=0,\omega _{heat}=0`$). The dispersion relation now takes the form (again setting $`\gamma =5/3`$)

$`{\displaystyle \frac{3}{2}}\omega ^3+\omega ^2\left({\displaystyle \frac{1}{2}}i\omega _{cool}+ik^2c_s^2\tau _{cond}\right)+\omega \left({\displaystyle \frac{5}{2}}k^2c_s^2+{\displaystyle \frac{3}{2}}\omega _J^2\right)+`$
$`+{\displaystyle \frac{3}{2}}i\omega _{cool}k^2c_s^2+{\displaystyle \frac{1}{2}}i\omega _{cool}\omega _J^2ik^4c_s^4\tau _{cond}+ik^2c_s^2\tau _{cond}\omega _J^2=0`$ (21)
Analysis of this dispersion relation shows that an instability is present only when the condition

$$k^2c_s^2<\frac{1}{2}\left(\omega _J^2+\frac{3}{2}\frac{\omega _{cool}}{\tau _{cond}}+\sqrt{\omega _J^4+5\frac{\omega _J^2\omega _{cool}}{\tau _{cond}}+\frac{9}{4}\frac{\omega _{cool}^2}{\tau _{cond}^2}}\right)$$
(22)
is satisfied.

We see that, as expected (e.g., Medvedev & Narayan 2002), conduction eliminates the thermal instability for short wavelength modes. Physically, the reason is that the shorter the wavelength the larger the temperature gradient for a given mode amplitude; since conduction tries to iron out temperature fluctuations and the conductive flux is proportional to the temperature gradient, short wavelength modes are stabilized. For long wavelengths, however, sound waves can no longer stabilize the gravitational Jeans instability, and (finite) conduction cannot stabilize the thermal instability. This is the physical meaning of the upper limit on $`k`$ given in equation (22).

Let us denote the critical wavevector for instability by $`k_{crit}`$ and the corresponding wavelength by $`\lambda _{crit}=2\pi /k_{crit}`$. The relevant question is the following: can an unstable mode with wavelength $`\lambda _{crit}`$ or larger fit within the system? Recall that our modes are assumed to be in the tangential direction and we carried out the analysis under a local WKB approximation. For this analysis to be meaningful, the wavelength of the mode must be substantially smaller than the circumference $`2\pi r`$. A given cluster model is thus stable if $`\lambda _{crit}`$ is greater than the largest wavelength for which the WKB analysis is valid, which is probably around $`\pi r`$.

Figure 4 shows the variation of $`\lambda _{crit}/\pi r`$ with $`r`$ for the 5 clusters for which we obtained satisfactory models with conduction (§3). We see that all 5 clusters are thermally stable according to the criterion described above. One interesting feature is that all the clusters show very similar minimal values of $`\lambda _{crit}/\pi r`$ which appear to be fairly close to the threshold of instability. Such a situation would be natural if the clusters, to begin with, had more hot gas than they have now and were thermally unstable. One expects gas to have cooled and dropped out early in the life of the cluster so that the system naturally reached, and has since maintained, a marginally stable configuration. A global stability analysis would need to be carried out to check this idea in more detail.

## 5 Discussion

NM01 found that thermal conduction with $`f\kappa /\kappa _{Sp}`$ of order a few tenths can explain the gross energetics of hot gas in galaxy clusters. They showed that the heat flux transported to the cooler central regions of the cluster from further out is roughly consistent with what is needed to replace the energy lost through radiative cooling. Gruzinov (2002) and Fabian et al. (2002) confirmed this conclusion.

In this paper, we have carried out a further test of the conduction model by investigating the shapes of the density and temperature profiles, $`n_e(r)`$ and $`T(r)`$, in an equilibrium cluster. We solve for these profiles self-consistently, as compared to previous authors, e.g., Bertschinger & Meiksin (1986); Bregman & David (1988); Voigt et al. (2002) and Fabian et al. (2002), who either used the observations directly or employed simple analytical expressions for the shapes of $`n_e(r)`$ and/or $`T(r)`$.

Among the 10 clusters that we have studied, 5 clusters, namely, A1795, A1835, A2199, A2390 and RXJ1347.5-1145, can be fitted well with a pure conduction model and with values of $`f0.20.4`$ (Table 2, Figs. 1, 2, 3). Since the model involves no fitting parameters other than $`f`$, which is used primarily to set the normalization of $`n_e`$, and $`r_c`$, which is not really a fitting parameter but rather is given one of two values (§3), the good agreement between the model profiles of $`n_e(r)`$ and $`T(r)`$ and the data is very encouraging. For these five clusters, the best-fit values of $`f`$ lie within the range expected on theoretical grounds (NM01). We believe that these are strong arguments in support of the conduction model. We should note, however, that the theoretical estimate $`f0.10.4`$ obtained by NM01 is based on a specific (plausible) model of the magnetic field topology, suggested by the work of Goldreich & Sridhar (1995, 1997), but the answer might change substantially if the field topology (which is poorly understood in clusters) is very different.

The other 5 clusters, namely A2052, A2597, Hydra A, Ser 159-03 and 3C295, require larger values of $`f`$, of order unity or greater, and we consider them to be inconsistent with a pure conduction model. It is interesting that these five clusters all have active nuclei with strong radio activity, so that nuclear activity might well supply the additional heat required by our models. It is unclear exactly how an AGN would heat the cluster medium. Turbulent mixing of hot gas in a jet with the surrounding cooler gas is a possibility (Churazov et al., 2002; Brüggen & Kaiser, 2002), though the effect of magnetic fields on this mixing is poorly understood. Heating of the gas by cosmic rays (Loewenstein et al., 1991) or hard X-rays (Ciotti & Ostriker, 2001) from the AGN are other possibilities.

In §4 we showed by means of a local linear stability analysis that the five clusters for which we have obtained good fits with a pure conduction model are all thermally stable. This is an interesting result because, without conduction, cooling gas in clusters is generally thermally unstable. We find it reassuring that a single mechanism, namely conduction, is able both to supply the heat lost via radiative cooling and to control the thermal instability in these clusters. In the other five clusters, which require an additional source of energy such as a central AGN, the thermal stability of the gas is not assured. One way of ensuring stability is to have a heat source with a strong non-linear dependence on density, e.g., $`\alpha >1.5`$ in equation (16); this appears somewhat unnatural. Alternatively, conduction, though not an important energy source in these clusters, might still play a role in controlling the instability.

Based on the above results, we suggest that AGN heating and conduction both are important in clusters. In some clusters, AGN heating dominates and conduction plays a secondary, though still important, role in helping to stabilize the system. In other clusters, AGN activity is weak, and conduction supplies the energy for cooling as well as prevents the gas from becoming unstable.

While our work suggests that conduction is effective over much of the volume of a cluster, we note that Chandra observations have revealed sharp temperature jumps in some clusters (e.g., Markevitch et al. 2000; Vikhlinin, Markevitch & Murray 2001ab). These cold fronts are associated with the interface between the intracluster medium and an infalling galaxy or group. The observations clearly indicate that conduction is strongly suppressed across the fronts. Vikhlinin et al. (2001b) have argued that the magnetic field is stretched out parallel to the cold front by the motion of the infalling galaxy/group and that this explains the very low conductivity. Mazzotta et al. (2002) present interesting evidence for a possible Kelvin-Helmholtz instability in the cold front in A3667, confirming that the fronts are probably transient features with lifetimes of only a couple of dynamical times. The presence of these fronts is thus unlikely to affect the strong conduction that we have hypothesized over the bulk of the cluster. If cold fronts survive for longer than the cooling time of the gas, then the thermally isolated cool gas would experience runaway cooling. Since such runaway cooling is ruled out by the absence of substantial mass dropout, the lifetimes of cold fronts cannot be longer than the cooling time.

The authors thank Bill Forman, Christine Jones Forman and Larry David for useful discussions, and the referee for helpful comments on the manuscript. RN thanks the National Science Foundation for support under grant AST 9820686.

Table 1.

Observational and cosmological parameters of the clusters.

| Name | $`T_{in}`$ (keV) | $`T_{out}`$ (keV) | $`M_{vir}`$ ($`M_{\mathrm{}}`$) | $`r_s`$ (kpc) | $`M_0`$ ($`M_{\mathrm{}}`$) | ref |
| --- | --- | --- | --- | --- | --- | --- |
| Abell 1795 | 2 | 7.5 | $`1.2\times 10^{15}`$ | 460 | $`6.6\times 10^{14}`$ | (1,2) |
| Abell 1835 | 5 | 11 | $`1.7\times 10^{15}`$ | 520 | $`1.1\times 10^{15}`$ | (3) |
|  | 4 | 8 | $`1.1\times 10^{15}`$ | 470 | $`8.8\times 10^{14}`$ | (4,5) |
| Abell 2052 | 1.3 | 3.5 | $`3.9\times 10^{14}`$ | 340 | $`2.3\times 10^{14}`$ | (6) |
| Abell 2199 | 1.6 | 5 | $`6.6\times 10^{14}`$ | 390 | $`3.8\times 10^{14}`$ | (7) |
| Abell 2390 | 4 | 11 | $`1.7\times 10^{15}`$ | 520 | $`1.1\times 10^{15}`$ | (8) |
| Abell 2597 | 1 | 4 | $`4.6\times 10^{14}`$ | 360 | $`2.8\times 10^{14}`$ | (9) |
| Hydra A | 3 | 4.1 | $`4.9\times 10^{14}`$ | 370 | $`2.9\times 10^{14}`$ | (10) |
| RXJ 1347.5-1145 | 6 | 16 | $`2.5\times 10^{15}`$ | 600 | $`1.8\times 10^{15}`$ | (11) |
| Sersic 159-03 | 2 | 2.7 | $`2.6\times 10^{14}`$ | 310 | $`1.6\times 10^{14}`$ | (12) |
| 3C295 | 2 | 6 | $`5.8\times 10^{14}`$ | 420 | $`4.9\times 10^{14}`$ | (13) |

References: (1) Ettori et al. (2002); (2) Tamura et al. (2001); (3) Schmidt et al. (2001); (4) Markevitch (2002); (5) Majerowicz et al. (2002); (6) Blanton et al. (2001); (7) Johnstone et al. (2002); (8) Allen et al. (2001b); (9) McNamara et al (2001); (10) David et al. (2001); (11) Allen et al. (2002); (12) Kaastra et al. (2001); (13) Allen et al. (2001a).

Notes: Observed parameters ($`T_{in}`$ is the temperature in the innermost region of the cluster, $`T_{out}`$ is the temperature outside the cooling region) and inferred parameters of the NFW profile (the virial mass of the cluster, $`M_{vir}`$; the scale radius $`r_s`$ and the characteristic mass $`M_0`$); see §2 on how the NFW parameters were obtained. The last column contains references to the original publications for the X-ray data on the gas temperature and electron density profiles.

Table 2.

Best-fit parameters of the conduction model.

| Name | f (best-fit) | $`n_e`$(0) (cm-3) | $`r_c`$ (kpc) |  |
| --- | --- | --- | --- | --- |
| Abell 1795 | 0.2 | 0.049 | $`r_s/20`$ | (see Fig. 1) |
| Abell 1835 | 0.4 | 0.17 | $`r_s/20`$ |  |
| Abell 2052 | 0.6 | 0.089 | 0 |  |
| Abell 2199 | 0.4 | 0.074 | 0 | (see Fig. 2) |
| Abell 2390 | 0.3 | 0.069 | $`r_s/20`$ | (see Fig. 3) |
| Abell 2597 | 2.4 | 0.13 | $`r_s/20`$ |  |
| Hydra A | 1.5 | 0.020 | 0 |  |
| RXJ 1347.5-1145 | 0.3 | 0.11 | $`r_s/20`$ |  |
| Sersic 159 | 5.6 | 0.042 | 0 |  |
| 3C295 | 1 | 0.084 | $`r_s/20`$ |  |

Notes: $`f=\kappa /\kappa _{Sp}`$ is the value of the conductivity coefficient required to reproduce the observed electron density, $`n_e`$(0) is the corresponding model electron density in the center, and $`r_c`$ is the core radius adjusted to fit the observed temperature profile, for which we only consider two values (see §3).


text:
# Iterating Random Functions on a Finite Set

## 1 Introduction

Let $`f_1,f_2,f_3,\mathrm{}`$ be a sequence of functions chosen independently and uniformly randomly from the $`n^n`$ functions on $`[n]`$. Let $`g_1=f_1`$, and for $`t>1`$ let $`g_t=f_tg_{t1}`$ be the composition of the first $`t`$ random functions. Define $`T(f_i_{i=1}^{\mathrm{}})`$ to be the smallest $`t`$ for which $`g_t`$ is a constant function. (i.e. $`g_t(i)=g_t(j)`$ for all $`ij.`$) This manuscript contains a simple derivation of the asymptotic distribution of $`T`$. We had originally intended to publish it in a journal, but we recently learned that the the asymptotic distribution of $`T`$ is “well known”. It was apparently known to Kingman twenty years ago ,,, and stronger results are fully proved in Donnelly. There is a lot of related work by Mőhle and others,e.g. ,, and .

For $`m>1`$, let $`\tau _m=\left|\{t:|Range(g_t)|=m\}\right|`$ be the the number of iterates for which the range has exactly $`m`$ elements. Thus $`T=_{m=2}^n\tau _m`$. The random variables $`\{\tau _m\}_{m=2}^n`$ are not independent. They are however conditionally  independent once we specify the set of visited states. Fortunately this set is well behaved, has some convenient properties that enable us to do computations.

Let $`\xi =\mathrm{log}\mathrm{log}n`$, and decompose $`T`$ as $`T=T_1+T_2`$, where $`T_1=\underset{m=2}{\overset{\xi }{}}\tau _m`$ and $`T_2=\underset{m=\xi +1}{\overset{n}{}}\tau _m`$. Let $`𝒜=\underset{m=1}{\overset{\xi }{}}[\tau _m>0]`$. The following facts from  will be needed (See also Theorem 5 of ):

###### Theorem 1

$`\mathrm{Pr}(𝒜)=1o(1),`$ and $`E(T_2)=o(n).`$

## 2 Characteristic Function

Let $`\lambda _k=\underset{j=1}{\overset{k1}{}}(1\frac{j}{n}).`$ Then we have

###### Theorem 2

$`E(e^{itT_1}|𝒜)=e^{it(\xi 1)}\underset{k=2}{\overset{\xi }{}}\frac{(1\lambda _k)}{1\lambda _ke^{it}}.`$

Proof:

Suppose $`g_{t1}`$ has an $`m`$ element range $`R=\{r_1,r_2,\mathrm{},r_m\}`$. What is the chance that the next function $`g_t`$ still has an $`m`$ element range? On $`R`$ we have $`n`$ choices for $`f_t(r_1)`$, then $`n1`$ choices for $`f_t(r_2)`$ etc. For $`xR`$, $`f_t(x)`$ can be chosen arbitrarily. Hence the number of functions $`f_t`$ for which $`g_t=f_tg_{t1}`$ has an an $`m`$ element range is $`n^{nm}\underset{j=0}{\overset{m1}{}}(nj).`$ Hence

$$\mathrm{Pr}(\tau _m=k|\tau _m>0)=\lambda _m^{k1}(1\lambda _m),$$
(1)
and consequently

$$E(e^{it\tau _m}|\tau _m>0)=\underset{k=1}{\overset{\mathrm{}}{}}\lambda _m^{k1}(1\lambda _m)e^{ikt}=\frac{(1\lambda _m)e^{it}}{1\lambda _me^{it}}$$

Now let $`\varphi _n(t)=E(e^{itT_1/n}|𝒜)`$ be the characteristic function of the normalized random variable $`T_1/n`$ on $`𝒜.`$ Then the following corollary follows immediately from Theorem 2.

###### Corollary 3

$`\varphi _n(t)=\underset{m=2}{\overset{\xi }{}}\frac{(1\lambda _m)e^{it/n}}{1\lambda _me^{it/n}}=e^{it(\xi 1)/n}\underset{m=2}{\overset{\xi }{}}\frac{(1\lambda _m)}{1\lambda _me^{it/n}}`$

###### Lemma 4

$`\varphi _n(t)=\underset{m=2}{\overset{\mathrm{}}{}}\frac{\left(\genfrac{}{}{0pt}{}{m}{2}\right)}{\left(\genfrac{}{}{0pt}{}{m}{2}\right)it}+o(1).`$

Proof: Note that, for $`m\xi `$,

$$1\lambda _m=\frac{1}{n}\left(\genfrac{}{}{0pt}{}{m}{2}\right)+O(\frac{\xi ^4}{n^2})$$
(2)
and

$$1\lambda _me^{it/n}=\frac{1}{n}\left(\left(\genfrac{}{}{0pt}{}{m}{2}\right)it\right)+O(\frac{\xi ^4}{n^2})$$
(3)
Therefore

$$\frac{(1\lambda _m)e^{it/n}}{1\lambda _me^{it/n}}=\frac{\left(\genfrac{}{}{0pt}{}{m}{2}\right)+O(\frac{\xi ^4}{n})}{\left(\genfrac{}{}{0pt}{}{m}{2}\right)it+O(\frac{\xi ^4}{n})}=\frac{\left(\genfrac{}{}{0pt}{}{m}{2}\right)}{\left(\genfrac{}{}{0pt}{}{m}{2}\right)it}\left(1+O(\frac{\xi ^4}{n})\right).$$
(4)
Therefore

$$\varphi _n(t)=\left(1+O(\frac{\xi ^4}{n})\right)^\xi \underset{m=2}{\overset{\xi }{}}\frac{\left(\genfrac{}{}{0pt}{}{m}{2}\right)}{\left(\genfrac{}{}{0pt}{}{m}{2}\right)it}$$

$$=(1+o(1))\underset{m=2}{\overset{\xi }{}}\frac{\left(\genfrac{}{}{0pt}{}{m}{2}\right)}{\left(\genfrac{}{}{0pt}{}{m}{2}\right)it}.$$

Finally, note that the infinite product $`\underset{m=2}{\overset{\mathrm{}}{}}\frac{\left(\genfrac{}{}{0pt}{}{m}{2}\right)}{\left(\genfrac{}{}{0pt}{}{m}{2}\right)it}=\underset{m=2}{\overset{\mathrm{}}{}}\frac{1}{1\frac{it}{\left(\genfrac{}{}{0pt}{}{m}{2}\right)}}`$ converges since $`\left(\genfrac{}{}{0pt}{}{m}{2}\right)^1`$ is convergent.

## 3 Simplification

To facilitate inversion, we reexpress the characteristic function $`\varphi `$ in a more convenient form. Working with the reciprocal, we have

$$\frac{1}{\varphi _n(t)+o(1)}=\underset{k1}{}\left(1\frac{2it}{k(k+1)}\right)=\underset{k1}{}\frac{(k\alpha )(k\beta )}{k(k+1)},$$
(5)
where $`\alpha =\frac{1\sqrt{1+8it}}{2},\beta =\frac{1+\sqrt{1+8it}}{2}.`$ It is well known  that

$$\frac{1}{\mathrm{\Gamma }(z)}=ze^{\gamma z}\underset{n1}{}\left((1+\frac{z}{n})e^{z/n}\right).$$

Since $`\alpha +\beta =1`$, the right side of equation (5) becomes

$$\underset{k1}{}\frac{(1\alpha /k)e^{\alpha /k}(1\beta /k)e^{\beta /k}}{(1+\frac{1}{k})e^{1/k}}=\frac{1}{\alpha \beta \mathrm{\Gamma }(\alpha )\mathrm{\Gamma }(\beta )}=\frac{\mathrm{cos}(\frac{\pi }{2}\sqrt{1+8it})}{2\pi it}.$$

Hence

$$\varphi _n(t)=\frac{2\pi it}{\mathrm{cos}(\frac{\pi }{2}\sqrt{1+8it})}+o(1).$$
(6)
## 4 Fourier Inversion

Inverting, we get the conditional density function $`f_n`$ for $`T_1/n`$:

$$f_n(x)=\frac{1}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}e^{itx}(\frac{2\pi it}{\mathrm{cos}(\frac{\pi }{2}\sqrt{1+8it})}+o(1))𝑑t.$$
(7)
Note that $`\frac{2\pi it}{\mathrm{cos}(\frac{\pi }{2}\sqrt{1+8it})}`$ has simple poles at $`t=i\left(\genfrac{}{}{0pt}{}{k}{2}\right),`$ for $`k=2,3,4,\mathrm{}`$ Since the residue at $`t=i\left(\genfrac{}{}{0pt}{}{k}{2}\right)`$ is $`i(1)^k(2k1)\left(\genfrac{}{}{0pt}{}{k}{2}\right)e^{\left(\genfrac{}{}{0pt}{}{k}{2}\right)x}`$, contour integration yields, for $`x>0,`$ $`f_n(x)=f(x)+o(1)`$ where

$$f(x)=\underset{k2}{}(1)^ke^{\left(\genfrac{}{}{0pt}{}{k}{2}\right)x}\left(\genfrac{}{}{0pt}{}{k}{2}\right)(2k1).$$

## 5 Main Result

For $`x>0`$, let $`F(x)=\underset{0}{\overset{x}{}}f(t)𝑑t.`$ Our main result is

###### Theorem 5

For any $`x>0`$, $`\underset{n\mathrm{}}{lim}\mathrm{Pr}(T/nx)=F(x).`$

Proof: For any $`x`$,

$$\mathrm{Pr}(T/nx)\mathrm{Pr}(T_1/nx)$$

$$=\mathrm{Pr}(T_1/nx|𝒜)\mathrm{Pr}(𝒜)+\mathrm{Pr}(T_1/nx|𝒜^c)\mathrm{Pr}(𝒜^c)$$

$$=\mathrm{Pr}(T_1/nx|𝒜)(1+o(1))+o(1)$$

$$=F(x)(1+o(1))+o(1).$$

In the other direction, let $`ϵ`$ be a fixed but arbitrarily small positive number. Then

$$\mathrm{Pr}(T/nx)\mathrm{Pr}(T_1/nxϵ\mathrm{and}T_2/nϵ)$$

$$\mathrm{Pr}(T_1/nxϵ)\mathrm{Pr}(T_2/n>ϵ)$$

$$\mathrm{Pr}(T_1/nxϵ)\frac{E(T_2/n)}{ϵ}$$

$$=F(xϵ)+o(1).$$

The theorem follows from this and the fact that $`F`$ is continuous.

## 6 Discussion

Although the ultimate behaviour of our chain is like Kingman’s coalescent , there are differences. In that process every state is visited, whereas in our process few of the high numbered states are visited. Let $`N=\underset{m=2}{\overset{n}{}}I_{[\tau _m>0]},`$ the number of states visited. In an earlier version of this manuscript, we conjectured that $`E(N)\sqrt{2\pi n}.`$ Robin Pemantle recently proved our conjecture and the corresponding central limit theorem. He may also be able to prove stronger and more general versions of this result, e.g. a functional limit theorem. .

Acknowledgement There is a large literature in applied probability that can be traced back to Kingman’s work, and it was not immediately obvious to us what is relevant. We are grateful to Simon Tavare for pointing our way to the the work of Donnelly and Mőhle.