text:
# 𝑅-modes of neutron stars with a solid crust

## 1 Introduction

Since the discovery of the gravitational radiation driven instability of the $`r`$-modes by Andersson (1998) and Friedman & Morsink (1998), a large number of studies on the properties of $`r`$-modes and inertial modes of rotating stars have been done to prove their possible importance in astrophysics (for a recent review see, e.g., Friedman & Lockitch 1999). Although early investigations of the $`r`$-mode instability were mainly applied to young hot neutron stars (e.g., Lindblom, Owen, & Morsink 1998; Owen et al. 1998), recently some authors have also discussed possible roles of the $`r`$-mode instability in old and cool neutron stars with a solid crust and a magnetic field, as those found in low mass X–ray binaries (LMXBs) (e.g., Andersson, Kokkotas, & Stergioulas 1999; Bildsten & Ushomirsky 2000; Rezzolla, Lamb, & Shapiro 2000). For example, Bildsten & Ushomirsky (2000) have suggested that the $`r`$-mode instability can be largely weakened by the effects of viscous damping in the boundary layer at the interface between the solid crust and the fluid core. Rezzolla et al. (2000) have also suggested that amplitude of $`r`$-mode oscillations tends to reduce due to the coupling to the magnetic field of a star.

It is well known that neutron stars with a solid crust show a rich spectrum of nonradial oscillation modes, for example, nonradial modes associated with the solid crust and those associated with the fluid core and the surface fluid ocean (McDermott et al. 1988). If the effects of rotation are taken into account, it is quite common to find avoided crossings, as functions of $`\mathrm{\Omega }`$, between various oscillation modes of the neutron star model with a solid crust (Lee & Strohmayer 1996). This is also the case for $`r`$-modes. When the angular rotation frequency $`\mathrm{\Omega }`$ of the star is small, the $`r`$-mode frequency observed in the corotating frame may be given by

$$\omega _r(\mathrm{\Omega })\frac{2m\mathrm{\Omega }}{l^{}(l^{}+1)},$$
(1)
where $`l^{}`$ and $`m`$ are the indices of the spherical harmonic function representing the dominant toroidal component of the displacement vector. On the other hand, the crustal toroidal mode frequency may be given by (e.g., Strohmayer 1991)

$$\omega _t(\mathrm{\Omega })\omega _t(0)+\frac{m\mathrm{\Omega }}{l^{}(l^{}+1)},$$
(2)
where $`\omega _t(0)`$ is the oscillation frequency of the mode at $`\mathrm{\Omega }=0`$. The frequencies of the two modes get close to each other at

$$\mathrm{\Omega }_{cross}\frac{l^{}(l^{}+1)}{m}\omega _t(0).$$
(3)
As shown by McDermott et al. (1988) and Lee & Strohmayer (1996), the oscillation frequency of the crustal toroidal mode that has no radial nodes of the eigenfunction is of order of $`\omega _t(0)/\sqrt{GM/R^3}10^2`$ for $`l^{}=2`$, and hence we have $`\mathrm{\Omega }_{cross}/\sqrt{GM/R^3}10^2\times (l^{}+1)`$ for the modes with $`l^{}=m`$, where $`M`$ and $`R`$ are the mass and radius of the neutron star, and $`G`$ is the gravitational constant. This suggests that the mode crossing can happen between the two modes at sufficiently small rotation rates. Almost all previous discussions on the $`r`$-mode instability, however, did not take account of the effects of mode crossing between the two different kinds of modes.111After our submission of this paper, a paper by Levin & Ushomirsky (2000) concerning a similar problem appeared in a preprint server “astro-ph”. They treated $`r`$-modes in a constant density neutron star with a constant shear modulus crust. In fact, the discussions have been based on the assumptions that the modal property of the $`r`$-modes is almost independent of the rotation frequency $`\mathrm{\Omega }`$ and hence that an extrapolation formula, obtained in the limit of $`\mathrm{\Omega }0`$, is applicable to rapidly rotating neutron stars with $`\mathrm{\Omega }/\sqrt{GM/R^3}1`$. Obviously these two assumptions are not necessarily correct for the $`r`$-modes that experience mode coupling with other kinds of oscillation modes.

In this paper, we examine the $`r`$-mode instability of neutron stars that have a solid crust, by considering the effects of the mode crossing between the $`r`$-modes and the crustal toroidal modes. The plan of this paper is as follows. In §2, we show the modal properties of the $`r`$-modes of neutron stars with a solid crust. In §3, we discuss the dissipation timescales of the $`r`$-modes. §4 is for discussions and conclusions.

## 2 $`r`$-Modes of Neutron Stars with a Solid Crust

Neutron star models we use in this paper are the same as those used in the modal analysis by McDermott et al. (1988). The models are taken from the evolutionary sequences for cooling neutron stars calculated by Richardson et al. (1982), where the envelope structure is constructed by following Gudmundsson, Pethick & Epstein (1983). These models are composed of a fluid core, a solid crust and a surface fluid ocean, and the interior temperature is finite and is not constant as a function of the radial distance $`r`$. The models are not isentropic and the Schwarzschild discriminant $`A`$ has finite values in the interior (see, e.g., Yoshida & Lee 2000b).

The method of calculation of nonradial oscillations of rotating neutron stars with a solid crust is the same as that used by Lee & Strohmayer (1996), in which the Cowling approximation is employed. The terms due to the Coriolis force are included in the perturbation equations but no effects of the centrifugal force are considered. This approximation may be justified for low frequency modes satisfying $`\left|2\mathrm{\Omega }/\omega \right|1`$ and $`\mathrm{\Omega }^2/(GM/R^3)1`$, where $`\omega `$ is the oscillation frequency observed in the corotating frame of the star (Unno et al. 1989). The eigenfunctions are expanded in terms of spherical harmonic functions $`Y_l^m(\theta ,\phi )`$ with different values of $`l`$ for a given $`m`$. Here spherical polar coordinates $`(r,\theta ,\phi )`$ are used. For example, the Lagrangian displacement vector, $`\xi _i`$ is expanded as

$$\xi _r=r\underset{l|m|}{\overset{\mathrm{}}{}}S_l(r)Y_l^m(\theta ,\phi )e^{i\sigma t},$$
(4)
$$\xi _\theta =r\underset{l,l^{}|m|}{\overset{\mathrm{}}{}}\left\{H_l(r)\frac{Y_l^m}{\theta }+T_l^{}(r)\frac{1}{\mathrm{sin}\theta }\frac{Y_l^{}^m}{\phi }\right\}e^{i\sigma t},$$
(5)
$$\xi _\phi =r\underset{l,l^{}|m|}{\overset{\mathrm{}}{}}\left\{H_l(r)\frac{1}{\mathrm{sin}\theta }\frac{Y_l^m}{\phi }T_l^{}(r)\frac{Y_l^{}^m}{\theta }\right\}e^{i\sigma t},$$
(6)
where $`l=|m|+2k`$ and $`l^{}=l+1`$ for even modes and $`l=|m|+2k+1`$ and $`l^{}=l1`$ for odd modes where $`k=0,1,2\mathrm{}`$ (Lee & Strohmayer 1996, see also Lee & Saio 1986). The symbol $`\sigma `$ denotes the oscillation frequency observed in an inertial frame, and we have $`\omega =\sigma +m\mathrm{\Omega }`$. In this paper, to obtain a good convergence of the eigenfunctions and eigenfrequencies at a given $`\mathrm{\Omega }`$, we keep an enough number of terms in the series expansion of the eigenfunctions.

We computed frequency spectra of $`r`$-modes, inertial modes, and crustal toroidal modes for the three-component neutron star models called NS05T7, NS05T8, and NS13T8 (see, McDermott et al. 1988). The obtained mode spectra for the three models are qualitatively the same. In this paper, we therefore show the mode spectrum for the model NS13T8, because this is the most compact model among the three. The mass $`M`$ and the radius $`R`$ of the model NS13T8 are $`M=1.326M_{\mathrm{}}`$ and $`R=7.853\mathrm{km}`$, respectively, and the central temperature is $`T_c=1.05\times 10^8\mathrm{K}`$.

In Figure 1, scaled eigenfrequencies $`\kappa \omega /\mathrm{\Omega }`$ of $`r`$-modes, inertial modes, and a crustal toroidal mode of the model NS13T8 are given as functions of $`\overline{\mathrm{\Omega }}\mathrm{\Omega }/\sqrt{GM/R^3}`$ for $`m=2`$, where $`\omega `$ is the frequency observed in the corotating frame of the star. Here only the fundamental $`r`$-modes with $`l^{}=m=2`$ are considered since they are most important for the $`r`$-mode instability of neutron stars (see Lockitch & Friedman 1999; Yoshida & Lee 2000a, 2000b). The crustal toroidal mode shown in Figure 1 is the fundamental mode with no radial nodes of the eigenfunction and has the smallest oscillation frequency of the mode of this kind for a given $`l^{}`$ at $`\mathrm{\Omega }=0`$ (see McDermott et al. 1988). Note that at $`\overline{\mathrm{\Omega }}0`$ it is practically impossible to correctly calculate rotational modes because of their coupling with high overtone $`g`$-modes. In this figure, we have used the notation given by $`{}_{l^{}}{}^{}r_{n}^{c}`$ and $`{}_{l^{}}{}^{}r_{n}^{s}`$ for the $`r`$-modes, $`{}_{l^{}}{}^{}t_{n}^{}`$ for the crustal toroidal modes, and $`{}_{l_0}{}^{}i`$ for the inertial modes, where $`l^{}`$ denotes the index of the spherical harmonic function associated with the dominant toroidal component of the displacement vector, $`n`$ is the number of radial nodes of the eigenfunction, and $`l_0`$ is the number introduced to classify the inertial modes (see, e.g., Yoshida & Lee 2000a). Note that the mode classification has been done at sufficiently small values of $`\overline{\mathrm{\Omega }}`$. Since there appear surface $`r`$-modes and core $`r`$-modes for the three-component neutron star models, we have also introduced the superscripts $`s`$ and $`c`$ to distinguish the two kinds of $`r`$-modes. The oscillation energy of the core $`r`$-mode is predominantly confined in the fluid core, and that of the surface $`r`$-mode is confined in the surface ocean. A typical eigenfunction of a surface $`r`$-mode is given in Lee & Strohmayer (1996). The modes in the figure are all odd modes (with $`l^{}=|m|`$ for the toroidal modes) and they are retrograde modes propagating in the opposite direction to that of stellar rotation.

Figure 1 clearly shows that the mode crossing between the $`{}_{2}{}^{}r_{0}^{c}`$-mode and the $`{}_{2}{}^{}t_{0}^{}`$-mode, as predicted by equations (1) to (3), leads to an avoided crossing of their frequencies as functions of $`\overline{\mathrm{\Omega }}`$. At the avoided crossing, the modal properties of the two modes are exchanged. This is understood by examining the eigenfunctions of the two modes before and after the crossing, which are indicated by the filled circles in Figure 1, where the labels (a), (c), and (d) attaching to the filled circles are corresponding to the panels (a), (c), and (d) in Figure 2. Note that the $`{}_{2}{}^{}t_{0}^{}`$-mode at $`\overline{\mathrm{\Omega }}=4\times 10^3`$, which should carry the label (b), is not shown in Figure 1, because $`\kappa 1`$. In Figure 2, the expansion coefficients $`S_3`$, $`H_3`$, and $`iT_2`$ of the two modes with $`l^{}=m=2`$ at $`\overline{\mathrm{\Omega }}=4\times 10^3`$ in panels (a) and (b) and at $`\overline{\mathrm{\Omega }}=4\times 10^2`$ in panels (c) and (d) are given as functions of the fractional radius $`r/R`$, where the eigenfunctions are normalized as $`iT_l^{}(r_{bc})=1`$ with $`r_{bc}`$ being the bottom of the solid crust of the model. When $`\overline{\mathrm{\Omega }}`$ is sufficiently small, the $`{}_{2}{}^{}r_{0}^{c}`$-mode in panel (a) and the $`{}_{2}{}^{}t_{0}^{}`$-mode in panel (b) have their typical eigenfunctions and the influence of other modes on the eigenfunction is negligible. Through the avoided crossing, however, the modal properties of the two modes are exchanged along the frequency curves as $`\overline{\mathrm{\Omega }}`$ increases. As shown by Figure 2, the mode in panel (d), which is corresponding to the mode (d) on the frequency curve of $`{}_{2}{}^{}t_{0}^{}`$ in Figure 1, has in the fluid core the eigenfunction characteristic of the $`{}_{2}{}^{}r_{0}^{c}`$-mode. On the other hand, the mode in panel (c), which corresponds to the mode (c) on the frequency curve of $`{}_{2}{}^{}r_{0}^{c}`$, has in the crust the eigenfunction characteristic of the $`{}_{2}{}^{}t_{0}^{}`$-mode. However, this avoided crossing between the two modes is not sharp in the sense that the eigenfunction of each of the two modes suffers the contamination of the other mode after the avoided crossing. For example, the mode (d), which may be regarded as a core $`r`$-mode, has also oscillation amplitudes in the solid crust. This is in a sharp contrast with avoided crossings, for example, between $`{}_{2}{}^{}r_{0}^{s}`$ and $`{}_{2}{}^{}r_{0}^{c}`$ and between $`{}_{5}{}^{}i`$ and $`{}_{2}{}^{}r_{0}^{c}`$, for which the mode properties are almost completely exchanged through the crossings. Because of the avoided crossing between the $`{}_{2}{}^{}r_{0}^{c}`$-mode and the crustal toroidal modes, the $`{}_{2}{}^{}r_{0}^{c}`$-mode will lose its simple modal property of the eigenfrequency and the eigenfunction as a function of $`\overline{\mathrm{\Omega }}`$. We note that the fundamental $`r`$-mode with $`l^{}=m`$ of fluid stars has such a simple modal property as a function of $`\overline{\mathrm{\Omega }}`$ (e.g., Yoshida & Lee 2000b; Karino et al. 2000).

Avoided crossing between the $`{}_{l^{}}{}^{}r_{0}^{c}`$ and $`{}_{l^{}}{}^{}t_{0}^{}`$ modes with $`l^{}=m`$ may be common for neutron stars with a solid crust. To show this, we calculate the $`{}_{l^{}}{}^{}r_{0}^{c}`$ modes with $`l^{}=m=2`$ for $`n=1`$ polytropic neutron star models in which a solid crust is artificially embedded and isentropic structure is assumed. We assume the same physical parameters, such as the mass, radius, and the crust thickness, as those of the model NS13T8. For the shear modulus $`\mu _0`$, we assume that the ratio $`\mu _0/\rho ^{4/3}`$ is constant in the crust where $`\rho `$ is the mass density (see e.g., McDermott et al. 1988). In Table 1, we tabulated the values of $`\mu _0/\rho ^{4/3}`$ in the crust and the frequencies $`\omega _t(0)`$ of the fundamental $`{}_{2}{}^{}t_{0}^{}`$ mode at $`\overline{\mathrm{\Omega }}=0`$. We note that $`\omega _t(0)`$ of the polytropic model is approximately proportional to the value of $`\sqrt{\mu _0/\rho ^{4/3}}`$. Figure 3 gives the scaled frequencies $`\kappa =\omega /\mathrm{\Omega }`$ of the $`{}_{2}{}^{}r_{0}^{c}`$-mode as functions of $`\overline{\mathrm{\Omega }}`$ for the polytropic models. The bends of the frequency curves of the $`r`$-modes as $`\mathrm{\Omega }`$ increases are caused by the avoided crossing with the crustal toroidal modes $`{}_{2}{}^{}t_{0}^{}`$, as found in Fig.1 for a realistic neutron star model. As Figure 3 indicates the avoided crossing between the $`{}_{2}{}^{}r_{0}^{c}`$-mode and the $`{}_{2}{}^{}t_{0}^{}`$-mode is a common phenomenon and occurs at around $`\overline{\mathrm{\Omega }}_{cross}`$, which is determined by $`\omega _t(0)`$ for given $`l^{}`$ and $`m`$ (see equation (3)). The apparent difference in behavior of the $`r`$-mode frequencies as a function of $`\mathrm{\Omega }`$ between the three polytropic models is therefore attributable to the difference in the location $`\mathrm{\Omega }_{cross}`$ of the avoided crossing, that is, to the difference in the property $`\mu _0/\rho ^{4/3}`$ of the solid crust.

Here, it is worth while to note that, as shown by Yoshida & Lee (2000b), the deviation of the model from isentropic structure is not important for the fundamental $`r`$-modes with $`l^{}=m`$ and inertial modes, if the Schwarzschild discriminant $`A`$ of the nonisentropic model has sufficiently small absolute values as is expected in neutron stars. As a matter of fact, we can obtain similar results to those shown in Figure 1, even if we artificially set $`A=0`$ in the model NS13T8.

## 3 Dissipation Timescales of Core $`r`$-Modes

For comparison with other studies on the $`r`$-mode instability, let us derive an extrapolation formula for the damping timescale of the fundamental core $`r`$-mode with $`l^{}=m=2`$ by calculating the $`r`$-modes at sufficiently small values of $`\overline{\mathrm{\Omega }}`$, for which the mode coupling with the crustal toroidal mode is negligibly weak. Then, using the same shear and bulk viscosity coefficients as those employed in Yoshida & Lee (2000a) and taking account of the viscous boundary layer (VBL) damping effects following the formulation given in Bildsten & Ushomirsky (2000) (see also Andersson et al (2000) for a correction), we may approximately express the formula as (see e.g., Lindblom et al. 1998)

$`{\displaystyle \frac{1}{\tau }}`$ $`=`$ $`{\displaystyle \frac{1}{\stackrel{~}{\tau }_S}}\left({\displaystyle \frac{10^8\mathrm{K}}{T_c}}\right)^2+{\displaystyle \frac{1}{\stackrel{~}{\tau }_B}}\left({\displaystyle \frac{T_c}{10^8\mathrm{K}}}\right)^6\left({\displaystyle \frac{\mathrm{\Omega }^2}{\pi G\overline{\rho }}}\right)`$ (7)
$`+`$ $`{\displaystyle \frac{1}{\stackrel{~}{\tau }_{VBL}}}\left({\displaystyle \frac{10^8\mathrm{K}}{T_c}}\right)\left({\displaystyle \frac{\mathrm{\Omega }^2}{\pi G\overline{\rho }}}\right)^{1/4}+{\displaystyle \frac{1}{\stackrel{~}{\tau }_{J,2}}}\left({\displaystyle \frac{\mathrm{\Omega }^2}{\pi G\overline{\rho }}}\right)^3,`$
where $`\overline{\rho }`$ is the average density of the star, and $`T_c`$ is the central temperature. The same viscosity coefficients are used both in the fluid and solid regions, and no effects of superfluidity in the core and in the inner crust are considered. Here, the first, second, third and fourth terms in the right-hand side of equation (7) are contributions from the shear viscosity, the bulk viscosity, the shear viscosity in VBL and the current quadrupole radiation, respectively. The factors containing the ratio $`T_c/10^8\mathrm{K}`$ have been introduced in equation (7) to extrapolate the damping rates with respect to the central temperature $`T_c`$.

In Table 2, the various dissipation timescales $`\stackrel{~}{\tau }`$ of the $`{}_{2}{}^{}r_{0}^{c}`$ modes with $`l^{}=m=2`$ are tabulated for the model NS13T8 and a simple $`n=1`$ polytropic model with a static solid crust. For the polytropic model, we employ $`M=1.4M_{\mathrm{}}`$ and $`R=12.5\mathrm{km}`$ and assume that a solid crust forms in the density region of $`\rho 1.5\times 10^{14}\mathrm{g}/\mathrm{cm}^3`$. Note that no wave propagation is assumed in the static crust of the polytropic model employed in this section. Bildsten & Ushomirsky (2000) (see also Anderson et al. 2000, and Rieutord 2000 for more detailed treatment of VBL) used the similar polytropic model with a static solid crust to estimate the VBL damping effects on the $`r`$-modes. To compute the $`r`$-modes of the simple polytropic model, the method given in Yoshida & Lee (2000a), in which effects of the centrifugal force and perturbations of the gravitational potential are taken into account, has been employed. As shown by Table 2, all the dissipative timescales, except that due to the bulk viscosity, have similar values both for the simple polytropic model and for the model NS13T8. We note that although the destabilization due to the current quadrupole radiation for the model NS13T8 is by one order of magnitude stronger than that for the polytropic model, the VBL damping timescales, which are most important among the damping mechanisms of the $`r`$-modes, are almost the same for the two models. We consider that the difference in $`\tau _{J,2}`$ between the two models is caused by the difference in the compactness $`GM/(c^2R)`$ of the models, where $`c`$ denotes velocity of light. The huge difference in the damping timescales associated with the bulk viscosity between the two models may come from the difference in the treatment of the modes in the crust. Although we count the damping contributions in the crust for the model NS13T8, we do not count them for the polytropic model because we have assumed that the crust is static and allows no oscillations. In addition to this, the model NS13T8 in an evolutionary sequence of cooling neutron stars has the crustal temperature that is higher by a factor 7 than the core temperature. Since the bulk viscosity coefficient is proportional to $`T^6`$, the contribution of the bulk viscosity in the crust is largely enhanced for the damping of the $`r`$-modes.

## 4 Discussion and Conclusion

It is well known that retrograde oscillations whose frequency $`\omega `$ satisfies the condition $`0<\omega /\mathrm{\Omega }<m`$ are unstable to the gravitational radiation reaction (Friedman & Schutz 1978). Using equation (2), we can show that the retrograde crustal toroidal modes satisfy this condition and become unstable to the gravitational radiation reaction when

$$\mathrm{\Omega }>\frac{l^{}(l^{}+1)}{m(l^2+l^{}1)}\omega _t(0),$$
(8)
where $`\omega _t(0)`$ denotes the oscillation frequency of the mode at $`\mathrm{\Omega }=0`$. This instability can happen at slow rotation rates of the star if $`\omega _t(0)`$ is sufficiently small. Note that all crustal toroidal modes shown in Figure 1 are unstable to the gravitational radiation reaction. Our preliminary stability calculation, however, shows that the instability of the crustal toroidal mode is very weak and may not be important.

In this paper, we have investigated the properties of $`r`$-modes of slowly rotating neutron stars with a solid crust. We have found that the mode property of the core $`r_0`$ mode with $`l^{}=m=2`$, which is most important for the $`r`$-mode instability of neutron stars, is strongly affected by mode coupling with the crustal $`{}_{2}{}^{}t_{0}^{}`$ modes as $`\mathrm{\Omega }`$ increases. This means that we cannot assume a simple mode property of the core $`r_0`$ modes with $`l^{}=m`$ as a function of $`\overline{\mathrm{\Omega }}`$. To discuss the $`r`$-mode instability of rotating neutron stars it has been common to derive an extrapolation formula for the damping timescale of the $`r`$-modes (e.g., Lindblom et al. 1998). The extrapolation is usually carried out about the interior temperature $`T`$ and the rotation frequency $`\mathrm{\Omega }`$. The extrapolation with respect to the temperature may be justified if the stars have an isothermal structure characterized, for example, by the central temperature $`T_c`$. On the other hand, since the mode properties of the $`r`$-mode are well known only in the limit of $`\overline{\mathrm{\Omega }}0`$, the damping timescales are calculated in this limit and are extrapolated to rapidly rotating neutron stars with $`\overline{\mathrm{\Omega }}1`$. This process of the extrapolation with respect to $`\mathrm{\Omega }`$ is justified if the mode frequency $`\omega `$ is well approximated by $`2m\mathrm{\Omega }/l^{}(l^{}+1)`$ and the eigenfunction has the dominant toroidal component, independent of $`\mathrm{\Omega }`$. We know that the above two assumptions are satisfied for the fundamental $`r`$-modes with $`l^{}=m`$ for fluid stars (Yoshida & Lee 2000b; Yoshida et al. 2000; Karino et al. 2000). In the case of neutron stars with a solid crust, however, the two assumptions are not satisfied except when $`\overline{\mathrm{\Omega }}0`$ and no mode couplings are expected. For example, when the $`r`$-mode is near the avoided crossing with the crustal $`{}_{2}{}^{}t_{0}^{}`$ mode, the eigenfunctions and eigenfrequencies of the $`r`$-mode behave peculiarly and change rapidly as functions of $`\mathrm{\Omega }`$. This makes it difficult to predict reliably the evolution of a neutron star driven by the $`r`$-mode instability through the avoided crossing by simply extrapolating the damping timescales obtained at small $`\mathrm{\Omega }`$. In addition to this, the eigenfunctions of the $`r`$-modes are contaminated by the crustal torsional modes at large $`\mathrm{\Omega }`$, and the extrapolated damping timescales are not necessarily accurate enough to give the useful estimations. These difficulties mean that the extrapolation formula of the damping timescales of the $`r`$-modes, derived in the limit of $`\overline{\mathrm{\Omega }}0`$, is not applicable to the discussion on the $`r`$-mode instability of the rapidly rotating neutron stars with a solid crust unless the effects of the avoided crossings with and the contamination of the eigenfunctions by the crustal torsional modes are correctly taken into account. In this sense, discussions based on the extrapolation formula of the damping timescales are not fully justified for neutron stars with a solid crust.


format: text

text:
# Extremal Paths on a Random Cayley Tree

## I Introduction

Optimization lies at the heart of a vast number of phenomena: Particles “seek” paths with minimal action, species “try” to maximize fitness, companies minimize cost. Many optimization problems, e.g., the traveling salesman problem, are notoriously hard. Additionally, optimization problems often involve randomness which makes them even more complicated. When the number of random entries is very large, however, optimization problems might become simpler thanks to the “concentration of the measure”. This means that the probability distribution of some random variable becomes highly localized, almost like a delta function concentrated around its average value. This phenomenon is well-known in probability theory where, for instance, it accounts for the fact that after flipping a coin $`N`$ times, the probability that the number of heads exceeds $`N/2`$ by more than, say, $`100\sqrt{N}`$ is about $`10^{800}`$. Below, we shall investigate an optimization problem in random media which demonstrates this concentration of the measure in the strongest sense (the variance of the optimized quantity is finite) and displays an unbinding phase transition. A remarkable hidden connection with traveling wave phenomena allows to explain both the concentration of the measure and the phase transition.

The optimization problem considered in this work can be formulated as follows: Take a rooted tree whose bonds have random lengths and find descending paths of extremal total length. We assume that a random length (or energy) is assigned to each bond of the tree, and that lengths are independent and chosen from the same probability distribution $`\rho (l)`$. The total length of a path is the sum of lengths of the bonds along the path, and we want to determine the minimal (maximal) length among the paths from the root to the bottom of the tree. We focus on a rooted Cayley tree where each node, except the root, has coordination number 3; the coordination number of the root is 2. However, our main results can be easily generalized to the case of a tree with arbitrary coordination number, and even to trees where the coordination number at any node is random.

Similar problems arise in numerous fields ranging from computer science to condensed matter physics, where it is known as the problem of directed polymers in random media. Indeed, the minimal path problem can be considered as the zero temperature limit of the problem of directed polymer on Cayley trees. Directed polymers on Cayley trees have been investigated in the past, and recently this problem has resurfaced in a surprisingly large number of apparently unrelated problems. However, the emphasis of Refs. was on the spin glass like transition occurring at a finite temperature. Here, we consider exactly zero temperature and focus on the unbinding transition driven by the parameter $`p`$ of the bimodal distribution, $`\rho (l)=p\delta _{l,1}+(1p)\delta _{l,0}`$. Recently, similar unbinding (or “depinning”) phase transitions have been found in a number of non-equilibrium processes without quenched disorder.

Let us first outline our main results. We derive exactly the statistics of both the minimal and maximal path on a random Cayley tree for arbitrary distribution $`\rho (l)`$ of the edge lengths. We find a class of distributions for which the minimal path undergoes an unbinding phase transition from a localized phase to a moving phase as a parameter of the distribution is varied. A particular example, which we study in detail in this paper, is the bimodal distribution. In the bimodal case, length of each bond can be either $`1`$ with probability $`p`$ or $`0`$ with probability $`(1p)`$. For any path starting at the root and moving down, the length of the path is the sum of the lengths of the bonds along the path, i.e., the number of $`1`$’s along the path. There are $`2^n`$ paths from the root to the $`n^{\mathrm{th}}`$ level. We show that the minimal length exhibits an unbinding transition as the parameter $`p`$ varies between 0 and 1. There exists a critical value $`p_c=1/2`$ such that the average minimal length grows linearly with $`n`$ when $`p>1/2`$, but gets pinned, i.e., remains finite for large $`n`$ when $`p<1/2`$. More specifically, we find that for large $`n`$, the average minimal length $`L_n^{\mathrm{min}}`$ behaves as

$$L_n^{\mathrm{min}}\{\begin{array}{cc}v_{\mathrm{min}}(p)n\hfill & p>1/2\text{,}\hfill \\ (\mathrm{ln}2)^1\mathrm{ln}\mathrm{ln}n\hfill & p=1/2\text{,}\hfill \\ \mathrm{finite}\hfill & p<1/2\text{.}\hfill \end{array}$$
(1)
For $`p>1/2`$, we shall later compute the minimal “velocity” $`v_{\mathrm{min}}(p)`$ exactly. The average maximal length, on the other hand, behaves as

$$L_n^{\mathrm{max}}v_{\mathrm{max}}(p)n\text{ for all p},$$
(2)
where $`v_{\mathrm{max}}(p)+v_{\mathrm{min}}(1p)=1`$. This duality relation is very general and valid for directed paths on arbitrary graphs, i.e., for finite-dimensional lattices. We also derive a generalized duality relation between minimal and maximal velocities for a certain class of bounded distributions.

The minimal length problem on a rooted Cayley tree, where the coordination number at each node is random, also displays an unbinding transition at $`p_c=1b^1`$, where $`b`$ is the average number of branches per node. At the critical point, the average minimal length $`L_n^{\mathrm{min}}`$ exhibits the same double logarithmic growth for any $`b>1`$, and above $`p_c`$ the minimal length grows linearly with $`n`$ with velocity $`v_{\mathrm{min}}(p,b)`$.

The rest of this paper is organized as follows. In Section II, we study the unbinding transition of the minimal path for the special bimodal distribution in detail. In Section III, we study the maximal path for the bimodal case and derive the duality relation. In Section IV, we generalize the results to a Cayley tree with random branching. In Section IV, we go beyond the bimodal distribution and obtain generalized results for arbitrary distribution of the edge length. Finally we conclude in Section VI with a brief summary and outlook.

## II Minimal Path(s)

Instead of studying $`P_n^{\mathrm{min}}(x)=\mathrm{Prob}(L_n^{\mathrm{min}}=x)`$, it proves convenient to consider the cumulative distribution, $`P_n(x)=\mathrm{Prob}(L_n^{\mathrm{min}}x)`$. Clearly, $`P_n(x)`$ is the probability that all possible paths to the $`n^{\mathrm{th}}`$ level have lengths $`x`$. Once we get $`P_n(x)`$, the distribution of the minimal length is found from the obvious identity $`P_n^{\mathrm{min}}(x)=P_n(x)P_n(x+1)`$.

It is easy to see that $`P_n(x)`$ satisfies the following recursion relation

$$P_{n+1}(x)=\left[(1p)P_n(x)+pP_n(x1)\right]^2.$$
(3)
This relation can be derived by analyzing various possibilities for the lengths of the two edges issuing from the root and taking into account that two subsequent daughter trees are statistically independent. For instance, both lengths are equal to zero with probability $`(1p)^2`$, and paths in the subsequent trees have lengths $`x`$ with probability $`P_n^2(x)`$; this provides the contribution $`(1p)^2P_n^2(x)`$ to $`P_{n+1}(x)`$. Similarly, one gets the contributions $`2p(1p)P_n(x1)P_n(x)`$ and $`p^2P_n^2(x1)`$, which sum up to Eq. (3).

We must solve Eq. (3) subject to the initial condition

$$P_0(x)=\{\begin{array}{cc}1\hfill & x0\text{,}\hfill \\ 0\hfill & x>0\text{.}\hfill \end{array}$$
(4)
Clearly, the minimal length $`L_n^{\mathrm{min}}`$ is a random variable which takes values between $`0`$ and $`n`$, and Eqs. (3)–(4) indeed show that $`P_n(x)=1`$ for $`x0`$ and $`P_n(x)=0`$ for $`x>n`$.

When $`n`$ grows, $`L_n^{\mathrm{min}}`$ increases as well and thus it should reach a limit. For $`p<1/2`$, this limit is finite, i.e., one can find an infinite path with only finite number of edges of length 1. Mathematically, it means that $`P_n(x)`$ approaches a stationary distribution, $`P_n(x)P(x)`$, which satisfies

$$P(x)=\left[(1p)P(x)+pP(x1)\right]^2.$$
(5)
Starting from $`P(0)=1`$, one computes $`P(x)`$ recursively: $`P(1)=p^2(1p)^2`$, etc. $`P(x)`$ vanishes extremely fast in the large $`x`$ limit: $`\mathrm{ln}P(x)2^x\mathrm{ln}p`$.

For $`p1/2`$, the minimal length grows as $`n`$ increases. To understand why this is so, recall that $`P_n(0)1`$ and look at $`P_n(1)`$. By inserting $`P_n(1)=1q_n`$ into Eq. (3) we find $`q_{n+1}=2(1p)q_n(1p)^2q_n^2`$. Hence $`P_n(1)1`$ for $`p1/2`$. Proceeding this line of reasoning one can verify that $`P_n(x)1`$ when $`n\mathrm{}`$. Additionally, we see the difference between the cases of $`p>1/2`$ and $`p=1/2`$: In the former situation, $`P_n(x)1`$ exponentially fast, while in the latter situation the approach is algebraic. Therefore, the size of the region where $`P_n(x)`$ is close to 1 should grow linearly with $`n`$ when $`p>1/2`$. Thus the distribution $`P_n(x)`$ approaches the traveling wave form (see Fig. (1)),

$$P_n(x)𝒫^2(y),y=xL_n^{\mathrm{min}}.$$
(6)
In Eq. (6), we have chosen $`y=xL_n^{\mathrm{min}}`$ as the wave variable. One can justify this choice by noting that the relation $`P_n^{\mathrm{min}}(x)=P_n(x)P_n(x+1)`$ leads to identity

$$L_n^{\mathrm{min}}=\underset{x=0}{\overset{\mathrm{}}{}}x\left[P_n(x)P_n(x+1)\right]=\underset{x=1}{\overset{\mathrm{}}{}}P_n(x),$$
(7)
and hence $`L_n^{\mathrm{min}}`$ is indeed an appropriate characterization of the location of the front. By inserting the wave form (6) into Eq. (3) and using $`L_n^{\mathrm{min}}vn`$ we find,

$$𝒫(yv)=(1p)𝒫^2(y)+p𝒫^2(y1),$$
(8)
The solution of this equation looks like a $`[10]`$ waveform with $`P(y)1`$ as $`y\mathrm{}`$ and $`P(y)0`$ as $`y\mathrm{}`$. While it is very hard to solve this nonlinear, nonlocal equation exactly, the velocity can be found by analyzing the tail region, $`y\mathrm{}`$, where $`1𝒫(y)`$ is small. Linearizing Eq. (8) in this region and noting that it admits an exponential solution, $`1𝒫(y)e^{\lambda y}`$, one finds that the velocity $`v_\lambda `$ is related to the decay exponent $`\lambda `$ via

$$v_\lambda =\frac{\mathrm{ln}\left[2(1p)+2pe^\lambda \right]}{\lambda }.$$
(9)
When $`\lambda >\mathrm{ln}\frac{2p}{2p1}`$, the velocity $`v_\lambda `$ is positive. While any such $`\lambda `$ is in principle allowed, and thus the entire velocity range of is feasible, a particular value is actually selected. This is similar to the velocity selection in a large class of problems with a traveling wave front. It is well known that for a wide class of initial conditions, the extremum value is chosen. From this general front selection principle one can infer that in the present case the selected value $`\lambda =\lambda ^{}(p)`$ corresponds to the maximum of $`v_\lambda `$, and hence the selected “velocity” is $`v_{\mathrm{min}}(p)=v_\lambda ^{}`$. Thus, $`\lambda ^{}`$ is a solution of the equation

$$\mathrm{ln}\left[2(1p)+2pe^\lambda ^{}\right]+\frac{p\lambda ^{}e^\lambda ^{}}{1p+pe^\lambda ^{}}=0,$$
(10)
and the selected velocity is

$$v_{\mathrm{min}}(p)=\frac{pe^\lambda ^{}}{1p+pe^\lambda ^{}}.$$
(11)
Although it seems impossible to provide an explicit expression for $`v_{\mathrm{min}}(p)`$, one can easily probe its behavior in the limiting cases $`p\frac{1}{2}`$ and $`p1`$. In these respective limits, one gets

$$v_{\mathrm{min}}(p)\frac{2p1}{\mathrm{ln}\frac{1}{2p1}}$$
(12)
and

$$1v_{\mathrm{min}}(p)\frac{\mathrm{ln}2}{\mathrm{ln}\frac{1}{1p}}.$$
(13)
One can also derive the leading correction to the dominant linear behavior of $`L_n^{\mathrm{min}}`$,

$$L_n^{\mathrm{min}}v_{\mathrm{min}}(p)n+\frac{3}{2\lambda ^{}}\mathrm{ln}n.$$
(14)
The logarithmic correction becomes especially important when $`p\frac{1}{2}`$, as can be seen from the asymptotics

$$L_n^{\mathrm{min}}\frac{ϵn}{\mathrm{ln}(1/ϵ)}+\frac{3\mathrm{ln}n}{2\mathrm{ln}(1/ϵ)},ϵ=2p1.$$
(15)
The logarithmic correction to the front position was first derived by Bramson in the context of a reaction-diffusion equation , and was subsequently re-derived and generalized by a number of authors. Our derivation of Eq. (14) follows an approach of Ref., and we only outline the main steps. First, we need to examine the exact equation (3) rather than the $`n\mathrm{}`$ limit, Eq. (8), which was sufficient to determine the velocity. However, the analysis is still performed in the region far behind the front, so we can use the linearized version of Eq. (3). Writing $`P_n(x)=1Q_n(x)`$, plugging it into Eq. (3), and ignoring quadratic terms gives

$$Q_{n+1}(x)=2(1p)Q_n(x)+2pQ_n(x1).$$
(16)
Now we seek a solution which has the form

$$Q_n(x)=n^\alpha G\left(yn^\alpha \right)e^{\lambda ^{}y},$$
(17)
with $`y=xL_n^{\mathrm{min}}`$ as previously and yet undetermined exponent $`\alpha `$ and the scaling function $`G(z)`$. The correction to the leading term in the front position should be incorporated as,

$$L_n^{\mathrm{min}}v_{\mathrm{min}}(p)n+c(n)$$
(18)
where the function $`c(n)`$ is yet undetermined. Plugging Eqs. (17)–(18) into Eq. (16) we find that different leading orders are compatible provided that $`\alpha =1/2`$ and $`c(n)=\beta \mathrm{ln}n`$ with some constant $`\beta `$. Additionally, the scaling function $`G(z)`$ satisfies a parabolic cylinder equation

$$\frac{d^2G}{dz^2}+z\frac{dG}{dz}+(2\beta \lambda ^{}1)G(z)=0.$$
(19)
Equation (19) should be solved subject to the boundary conditions $`G(z)0`$ as $`z\mathrm{}`$ as $`Q_n(x)`$ must vanish for $`x\mathrm{}`$ and $`G(z)z`$ for $`z0`$ to ensure that $`Q_n(x)`$ is independent on $`n`$ in the limit $`n\mathrm{}`$. The boundary condition in the large $`z`$ limit selects one of the two possible solutions: $`G(z)=A\mathrm{exp}(z^2/4)D_{2(\beta \lambda ^{}1)}(z)`$, where $`D_\nu `$ is the parabolic cylinder function with index $`\nu `$. The second boundary condition $`G(z)z`$ fixes the index of the parabolic cylinder function, $`2(\beta \lambda ^{}1)=1`$. Hence, $`\beta =3/2\lambda ^{}`$ thus completing the derivation of Eq. (14). Note that the sign of the logarithmic correction term in Eq. (14) is positive. This is different from the most earlier studied problems although the positive sign has been occasionally seen, see Ref..

In the the critical case of $`p=1/2`$, the convergence of the distribution $`P_n(x)`$ towards the asymptotic value is only algebraic, and therefore the front propagates extremely slowly. The simplest way to determine the rate of propagation is again to look at $`P_n(x)`$ far behind the front. Writing $`P_n(x)=1Q_n(x)`$ and plugging it into Eq. (3) yields

$`Q_{n+1}(x)`$ $`=`$ $`Q_n(x)+Q_n(x1)`$ (20)
$``$ $`\left[{\displaystyle \frac{Q_n(x)+Q_n(x1)}{2}}\right]^2.`$ (21)
$`Q_n(x)=0`$ for $`x0`$, so the first non-trivial $`Q_n`$’s are $`Q_n(1)q(n)`$. From Eq. (20), $`q(n+1)=q(n)\frac{q^2(n)}{4}`$. In the large $`n`$ limit, we can employ the continuum approximation to get $`\frac{dq}{dn}=\frac{q^2}{4}`$ and thence $`Q_n(1)4n^1`$. Similarly, we find the asymptotics $`Q_n(2)4n^{1/2}`$. Generally, Eq. (20) implies $`Q_n(x)2\sqrt{Q_n(x1)}`$, from which

$$Q_n(x)4n^{2^{(x1)}}=4\mathrm{exp}\left(2e^{\mathrm{ln}\mathrm{ln}nx\mathrm{ln}2}\right).$$
(22)
This demonstrates that at the critical point, $`Q_n(x)`$ also has a traveling wave form, $`Q_n(x)=g(xx_n)`$, with the front at $`x_n(\mathrm{ln}2)^1\mathrm{ln}\mathrm{ln}n`$.

Equation (22) formally applies for $`xx_n`$. To investigate $`Q_n(x)`$ in the entire range of $`x`$, we again use the fact that the distribution $`Q_n(x)`$ should approach a traveling wave form, $`Q_n(x)Q(y)`$ with $`y=n^{2^{(x1)}}`$, as $`n\mathrm{}`$. By inserting this into (20) and taking $`n\mathrm{}`$ limit, we arrive at

$$Q(y)=2\sqrt{Q(y^2)}Q(y^2).$$
(23)
From this equation, one finds

$$Q(y)=4y4y^22y^3+\mathrm{}\mathrm{for}y0.$$
(24)
The first term in this expansion indeed agrees with Eq. (22). Similarly, one gets

$$1Q(y)\mathrm{exp}\left[\frac{\mathrm{const}}{1y}\right]\mathrm{for}y1.$$
(25)
To provide a more rigorous derivation of the growth law for $`L_n^{\mathrm{min}}`$, note that $`L_n^{\mathrm{min}}=_{x1}P_n(x)=n_{x1}Q_n(x)`$. Inserting $`Q_n(x)=Q(y)`$ and replacing the sum by an integral (valid for large $`n`$), we get,

$$L_n^{\mathrm{min}}n\frac{1}{\mathrm{ln}2}_{2^{n1}\mathrm{ln}n}^{\mathrm{ln}n}Q(e^z)\frac{dz}{z}.$$
(26)
In the limit of large $`n`$, the most important contribution in the integral comes from its lower limit. Using $`Q(y)1`$ as $`y1`$ (see Eq. (25)), we arrive at the desired result,

$$L_n^{\mathrm{min}}(\mathrm{ln}2)^1\mathrm{ln}\mathrm{ln}n.$$
(27)
## III Maximal Path(s)

We now turn to the statistics of $`L_n^{\mathrm{max}}`$, the length of the maximal path for the bimodal distribution. It is now convenient to define $`R_n(x)=\mathrm{Prob}(L_n^{\mathrm{max}}x)`$. By proceeding as in the minimal length problem, it is easy to show that $`R_n(x)`$ evolves according to the same equation (3) with $`P`$ replaced by $`R`$, i.e.,

$$R_{n+1}(x)=\left[(1p)R_n(x)+pR_n(x1)\right]^2.$$
(28)
The only difference from the minimal case is in the initial condition. Instead of Eq. (4), we now have

$$R_0(x)=\{\begin{array}{cc}1\hfill & x0\text{,}\hfill \\ 0\hfill & x<0\text{.}\hfill \end{array}$$
(29)
This difference in the initial condition, however, leads to a different behavior of the maximal front as demonstrated below.

Not surprisingly, equation (28) admits a traveling wave solution, $`R_n(x)=R(y=xv_{\mathrm{max}}n)`$ for all $`p`$ (see Fig. (2)). However, in contrast to the $`[10]`$ waveform of the minimal case, the function $`R(y)`$ now looks like a $`[01]`$ waveform with $`R(y)0`$ as a double exponential behind the front ($`y\mathrm{}`$) and $`R(y)1`$ exponentially ahead the front ($`y\mathrm{}`$). We then analyze the wavefront near the ‘forward’ tail, as opposed to the ‘backward’ tail of the minimal problem. Substituting $`1R(y)e^{\mu y}`$, we find that there exists a family of traveling wave solutions with velocity $`v_\mu `$ parametrized by $`\mu `$,

$$v_\mu =\frac{\mathrm{ln}\left[2(1p)+2pe^\mu \right]}{\mu }.$$
(30)
This velocity $`v_\mu `$ has a unique minimum at $`\mu =\mu ^{}`$. By the front selection mechanism, this minimum value $`v_\mu ^{}`$ is selected as the velocity $`v_{\mathrm{max}}(p)`$ of the front. Thus, we have

$$\mathrm{ln}\left[2(1p)+2pe^\mu ^{}\right]\frac{p\mu ^{}e^\mu ^{}}{1p+pe^\mu ^{}}=0,$$
(31)
and the selected maximal velocity is

$$v_{\mathrm{max}}(p)=\frac{pe^\mu ^{}}{1p+pe^\mu ^{}}.$$
(32)
Using these results together with Eqs. (10) and (11), one immediately finds

$$v_{\mathrm{min}}(p)+v_{\mathrm{max}}(1p)=1.$$
(33)
This is the duality relation. Unlike the behavior of the minimal length, the maximal length does not undergo an unbinding transition in the sense that it increases linearly with $`n`$ for large $`n`$ for all $`0p1`$. However, there is still a transition at $`p=1/2`$: The velocity $`v_{\mathrm{max}}(p)`$ increases from $`0`$ to $`1`$ as $`p`$ increases from $`0`$ to $`1/2`$ and then sticks to $`1`$ for $`p1/2`$ (see Fig. (3)).

For the maximal front also, one can easily compute the sub-leading logarithmic correction to the front position. Proceeding as in the minimal case, we find

$$L_n^{\mathrm{max}}v_{\mathrm{max}}(p)n\frac{3}{2\mu ^{}}\mathrm{ln}n.$$
(34)
Note the negative sign in the correction term, as opposed to the positive sign for the minimal case in Eq. (14). This is also consistent with the duality relation in Eq. (33). Thus, while the coefficient $`3/2`$ of the logarithmic correction term seems to be universal not just in the present problem but in many other cases, the sign of the correction term is not and can either be positive or negative.

The duality relation (33) can also be derived by a general argument which does not depend on the structure of the underlying tree and therefore is valid for extremal directed paths on arbitrary graphs, in particular, on arbitrary lattices. The argument follows from the observation that for bimodal distribution, if one replaces the $`0`$’s by $`1`$’s and the $`1`$’s by $`0`$’s on the minimal path (and thereby replacing $`p`$ by $`1p`$), then the minimal path becomes the maximal path. Let $`n_0`$ and $`n_1`$ denote respectively the number of $`0`$’s and $`1`$’s on the minimal path. Evidently, $`n_0+n_1=n`$. Then clearly, $`L_n^{\mathrm{min}}(p)=n_1`$ and by duality, $`L_n^{\mathrm{max}}(1p)=n_0`$. Adding the two quantities, we get $`L_n^{\mathrm{min}}(p)+L_n^{\mathrm{max}}(1p)=n`$. Dividing by $`n`$, we immediately get the duality relation Eq. (33). Below, we shall derive a more generalized version of this duality relation which goes beyond the bimodal distribution.

## IV Generalization to Trees with Random Branching

The generalization of the above analysis to the case of the rooted Cayley tree with a random number of branches is straightforward. Let $`b_m`$ is the probability that the number of branches is equal to $`m`$ and $`b`$ is the average number of branches, $`b=_{m1}mb_m`$. Equation (3) is replaced by

$$P_{n+1}(x)=\underset{m=1}{\overset{\mathrm{}}{}}b_m\left[(1p)P_n(x)+pP_n(x1)\right]^m,$$
(35)
and the subsequent analysis repeats the steps detailed in the case of the binary tree. In particular, $`P_n(x)`$ approaches the stationary distribution for $`p<p_c`$, while above $`p_c`$ the minimal length grows according to Eq. (14) with $`\lambda ^{}`$ a solution to equation

$$\mathrm{ln}\left[b(1p)+bpe^\lambda ^{}\right]+\frac{p\lambda ^{}e^\lambda ^{}}{1p+pe^\lambda ^{}}=1,$$
(36)
and $`v_{\mathrm{min}}(p)`$ given by the same relation (11). At the critical point $`p_c=1b^1`$, the minimal length is given by the same expression (27) for any tree. Also, the velocity of the maximal path can be determined from the general duality relation in Eq. (33).

## V Generalization to Arbitrary Distributions

We now consider the statistics of extreme path(s) on a binary tree for arbitrary distribution $`\rho (l)`$ of the edge lengths. Let us first consider the minimal path(s). The velocity of the minimal path for an arbitrary distribution can be extracted from earlier results by Derrida and Spohn which we re-derive below for the sake of completeness. We consider again $`P_n(x)=\mathrm{Prob}(L_n^{\mathrm{min}}x)`$. Proceeding as in the case of the bimodal distribution, one finds that the appropriate generalization of Eq. (3) is given by

$$P_{n+1}(x)=\left[_0^{\mathrm{}}𝑑l\rho (l)P_n(xl)\right]^2.$$
(37)
with the same initial condition, Eq. (4), as earlier. The distribution of the minimal length is then given by $`P_n^{\mathrm{min}}(x)=dP_n(x)/dx`$ and the average minimal length is $`L_n^{\mathrm{min}}=_0^{\mathrm{}}P_n(x)𝑑x`$.

Defining $`P_n(x)=G_n^2(x)`$, we recast Eq. (37) into

$$G_{n+1}(x)=_0^{\mathrm{}}𝑑l\rho (l)G_n^2(xl).$$
(38)
For any given $`n`$, $`G_n(x)1`$ as $`x\mathrm{}`$ and $`G_n(x)0`$ as $`x\mathrm{}`$. Substituting $`G_n(x)=1g_n(x)`$ in Eq. (38) and using normalization of $`\rho (l)`$ one finds,

$$g_{n+1}(x)=_{\mathrm{}}^x𝑑y\rho (xy)\left[2g_n(y)g_n^2(y)\right],$$
(39)
where $`g_n(x)0`$ as $`x\mathrm{}`$ and $`g_n(x)1`$ as $`x+\mathrm{}`$. We analyze the above equation in the ‘backward’ tail region, i.e., when $`x\mathrm{}`$. In this limit, one can neglect the nonlinear term in $`g_n`$ inside the integral in Eq. (39). The resulting linear equation admits a traveling wave solution of the form, $`g_n(x)=e^{\lambda (xv_\lambda n)}`$, with

$$v_\lambda =\frac{1}{\lambda }\mathrm{ln}\left[2_0^{\mathrm{}}𝑑ze^{\lambda z}\rho (z)\right].$$
(40)
For generic length distributions $`\rho (l)`$, this function has a maximum at $`\lambda =\lambda ^{}`$ and by the general velocity selection principle, this maximum velocity is selected, i.e., $`v_{\mathrm{min}}=v_\lambda ^{}`$.

The question we are interested in is for which class of distributions $`\rho (l)`$, an unbinding transition can occur for the minimal path. Such a transition will occur if the velocity $`v_\lambda `$ in Eq. (40) ceases to have a unique maximum. By analyzing Eq. (40) one sees that this can happen only if $`\rho (l)`$ has a nonzero delta-function weight at $`l=0`$, i.e., when $`\rho (l)=a\delta (l)+(1a)f(l)`$ with $`0<a<1`$ and $`f(l)`$ is some positive function with $`_0^{\mathrm{}}𝑑lf(l)=1`$. The unbinding transition occurs as the parameter $`a`$ is varied. Note that the positivity condition of velocity in Eq. (40) demands that $`2a<1`$. Thus the critical point always occurs at $`a=1/2`$. The bimodal distribution considered in the previous sections is a special case of this class of distributions with $`a=1p`$ and $`f(l)=\delta (l1)`$.

For $`a<1/2`$, the average minimal length increases linearly with $`n`$ for large $`n`$ with the velocity $`v_\lambda ^{}`$ obtained from Eq. (40). For $`a>1/2`$, the function $`P_n(x)`$ reaches a stationary distribution for large $`n`$ and hence the $`L_n^{\mathrm{min}}`$ saturates to a nonzero constant as $`n\mathrm{}`$. What happens at the critical point $`a=1/2`$ for generic distributions $`f(l)`$’s? For $`a=1/2`$, the Eq. (38) reduces to,

$$G_{n+1}(x)=\frac{1}{2}G_n^2(x)+\frac{1}{2}_0^{\mathrm{}}𝑑lf(l)G_n^2(xl).$$
(41)
Detailed analysis of the above equation reveals a rather remarkable universal property at the critical point. It turns out that for generic $`f(l)`$, there are only two possible behaviors of $`L_n^{\mathrm{min}}=_0^{\mathrm{}}𝑑xG_n^2(x)`$, depending on whether the function $`f(l)`$ has a gap or not near $`l=0`$. If $`f(l)`$ does not have a gap at $`l=0`$, it turns out that as $`n\mathrm{}`$, $`G_n(x)`$ in Eq. (41) tends to a stationary distribution and hence $`L_n^{\mathrm{min}}\mathrm{const}`$. In the opposite case, when $`f(l)`$ has a gap $`\mathrm{\Delta }`$ near $`l=0`$, it turns out that as $`n\mathrm{}`$,

$$L_n^{\mathrm{min}}\frac{\mathrm{\Delta }}{\mathrm{ln}2}\mathrm{ln}\mathrm{ln}n.$$
(42)
The bimodal distribution corresponds to the special case $`\mathrm{\Delta }=1`$.

These general results are best demonstrated by specific examples which can be worked out explicitly. Let us first consider an example which is gapless, such as the exponential distribution, $`f(l)=e^l`$. Remarkably, Eq. (41) can be solved exactly in this case. It is easy to verify that Eq. (41) admits an exponential solution, $`G_n(x)=A_ne^x`$, where $`A_{n+1}=(1+A_n^2)/2`$ with $`A_0=0`$. As $`n\mathrm{}`$, $`A_n`$ approaches to $`1`$ (more precisely, $`A_n12/n`$). Therefore, $`L_n^{\mathrm{min}}=𝑑xG_n^2(x)=A_n^2/2`$ approaches to 1/2 as $`n\mathrm{}`$. Next we consider a distribution with gap such as $`f(l)=\frac{1}{4}\delta (l1)+\frac{1}{4}\delta (l2)`$. This case also can be worked out explicitly following the same steps as used for the bimodal case. It turns out that essentially all the steps are identical to the bimodal case, except that $`n`$ in Eq. (27) gets replaced by $`n/2`$. But this does not change the leading asymptotic behavior for large $`n`$ which is still given by Eq. (27).

We now turn to the maximal path(s) for an arbitrary distribution $`\rho (l)`$. As in the case of the bimodal distribution, $`R_n(x)=\mathrm{Prob}(L_n^{\mathrm{max}}x)`$ satisfies the same equation (37) as $`P_n(x)`$, i.e.,

$$R_{n+1}(x)=\left[_0^{\mathrm{}}𝑑l\rho (l)R_n(xl)\right]^2,$$
(43)
the only difference is that the initial condition, Eq. (4), should be replaced by Eq. (29). Substituting $`R_n(x)=[1s_n(x)]^2`$, we get

$$s_{n+1}(x)=_{\mathrm{}}^x𝑑y\rho (xy)\left[2s_n(y)s_n^2(y)\right],$$
(44)
with the boundary conditions $`s_n(x)1`$ as $`x\mathrm{}`$ and $`s_n(x)0`$ as $`x\mathrm{}`$. We now have to analyze the ‘forward’ tail of the front, i.e., the behavior in the $`x\mathrm{}`$ limit. We first consider bounded distributions $`\rho (l)`$’s with an upper cut-off $`\mathrm{\Lambda }`$. In that case, the lower limit of the integration in Eq. (44) becomes $`x\mathrm{\Lambda }`$. In the tail region where $`x\mathrm{\Lambda }`$, one can again neglect the nonlinear term in $`s_n`$ inside the integral in Eq. (44) and the resulting linear equation admits a traveling wave solution of the form $`s_n(x)=e^{\mu (xv_\mu n)}`$, where

$$v_\mu =\frac{1}{\mu }\mathrm{log}\left[2_0^\mathrm{\Lambda }𝑑ze^{\mu z}\rho (z)\right].$$
(45)
For generic distribution, $`v_\mu `$ has a unique minimum at $`\mu =\mu ^{}`$ and via the front selection principle, this minimum velocity is chosen and one gets, $`v_{\mathrm{max}}=v_\mu ^{}`$. Comparing Eqs. (40) and (45) it becomes evident that for the bounded symmetric distributions, i.e., when $`\rho (l)=\rho (\mathrm{\Lambda }l)`$, one gets the general relation,

$$v_{\mathrm{min}}+v_{\mathrm{max}}=\mathrm{\Lambda }.$$
(46)
For unbounded distributions, however, one has to be careful since it is not obvious that one can neglect the nonlinear term in $`s_n`$ inside the integral in Eq. (44). Equation (45) for the maximal velocity still remains valid as long as $`\rho (l)`$ decays with $`l`$ exponentially or faster. To see this explicitly in an example, let us consider the exponential distribution, $`\rho (l)=e^l`$. In this case, one can transform the the integral equation (44) into the following difference-differential equation,

$$\frac{ds_{n+1}(x)}{dx}=2s_n(x)s_n^2(x)s_{n+1}(x).$$
(47)
In this differential form, it is clear that for large $`x`$, one can neglect the nonlinear term on the right-hand side of Eq. (47). By inserting $`s_n(x)e^{\mu (xv_\mu n)}`$ into resulting linear equation, one gets

$$v_\mu =\frac{1}{\mu }\mathrm{ln}\left(\frac{2}{1\mu }\right).$$
(48)
Note that the same formula is obtained if one substitutes $`\rho (l)=e^l`$ directly into the general expression (45). Minimizing with respect to $`\mu `$, one gets $`v_{\mathrm{max}}`$. The corresponding minimal velocity $`v_{\mathrm{min}}`$ for the exponential distribution can be obtained from the general formula in Eq. (40) with $`\rho (z)=e^z`$ and then maximizing with respect to $`\lambda `$. So, for the exponential distribution we finally get,

$$v_{\mathrm{min}}=0.231961\mathrm{}\mathrm{and}v_{\mathrm{max}}=2.678346\mathrm{}$$
(49)
Thus for unbounded distributions, there does not seem to exist any simple relation between minimal and maximal velocities.

## VI Conclusion

In this work, we have shown that the length of the minimal path on a Cayley tree exhibits an unbinding phase transition from a localized phase to a moving phase. This phase transition is driven by the parameter $`p`$ of the bimodal distribution. In the localized phase ($`p<p_c`$), the minimal length distribution approaches a stationary depth-independent form $`P^{\mathrm{min}}(x)`$. This distribution vanishes extremely sharply, as the double exponential, in the large $`x`$ limit. In the moving phase ($`p>p_c`$), the minimal length distribution $`P_n^{\mathrm{min}}(x)`$ approaches a traveling wave form with the front velocity relaxing algebraically, $`v_{\mathrm{min}}=v(p)+3/(2\lambda ^{}n)+𝒪\left(n^{3/2}\right)`$, in the large $`n`$ limit. Specifically, $`P_n^{\mathrm{min}}(x)`$ approaches the solitary traveling wave $`𝒫(y)𝒫(y+1)`$, with $`y=xL_n^{\mathrm{min}}`$, the front position $`L_n^{\mathrm{min}}`$ given by Eq. (14), and $`𝒫(y)`$ being a solution of Eq. (8). In the critical case, $`p=p_c`$, the minimal length distribution approaches the solitary traveling wave, although the front propagates extremely slowly, as an iterated logarithm. Given that the distribution of the minimal length is the soliton of a finite width, the variation of the minimal length is finite. In other words, $`L_n^{\mathrm{min}}`$ should not vary much from sample to sample.

We have also studied the complementary problem of the statistics of the maximal path. We have found that the maximal length always grows linearly with the height, with the velocity $`v_{\mathrm{max}}`$ again determined via the front selection mechanism. For the bimodal length distribution, the minimal and maximal velocities are connected via the duality relation, $`v_{\mathrm{min}}(p)+v_{\mathrm{max}}(1p)=1`$, which admits extensions to finite-dimensional situations and to arbitrary bounded symmetric length distributions.

We have found that the tail of the minimal length distribution is a double exponential, in contrast with a simple exponential that occurs in most traveling wave problems. This is not very surprising since similar tails were found in statistics of extremes. Note, however, that our exact results exemplify extreme value statistics for correlated random variables (see  for a recent review), while classical results correspond to the case of independent identically distributed random variables.

PLK acknowledges support from NSF, ARO, and CNRS.


text:
# Spectroscopic fingerprints of a surface Mott-Hubbard insulator: the case of SiC(0001)

## Abstract

We discuss the spectroscopic fingerprints that a surface Mott-Hubbard insulator should show at the intra-atomic level. The test case considered is that of the Si-terminated SiC(0001)$`\sqrt{3}\times \sqrt{3}`$ surface, which is known experimentally to be insulating. We argue that, due to the Mott-Hubbard phenomenon, spin unpaired electrons in the Si adatom dangling bonds are expected to give rise to a Si-2p core level spectrum with a characteristic three-peaked structure, as seen experimentally. This structure results from the joint effect of intra-atomic exchange, spatial anisotropy, and spin-orbit coupling. Auger intensities are also discussed.

Keywords: Semi-empirical models and model calculations; Photoemission; Auger electron spectroscopy; Magnetic phenomena; Surface electronic phenomena; Silicon carbide; Insulating surfaces.

The fractional adlayer structures on semiconductors show rich phase diagrams with potential instabilities of charge density wave (CDW) and spin density wave (SDW) type at low temperatures, as well as Mott insulating phases . The (0001)$`\sqrt{3}\times \sqrt{3}`$ surface of hexagonal SiC, (as well as the closely analogous (111)$`\sqrt{3}\times \sqrt{3}`$ of cubic SiC) expected from standard LDA calculations to be a 2D metal with a half-filled narrow band of surface states in the bulk energy gap  – is instead experimentally proven to be an insulator with a rather large ($`2`$ eV) band gap . It has been suggested that this system is a Mott insulator due to the large value of ratio $`U/W`$, where $`U`$ is the Coulomb interaction parameter, of the order of several eV, and $`W`$ is the surface bandwidth, calculated to be about $`0.35`$ eV . Very recent STM data further confirm this picture .

A recent LSDA+U calculation predicts an antiferromagnetic insulating state with $`3\times 3`$ magnetic periodicity and a indirect gap of about $`1.5`$eV, in reasonable agreement with the experimental value . An exchange coupling $`J`$ between neighboring Si adatoms of about $`30K`$ is also calculated. The LSDA+U magnetic state, with spin-collinearity enforced in the calculation, turns out to have a uniform magnetization $`m^z=1/3`$: one spin-down and two spin-up adatom surface bands are completely filled in the $`3\times 3`$ Brillouin zone. Crudely speaking, this corresponds to two adatom dangling-bond orbitals having spin-up electrons and one having spin-down electron in each triangular plaquette. When spin non-collinearity is allowed one finds, e.g., from Hartree-Fock Hubbard model calculations, that the state with spins lying on a plane – in a $`120^o`$ three-sublattice Néel structure, typical of the triangular lattice Heisenberg antiferromagnet – has slightly better energy than the collinear one, but the resulting energy gain is so small, about $`0.5`$ meV/adatom, that unrealistically precise ab-initio calculations would be necessary to decide which of the two is the actual ground state. Moreover, the bands for the two solutions are rather similar, apart from some extra splittings introduced by the collinear magnetic solution .

In reality, both the LSDA+U and the Hartree-Fock Hubbard calculations, being mean-field approaches are at best only cartoons for the actual Mott-Hubbard state, which is insulating not because it is magnetic, but rather due to strong correlations inhibiting electron hopping. Once this is appreciated, one realizes that – even in the absence of any true magnetic long-range order – finite magnetic moments on the Si adatoms should be experimentally detectable up to high temperatures, so long as there is an insulating gap. The low-temperature low-energy collective physics of these unpaired adatom spins is instead governed by a much smaller energy scale, the exchange coupling $`J`$. Is it conceivable that a genuine low-temperature spin long-range order (LRO) would be detectable in this quasi-2D system? It would seem unlikely, since: 1) The estimates for $`J`$ already suggest that LRO, if present, must be there only at very low temperature. Moreover, thermal fluctuations might destroy LRO altogether, if the continuous spin symmetry is not broken by anisotropy effects due to spin-orbit coupling; 2) It is at present unclear which spin ordering will be stabilized by the spin-orbit coupling; 3) If spins lie on a plane, with zero net magnetization, they are presumably hard to detect with LRO sensitive techniques. Therefore, without ruling out totally the possibility of low-temperature LRO, we believe that the SiC(0001)$`\sqrt{3}\times \sqrt{3}`$ surface will most likely be in an overall paramagnetic Mott insulating state, in spite of the existence of an on-site moment, at least down to liquid nitrogen temperature.

It is the aim of the present paper to discuss in some detail what kind of experimental signatures one should look for in a surface Mott-Hubbard insulator, apart from a gap in the photoemission spectrum. As mentioned previously, methods which require magnetic LRO are likely not viable. The situation could be much more promising at the intra-atomic level. We will argue below that the very fact that the system is in a Mott-Hubbard insulating state, independently of the magnetic LRO issue, brings about unpaired spins in the adatom dangling bonds, which in turn determine exchange splittings in the adatom core levels which should be large and detectable. These splittings may in fact have been already observed for SiC(0001) in the surface component of the Si-2p core level photoemission spectrum, although attributed so far to some unspecified charge inequivalence of two types of Si adatoms.

To make things more concrete, let us assume that, due to the Mott-Hubbard phenomenon, the $`3p_z`$-like dangling bond orbital of the Si adatom in the $`\sqrt{3}\times \sqrt{3}`$ SiC(0001) is singly occupied by an unpaired electron, say of spin up. In reality, the Wannier function for the surface state is not a pure $`3p_z`$, since it extends itself down into the first SiC bilayer: the weight $`\alpha `$ of the Si $`3p_z`$ orbital in the Wannier function is about $`0.5`$, the remaining lobe having a negligeable overlap with the Si-adatom core . The $`3s`$ and $`3p_{x,y}`$ orbitals of the Si-adatom are involved in the bonding with first-layer C atoms, and contain spin-paired electrons. Therefore, a Si-adatom has two distinct sources of anisotropy: i) the first comes about because the occupation $`\alpha `$ of the $`3p_z`$ orbital is not equal to the occupation, call it $`\beta `$, of the $`3p_x`$ or $`3p_y`$ orbitals; ii) the second, and most important, source of anisotropy is due to the fact that the electron in the $`3p_z`$ dangling bond is spin unpaired, as opposed to the standard singlet bonding configuration in $`3p_{x,y}`$. These anisotropies break the core level symmetry between $`2p_z`$ and $`2p_{x,y}`$ by an amount which we will argue to be not negligeable compared to spin-orbit. We will now estimate the expected $`2p`$ core level splittings, with and without a spin-orbit coupling, and show that while spatial anisotropy of type i) by itself would result in a very small effect, the spin induced anisotropy ii) has a large influence on the spectrum.

The calculation we perform is very simple, and is based on the Si $`2p3p`$ Slater exchange and Coulomb integrals, obtained from pure atomic Hartree-Fock calculations. Assume the Si adatom valence electrons to be distributed with weights $`\alpha `$ in $`3p_z`$, and $`\beta /2`$ in each $`3p_{x,y}`$ and $``$. Let us consider first the situation in absence of spin-orbit coupling, $`\lambda _{\mathrm{SO}}=0`$. The Hartree-Fock energy of the $`2p_z`$ states can be written as

$$ϵ_{2p_z}=\alpha [K_{2p_z,3p_z}J_{2p_z,3p_z}]+2\frac{\beta }{2}[2K_{2p_z,3p_x}J_{2p_z,3p_x}]+\mathrm{},$$
(1)
where the dots stand for one-body terms and two-body contributions due to other spin-paired occupied orbitals which will exactly cancel out in all the energy differences which we will calculate. Here $`K_{2p_z,3p_z(x)}=(2p_z,3p_{z(x)}|e^2/|𝐫_1𝐫_2||2p_z,3p_{z(x)})`$, and $`J_{2p_z,3p_z(x)}=(2p_z,3p_{z(x)}|e^2/|𝐫_1𝐫_2||3p_{z(x)},2p_z)`$, are the relevant $`2p3p`$ Coulomb and exchange Slater integrals. The factor $`2`$ in front of the $`\beta /2`$ term is due to equal contributions coming from $`3p_y`$ and $`3p_x`$. The factor $`2`$ in front of $`K_{2p_z,3p_x}`$ keeps into account the direct contribution due to both $``$ and $``$ $`3p_{x,y}`$ electrons. (The corresponding exchange term has no factor $`2`$ because only $`3p_{x,y}`$ $``$ electron contribute to the exchange gain of the $`2p_z`$ level.) A similar calculation for $`ϵ_{2p_x}=ϵ_{2p_y}`$ gives

$$ϵ_{2p_x}=\alpha [K_{2p_x,3p_z}J_{2p_x,3p_z}]+\frac{\beta }{2}[2K_{2p_x,3p_x}J_{2p_x,3p_x}]+\frac{\beta }{2}[2K_{2p_x,3p_y}J_{2p_x,3p_y}]+\mathrm{},$$
(2)
which can be reexpressed, due to rotational symmetry, in terms of the previously introduced $`K_{2p_z,3p_z(x)}`$ and $`J_{2p_z,3p_z(x)}`$. We now estimate the Slater integrals assuming pure atomic hydrogen-like $`2p3p`$ orbitals (the effect of intra-atomic hybridization in a full atomic HF calculation would change a bit the numbers: these effects will be estimated later on). Numerically, the two different Coulomb integrals appearing in Eqs. 1-2 are large and comparable ($`K_{2p_z,3p_z}30.5`$ eV, $`K_{2p_z,3p_x}28.2`$ eV), whereas the $`2p_z3p_z`$ exchange integral is smaller ($`J_{2p_z,3p_z}4.3`$ eV), but one order of magnitude larger than the $`2p_z3p_x`$ exchange integral ($`J_{2p_z,3p_x}0.5`$ eV). The large difference between $`J_{2p_z,3p_z}`$ and $`J_{2p_z,3p_x}`$ immediately implies that the exchange splitting between $`2p_z`$ and $`2p_z`$ is rather large

$`ϵ_{2p_z}ϵ_{2p_z}\alpha J_{2p_z,3p_z}2eV`$

whereas the correponding exchange splitting between $`2p_{x/y}`$ and $`2p_{x/y}`$ is one order of magnitude smaller:

$`ϵ_{2p_{x/y}}ϵ_{2p_{x/y}}\alpha J_{2p_z,3p_x}0.2eV`$

We can also calculate, using Eqs. 1-2, the splitting $`\mathrm{\Delta }=ϵ_{2p_{x/y}}ϵ_{2p_z}`$ between the $`2p_z`$ and the $`2p_{x/y}`$ levels,

$$\mathrm{\Delta }=\alpha [(J_{2p_z,3p_z}J_{2p_z,3p_x})(K_{2p_z,3p_z}K_{2p_z,3p_x})]+\frac{\beta }{2}[2(K_{2p_z,3p_z}K_{2p_z,3p_x})(J_{2p_z,3p_z}J_{2p_z,3p_x})].$$
(3)
Notice that the structure of the $`\alpha `$ contribution in Eq. 3 differs from that of the $`\beta `$ term because the $`3p_z`$ is assumed to be spin unpaired (anisotropy ii), above): If we were to assume that the occupancy $`\alpha `$ was due $`\alpha /2`$ $`3p_z`$ and $`\alpha /2`$ $`3p_z`$, the $`\alpha `$ term would be given by $`(\alpha /2)[(J_{2p_z,3p_z}J_{2p_z,3p_x})2(K_{2p_z,3p_z}K_{2p_z,3p_x})]`$ and would exactly cancel off the $`\beta `$ term when $`\alpha =\beta `$ (bulk case). Plugging in the numbers for the Coulomb and exchange integrals we get $`[2(K_{2p_z,3p_z}K_{2p_z,3p_x})(J_{2p_z,3p_z}J_{2p_z,3p_x})][2(2.3eV)3.8eV]=0.6eV`$, and $`[(J_{2p_z,3p_z}J_{2p_z,3p_x})(K_{2p_z,3p_z}K_{2p_z,3p_x})][3.8eV2.3eV]=1.5eV`$, which provide the desidered estimate for $`\mathrm{\Delta }\alpha [1.5eV]+(\beta /2)[0.6eV]`$, i.e., $`\mathrm{\Delta }1eV`$ for $`\alpha 0.5`$ and $`\beta 1`$. The effect of intra-atomic orbital hybridization which modifies the atomic Hartree-Fock orbitals with respect to pure hydrogen-like $`2p3p`$ atomic orbitals, is estimated, from a full Hartree-Fock $`Si^+`$ atomic calculation, to further reduce $`\mathrm{\Delta }`$ and the $`2p_{z/x}`$ exchange splittings by a factor of about $`2`$, so that we arrive at the crude estimate $`\mathrm{\Delta }0.5eV`$, $`ϵ_{2p_z}ϵ_{2p_z}1eV`$, and $`ϵ_{2p_x}ϵ_{2p_x}0.1eV`$. The overall picture for the Si-adatom $`2p`$ spectrum in absence of spin-orbit, but with all the previously described anisotropy effects included, is summarized in Fig. 1(b). When we include the spin-orbit interaction $`\lambda _{\mathrm{SO}}𝐋𝐒`$, the final core hole levels are obtained as in Fig. 2. For $`\lambda _{\mathrm{SO}}=0.4`$ eV (the value extracted from the bulk $`2p_{J=1/2}2p_{J=3/2}`$ experimental core level splitting) we find a roughly three-peaked structure, the broad central peak four times as strong as each side peak. This offers an alternative explanation of the experimental lineshape (see Fig. 2) in terms of a single exchange-split multiplet, rather than two chemically inequivalent sites $`S_1`$ and $`S_2`$  whose existence is otherwise not supported. Incidentally, a similar three-peaked structure for the Si-2p core level spectrum has been experimentally observed for a sub-monolayer Si segregated at the surface of iron.

Additional direct experimental evidence for the Mott-Hubbard state could be obtained by a careful study of adatom ($`3p_z`$, $`3p_z`$) Auger spectral intensities, which should be easily singled out owing to the large gaps. The probability of double occupancy of the adatom dangling bond orbital should be almost completely suppressed, dropping from the band value of $`1/4`$, to a value of order $`t/U`$ – where $`t`$ is the hopping matrix element and $`U`$ is the effective on-site Coulomb repulsion – which is one order of magnitude smaller. Even considering that only half of the orbital is adatom $`3p_z`$, this surface should show a ($`3p_z`$, $`3p_z`$) Auger intensity which, by comparison with the remaining $`(3p,3p)`$ values, is anomalously small, as a proof of its Mott-Hubbard state.

Summarizing, we have argued that the Mott-Hubbard insulating surface state in SiC(0001) should leave clear and detectable intra-atomic fingerprints. The Si $`2p`$ core hole will exhibit a $`6`$-fold multiplet resulting from the joint effect of intra-atomic exchange, asymmetry, and spin-orbit, with a characteristic three-peaked structure. Auger intensities should also be affected. More experimental and theoretical effort is clearly called for to check these strong correlations and related magnetic effects, possibly the first of this magnitude to be suggested for an sp-bonded, valence semiconductor surface.

We thank S. Modesti for useful discussions. Work at SISSA was partially supported by INFM through PRA LOTUS and HTSC, by MURST through COFIN97, and by the EU, through contract FULPROP ERBFMRXCT970155.

Figures


text:
# Fractal Electromagnetic Showers

## 1 Introduction

One of the most serious problems in the analysis of cosmic ray data is the complex and time-consuming nature of the codes used for shower simulation. In order to try to capture the detailed physics of the processes involved, it is customary to directly simulate  the multiplicative branching process whereby an initial particle gives rise to two or more secondary particles, each of which, in turn, initiates what is essentially its own shower, albeit now at lower energy.

Such a process can give rise to large fluctuations, and the final distributions of ground particles and their energies (as well as the longitudinal distribution of the shower as a whole) are difficult to model with simple parametrizations unless one is happy to settle for a description of the mean behaviour of the shower and forego knowledge of the fluctuations. Indeed, this is the leading reason that so much Monte Carlo time must be used for shower simulations: there are no simple analytical forms for the relevant distributions which can describe the fluctuations. The issue is a pressing one for experiments collecting large amounts of data which may be difficult to compare against theory in any form other than a large number of simulated events.

Here we report on the observation that electromagnetic showers display self-similar behaviour which can be described by a multifractal geometry and describe first steps towards formalizing this concept. Our eventual goal is to describe showers in terms of what we argue here is the relevant geometry: not one of smooth functions, but one which allows for irregular geometries which are better described in terms of fractals. We consider here only electromagnetic showers, but plan to study hadronic showers in future work.

## 2 Self-Similarity in Electromagnetic Showers

The idea that an electromagnetic shower should, in some sense, be a fractal is almost obvious. It is generated recursively from the two processes”

1. pair creation: $`\gamma \mathrm{e}^+\mathrm{e}^{}`$ in the electric field of a nucleus and;
2. Bremsstrahlung: $`\mathrm{e}^\pm \mathrm{e}^\pm \gamma `$ as an electron or positron is deflected by the electric field of a nucleus

This is illustrated in figure 1 which shows the particles making up a shower produced by a 100 GeV electron entering a block of aluminum 150 cm long (radiation length 8.9 cm) as simulated using the geant4 program .

Each final state particle from an interaction effectively initiates its own electromagnetic shower, and each process has a similar cross section to occur in matter. As long as the energies involved are large compared with the energy required to create an electron-positron pair (and thus also large compared to atomic processes such as ionization), each step is much the same as the one before it, but at a reduced energy.

Figure 2 shows a slice through the block right at the far end with the point of intersection of each particle with the slice shown as a black dot whose radius is independent of energy. Here one clearly sees the shower core, with a diminishing density of particles with distance from the centre.

## 3 Fractals and Multifractals

There are many ways to characterize self-similar objects, but the most common and well-known way is in terms of fractal dimensions. There are many different concepts of fractal dimension which are useful, and perhaps the most obvious is that of mass dimension, $`D_M`$. The idea here is to see how the total energy $`E_{TOT}(R)`$ (considered now as a sort of weight) within a disk of radius $`R`$ varies as $`R`$ is changed. If the distribution were one-dimensional (a line of uniform energy deposited in the plane), one would find

$$E_{TOT}(R)R^1$$
(1)
and one would take the exponent in the foregoing equation to be the dimension of the distribution.

If the energy were uniformly distributed over the whole plane, one would find

$$E_{TOT}(R)R^2$$
(2)
and conclude again that the exponent in the scaling law for the energy should be interpreted as the dimension of the distribution.

In the event that a scaling law of the form $`E_{TOT}(R)R^{D_M}`$ holds for a non-integer $`D_M`$, we call $`D_M`$ the “fractal mass dimension”. A plot of $`\mathrm{log}(E)`$ as a function of $`\mathrm{log}(R)`$ will then have a slope in the limit of small $`R`$ which is $`D_M`$.

Two points are important to keep in mind here: first that there are no true fractals in nature as there are always some smallest and largest value for variables in the problem beyond which scaling behaviour does not hold, and second that one must be careful to watch for systematic effects which can bias estimates of the dimension. Systematic effects which we have had to be wary of include the fact that early in the shower development the central core can contain particles which carry a large fraction of the initial energy and give the radial energy distribution a spike at small $`R`$ which does not correspond to scaling behaviour.

In the case of the electromagnetic shower with the slice taken at the end of the shower at 150 cm, we look at the summed energy (scaled so that the total energy is 1) as a function of the fraction of the radius out (scaled so that the maximum radius is 1). This quantity we denote as $`I(R|1)`$ for reasons which will become clear later in the text. Plotting logarithms against logarithms (base 10), we find the distribution shown in figure 3.

The first thing to notice is that the curve is reasonably approximated by a straight line at small radii. The second thing to notice is that the whole curve is not a straight line. At large radii we start to reach the physical boundaries of the shower and cannot expect scaling to hold.

In fact, even at very small radii, there is some anomalous structure which can be traced to the effects of very energetic particles very close to the core, which give an additional spike of energy to the distribution which cannot be expected to be a part of any overall scaling behaviour. This effect is more pronounced earlier in the shower.

The scaling properties of the shower are thus different in different parts of the plane, and in order to quantify this further, we study the scaling behaviour of cumulative moments of the energy distribution defined for $`q>0`$ by

$$I(R|q)=\frac{_{r<R}E_i^q}{_{\mathrm{all}i}E_i^q}$$
(3)
where $`E_i`$ are the energies contained in a disk going out to radius $`R`$ and the sum is taken over all particles within a distance $`r<R`$. What units are used is not important as we are only interested in the average scaling behaviour of the curves at small $`R0`$. (As discussed earlier in the text, the region of very small $`R`$ should be avoided for physical reasons, and we will avoid the subtleties of precise numerical analyses in this short communication.) For graphical purposes here, $`R`$ is normalized so that the particle with the largest radial distance out is at $`R=1`$ and the moments are defined so that their value at maximum radius is unity. We can then introduce an infinite family of fractal dimensions $`D_q`$ defined for $`q>0`$ by

$$D_q=\underset{R0}{lim}\frac{1}{q}\frac{\mathrm{log}I(R|q)}{\mathrm{log}R}$$
(4)
with the understanding that the limit must still lie in the scaling region in physical examples.

Figure 4 shows the scaling behaviour of moments of the electromagnetic shower corresponding to how the sums of the squares and cubes of the energy grow with distance.

For a homogeneous and uniform fractal structure we expect the $`D_q`$ to be equal. If not, then we describe the distribution as multifractal in that it requires more than one fractal dimension in order to characterize it. The associated $`D_q`$ for small $`q`$ estimated from finite differences in the scaling region are all approximately equal within the errors in the data here and approximately unity, suggesting a good degree of homogeneity. It is important to keep in mind that the results in this paper are presented for a full, realistic GEANT simulation, and include ionization, delta-ray, and other soft processes, so some care is needed in interpreting the results as if they corresponded to a pure electromagnetic shower generated only by pair creation and Bremsstrahlung (which is, of course, not realizable in nature).

The definition of fractal dimensions can also be continued to $`q0`$, but this has some subtleties involved with the fact that as $`q\mathrm{}`$ the highest energy particles contribute most, while as $`q\mathrm{}`$ the lower energy ones dominate. In particular, some care must be used with the $`D_q`$ for $`q<0`$ as they give high weights to softer particles which are not part of the hard shower process. These matters, as well as more precise results on dimensions including energy and material dependence will be presented elsewhere.

## 4 Further Work

Clearly space limitations make it impossible to cover the material as completely as one would like, but several points concerning work not discussed here are worth making. First of all, we expect fractal behaviour in all three dimensions, and in this discussion we have neglected the longitudinal scaling behaviour, where the full shower is made of many scaled and translated showers superimposed along the shower axis. In addition, there are clearly angular correlations and fluctuations, and studies can be made at a given fixed radius of the scaling behaviour of the shower as a function of the angular coordinate which we have integrated out in this discussion. The relation of these ideas to the concept of intermittency, especially as studied in hadronic jets has not escaped our notice and is currently under investigation.

One of the main goals of this work is to better understand the geometry of electromagnetic (and other) showers in order to try to parametrize them by the appropriate non-smooth basis functions, such as wavelets. Such a parametrization should allow the fast generation of showers without the attendant loss of information concerning large fluctuations which has so far been handled only by the use of enormous computational resources.

## 5 Acknowledgements

We would like to thank the US National Science Foundation and CONICET, Argentina for support. We would also like to thank our collaborators on the Pierre Auger Project, as well as on L3 and CMS for useful discussions on electromagnetic calorimetry.


text:
# Aging of the Zero-Field-Cooled Magnetization in Ising Spin Glasses: Experiment and Numerical Simulation

\[

## Abstract

A new protocol of the zero-field-cooled (ZFC) magnetization process is studied experimentally on an Ising spin-glass (SG) Fe0.50Mn0.50TiO3 and numerically on the Edwards-Anderson Ising SG model. Although the time scales differ very much between the experiment and the simulation, the behavior of the ZFC magnetization observed in the two systems can be interpreted by means of a common scaling expression based on the droplet picture. The results strongly suggest that the SG coherence length, or the mean size of droplet excitations, involved even in the experimental ZFC process, is about a hundred lattice distances or less.

PACS numbers: 75.10.Nr, 75.40.Gb, 75.50.Lk

\]

Spin glasses generally exhibit puzzling aging behavior such that their magnetic response in the spin-glass (SG) phase below its transition temperature $`T_\mathrm{c}`$ depends crucially on their thermal history within the phase. For example, an isothermal zero-field-cooled (ZFC) magnetization, $`m_{\mathrm{ZFC},T}(t_a,t_\mathrm{w})`$, is defined as the magnetization induced by a magnetic field $`h`$ which is applied after a system is rapidly quenched from above $`T_\mathrm{c}`$ to a temperature $`T`$ below $`T_\mathrm{c}`$ and is aged for a waiting time of $`t_\mathrm{w}`$, and $`t_a=\tau +t_\mathrm{w}`$ with $`\tau `$ being the time that the system ages at $`T`$ under $`h`$. With increasing $`\tau `$, $`m_{\mathrm{ZFC},T}(t_a,t_\mathrm{w})`$ grows slowly. In contrast to ideal equilibrium situations, however, it depends not only on $`\tau `$ but also on the waiting time $`t_\mathrm{w}`$, namely, the larger $`t_\mathrm{w}`$, the slower the growth of $`m_{\mathrm{ZFC},T}(t_a,t_\mathrm{w})`$. This was the first observation of aging phenomena in spin glasses . Since then most of the studies have focused on its $`t_\mathrm{w}`$-dependence.

In the present work, on the other hand, we have studied the ZFC magnetization focusing on its dependence on thermal processes after the field $`h`$ is applied. We denote it by $`m_{\mathrm{ZFC},T[\tau ]}(t_a,t_\mathrm{w})`$ where $`t_\mathrm{w}`$ symbolically represents a common thermal history before $`h`$ is applied, while $`T[\tau ]`$ does various heating processes after $`h`$ is switched on, e.g., with a fixed heating rate but with an intermittent stop(s) at a different temperature(s). We have carried out this new protocol of the ZFC magnetization measurement both experimentally on an Ising SG material Fe0.50Mn0.50TiO3, hereafter referred to as FeMnTiO, and numerically on the 3D Gaussian Edwards-Anderson (EA) Ising SG model. Quite interestingly, we have found that $`m_{\mathrm{ZFC},T[\tau ]}(t_a,t_\mathrm{w})`$ measured in the two systems are well interpreted by a unified way, or more explicitly, they are described by a common scaling expression which can be derived based on the droplet picture for the SG phase . The purpose of the present letter is to demonstrate this finding and to discuss its implications on the aging dynamics in spin glasses.

According to the droplet picture, the mean size, $`R(t_\mathrm{w})`$, of domains of an equilibrium state grows continuously when a spin glass relaxes with increasing $`t_\mathrm{w}`$. In fact its growth law in isothermal aging at a temperature $`T`$, denoted by $`R_T(t_\mathrm{w})`$, has been numerically studied by the various groups , and within the numerical accuracy all the results are consistent with each other. Here we quote our result,

$$R_T(t_\mathrm{w})(t_\mathrm{w}/\tau _0)^{b_\mathrm{f}T/T_\mathrm{c}},$$
(1)
where $`\tau _0`$ is the microscopic unit of time which is one Monte Carlo step (MCS) in the simulation, and $`b_\mathrm{f}0.17T_\mathrm{c}/J`$  $`0.153`$ with $`T_\mathrm{c}0.9J`$ , $`J`$ being the variance of the interactions. Furthermore we have shown  that the aging properties of the spin auto-correlation function are consistently explained when the same functional form is adopted for the maximum size of droplets which can be activated within a time scale of $`\tau `$, $`L_T(\tau )`$, i.e.,

$$L_T(\tau )(\tau /\tau _0)^{b_\mathrm{f}T/T_\mathrm{c}}$$
(2)
with the same constant $`b_\mathrm{f}`$. For an aging process in which the temperature is changed as in our new protocol, the maximum size of droplet excitations at time $`t_a=\tau +t_\mathrm{w}`$, denoted by $`L(\tau )`$, is evaluated by the successive use of eq.(2) as described by eqs.(4) and (5) below.

In terms of $`L(\tau )`$ introduced just above the common scaling expression we have found is written as

$$h^1m_{\mathrm{ZFC},T[\tau ]}(t_a,t_\mathrm{w})=\chi _{\mathrm{ZFC}}(L(\tau ),R(t_\mathrm{w})),$$
(3)
where $`R(t_\mathrm{w})`$ is the mean domain size at the instance when the field $`h`$ is switched on. For isothermal processes under a field within the linear response regime the droplet theory  has already predicted a similar scaling form, i.e., $`m_{\mathrm{ZFC},T}(t_a,t_\mathrm{w})`$/h is a function of $`R_T(t_\mathrm{w})`$ and $`L_T(\tau )`$. Equation (3) extends it to aging processes which involve temperature changes. Also, by the low temperature expansion of the droplet theory, $`m_{\mathrm{ZFC},T}(t_a,t_\mathrm{w})`$/h in the quasi-equilibrium regime with $`R_T(t_\mathrm{w})L_T(\tau )`$ is shown  not to depend explicitly on $`T`$. Equation (3) shares the same property but for the whole range of $`\tau `$. In relation to this it is worth pointing out that the field-cooled (FC) magnetization, $`m_{\mathrm{FC}}(T)`$, which is supposed to be lim$`{}_{t_a\mathrm{}}{}^{}m_{\mathrm{ZFC},T}^{}(t_a,t_\mathrm{w})`$, is nearly independent of $`T`$ for the Ising spin glasses such as the FeMnTiO and the EA Ising model.

The unified description of the ZFC magnetizations means that eqs.(3) and (2) are quantitatively consistent with the measured results on the FeMnTiO spin glass when $`\tau _010^{12}`$s and the same value of $`b_\mathrm{f}`$ as that of the EA-SG model are used. It is, to our knowledge, the first case in the study on aging in spin glasses that such a quantitative agreement is obtained in the direct comparison of the experiment with the numerical simulation.

We have measured the induced magnetization of the single crystal FeMnTiO sample with $`T_\mathrm{c}=`$21.5K  by applying the field parallel to the hexagonal $`c`$-axis of the sample (see  for further experimental details). The numerical simulations have been done on the 3D EA model with nearest neighbour interactions of mean zero and variance $`J`$ which is the unit of energy. We have used the heat-bath dynamics and have simulated several hundred samples with a linear size $`L=24`$ to make the error bars less than the size of the symbol in each plot. It is noted that the units of time and temperature used in the present work are respectively second and Kelvin in the experiment, and MCS and $`J`$, which we put unity below, in the simulation. Our processes measuring magnetizations in the experiment \[simulation\] are as follows. The system was first cooled from 50K to 5K \[2.5$`J`$ to 0.3$`J`$\] at a rate of 1.0 K/min \[0.05$`J`$/10MCS\] under zero field. Once the lowest temperature was reached the field of 100 Oe \[0.2$`J`$\] was applied and the sample was heated stepwise with a step of height 0.2K \[0.01$`J`$\] and length 100s \[100MCS\]. The last 13s \[10MCS\] of the latter were used to measure $`m_{\mathrm{ZFC},T[\tau ]}(t_a,t_\mathrm{w})`$ which we regard as the ZFC magnetization at that temperature, $`m_{\mathrm{ZFC}}(T)`$. At 50K, in the experiment, we turned to decrease the temperature and measured the FC magnetization, $`m_{\mathrm{FC}}(T)`$, on cooling with the same rate and in the same field. In the simulation $`m_{\mathrm{FC}}(T)`$ was measured on cooling the sample under $`h=0.2J`$ from 2.5$`J`$ down to 0.3$`J`$ with the same rate as for $`m_{\mathrm{ZFC}}(T)`$.

The ZFC and FC magnetizations as well as the one with an intermittent stop of the real and model Ising spin glasses, shown in Fig. 1, look qualitatively quite similar to each other. Particularly, $`m_{\mathrm{FC}}(T)`$ of both systems become nearly flat at $`T`$ sufficiently lower than $`T_\mathrm{c}`$. We can also observe the following common feature in $`m_{\mathrm{ZFC}}(T)`$ with an intermittent stop. During the stop $`m_{\mathrm{ZFC}}(T)`$ increases, and when the heating is restarted, it continues to increase. However it increases significantly much slowly than $`m_{\mathrm{ZFC}}(T)`$ without the stop, and is caught up by the latter at a certain $`T`$ below $`T_\mathrm{c}`$.

In the inset \[main frame\] of Fig. 2 we demonstrate the most interesting observation in our new protocol. It is the comparison of the explicit time evolutions of $`m_{\mathrm{ZFC},T[\tau ]}(t_a,t_\mathrm{w})`$ of the FeMnTiO spin glass \[the EA SG model\] with one intermittent stop at $`T=15.0`$K \[0.55$`J`$\] and two stops at $`T=14.6`$K by 1 hour and at $`T=15.0`$K \[$`T=0.50J`$ by 3000 MCS and at $`T=0.55J`$\]. Here we rewrite them as $`m_{\mathrm{ZFC}}(t)`$ since we are interested in their time evolution. The origin of the time $`t`$ is set at the instance when the system is heated up to 15.0K \[0.55$`J`$\] for the 1-stop process and to 14.6K \[0.50$`J`$\] for the 2-stop process. It is seen that, when the 1-stop curve (crosses) is shifted to the right by 2250 s \[2100 MCS\] (open circles), it almost perfectly collapses to the branch of $`T=15.0`$K \[0.55$`J`$\] of the 2-stop curve (solid circles). In other words, once $`m_{\mathrm{ZFC},T[\tau ]}(t_a,t_\mathrm{w})`$ of two different thermal processes coincide with each other, their subsequent evolutions under a common condition of $`T`$ and $`h`$ are also identical. Keeping in mind that these processes have a common $`R(t_\mathrm{w})`$, or a common thermal history before $`h`$ is applied, the results on both the real and model spin glasses are consistent with eq.(3) derived from the droplet picture.

To ascertain the scaling expression of eq.(3), and to demonstrate that the $`T`$-dependence of the scaling function $`\chi _{\mathrm{ZFC}}(L(\tau ),R(t_\mathrm{w}))`$ comes out only through that of $`L(\tau )`$, we have performed the simulation of isothermal ZFC magnetization processes at various temperatures but by choosing the waiting time $`t_\mathrm{w}`$ after the instantaneous quench such that $`R_T(t_\mathrm{w})`$ given by eq.(1) become the same for all temperatures. The results are shown in Fig. 3 where $`m_{\mathrm{ZFC},T}(t_a,t_\mathrm{w})`$ are plotted against the time $`\tau `$ after the field $`h=0.2`$ is switched on. The chosen $`R_T(t_\mathrm{w})`$ is $`2.1`$ and the corresponding $`t_\mathrm{w}`$ at each temperature are indicated in the figure.

The evolutions of $`m_{\mathrm{ZFC},T}(t_a,t_\mathrm{w})`$ in Fig. 3 depend on $`T`$ except at $`\tau 1`$. When, however, the same data are plotted against $`L_T(\tau )`$ evaluated by eq.(2), they are nicely superposed upon each other as shown in the inset of Fig. 3. This result strongly supports the scaling expression of eq.(3) in the following two respects: the thermal history before the field is applied is in fact renormalized to $`R(t_\mathrm{w})`$, and the size $`L_T(\tau )`$ of droplet excitations which are dominant in the time scale $`\tau `$ after the field is applied is properly given by eq.(2). We also note that, within the accuracy of the present numerical analysis, the scaling behavior is ascertained not only in the quasi-equilibrium regime, where $`L_T(\tau )<R_T(t_\mathrm{w})`$ or $`\tau t_\mathrm{w}`$, as predicted by the droplet theory , but also in the out-of-equilibrium, or aging regime, where $`L_T(\tau )>R_T(t_\mathrm{w})`$ or $`\tau t_\mathrm{w}`$. In the latter regime $`L_T(\tau )`$ plays a role of the mean domain size at time $`t_a`$ in place of $`R_T(t_\mathrm{w})`$.

Lastly let us extend the scaling analysis done above to our new experimental protocol in which the temperature is changed. For this purpose it is important to note that, so long as the system is in the SG phase, the mean domain size $`R(t_\mathrm{w})`$ of the EA model analyzed numerically ever grows continuously even when the temperature is changed discontinuously . It is then natural for us to expect that this is also the case for $`L(\tau )`$ under a fixed field $`h`$, so that we may successively use eq.(2) as follows. Suppose that temperature is changed as $`T_i=T_{i1}+\mathrm{\Delta }T_i`$ and a period of time $`d\tau _i`$ is spent at each temperature $`T_i`$. We assume that the length scale achieved after the step at $`T_i`$ is given as

$$L_i=(\tau _i^{\mathrm{eff}}+d\tau _i)^{b_\mathrm{f}T_i/T_\mathrm{c}},$$
(4)
where $`\tau _i^{\mathrm{eff}}`$ is an effective elapsed time to take into account the thermal history up to the instance just before the temperature is changed to $`T_i`$. It is the time needed to achieve the length scale $`L_{i1}`$ reached after the $`(i1)`$-th step by isothermal aging at $`T_i`$,

$$(\tau _i^{\mathrm{eff}})^{b_\mathrm{f}T_i/T_\mathrm{c}}=L_{i1}.$$
(5)
Combining these we obtain a recursion formula for $`L_i`$ and/or $`\tau _i^{\mathrm{eff}}`$ which is solved for any given temperature process with $`L(\tau =0)=0`$.

We have first applied the above-mentioned analysis to $`m_{\mathrm{ZFC},T[\tau ]}(t_a,t_\mathrm{w})`$ of the EA model in the processes with $`T`$-changes shown in Figs. 1 and 2. The results are shown in Fig. 4. All $`m_{\mathrm{ZFC},T[\tau ]}(t_a,t_\mathrm{w})`$ nicely lie on a universal curve, indicating that the scaling of eq.(3), combined with the above analysis on $`L(\tau )`$, does work well. In the figure we also draw one of the curves in the isothermal process shown in the inset of Fig. 3. It has the different thermal history before the field is applied, i.e., the different $`R(t_\mathrm{w})`$, from those of the $`T`$-change process. It significantly deviates from the universal curve for the latter as expected from eq.(3).

Surprisingly, the $`m_{\mathrm{ZFC},T[\tau ]}(t_a,t_\mathrm{w})`$ curves measured in the FeMnTiO spin glass and shown in the inset of Figs. 1 and 2 obey the same scaling law as those of the EA model. Indeed, their scaling plots demonstrated in Fig. 5 are obtained by making use of eq.(2) with the identical exponent of $`b_\mathrm{f}T/T_\mathrm{c}`$ but with $`\tau _0=10^9`$ s and $`10^{12}`$ s, instead of $`\tau _0=1`$ MCS for the EA model. The two portions with the dense data points at around $`m_{\mathrm{ZFC},T[\tau ]}(t_a,t_\mathrm{w})0.14`$ and $`0.20`$ represent the results with intermittent stops at 9.6K, 10.0K and at 14.6K, 15.0K, respectively. At a glance we see that the scaling with $`\tau _0=10^{12}`$ s is the better. We have not performed further analysis on the sensitivity of the collapse of the curves to the parameter values, because the evolution law of eq.(2) or (1) already involves certain uncertainty. For example, $`R_T(t_\mathrm{w})`$ extracted by the ac susceptibility in a sufficiently wide frequency range is better expressed by a function of ln$`t_\mathrm{w}`$ . Furthermore we are faced to similar but more serious problems in the temperature range close to $`T_\mathrm{c}`$, such as the crossover to the critical regime , whose analysis has been omitted in the present work. However we consider that the scaling results obtained above are enough for us to claim that the droplet picture expressed by eq.(3) is consistent with the aging behavior of the ZFC magnetization in the FeMnTiO Ising spin glass, and that the length scale involved in the process is about a hundred lattice distances or less. They also imply the continuous growth of $`L(\tau )`$ and so $`R(t_\mathrm{w})`$, without suffering from the chaotic effect , at least in the present experimental protocol.

Another interesting result of the present work is that the evolution law of $`L_T(\tau )`$ given by eq.(2), or the growth law of $`R_T(t_\mathrm{w})`$ given by eq.(1), almost coincides with the growth law of the SG coherence length, $`\xi _T(t)`$, which Joh et al.  has extracted from the thermoremanent magnetization (TRM) of the spin glasses CuMn, AgMn and CdIn0.3Cr1.7S4, though none of them is an Ising spin glass. Their results imply that the evolution law of $`L_T(\tau )`$ depends on the field at least in the the crossover region between the quasi-equilibrium and the aging regimes. In the present work, however, such a field effect has been neglected regarding it as a higher order effect.

To conclude we have demonstrated that the ZFC magnetizations measured in both real and model Ising spin glasses are described by a common scaling expression which is consistent with the droplet picture on the SG aging dynamics. The result indicates that the numerical simulation on the model spin glass supplies us not only qualitative but also quantitative information on aging phenomena in real spin glasses.

This work is supported by Grant-in-Aids for International Scientific Research Program (#10044064) and for Scientific Research Program (#10640362), from the Ministry of Education, Science, Sports and Culture, Japan. The present simulation has been performed by the computing facilities at Supercomputer Center, Institute for Solid State Physics, University of Tokyo.