text:
# Driving forces for Ag-induced periodic faceting of vicinal Cu(111)

## Abstract

Adsorption of submonolayer amounts of Ag on vicinal Cu(111) induces periodic faceting. The equilibrium structure is characterized by Ag-covered facets that alternate with clean Cu stripes. In the atomic scale, the driving force is the matching of Ag(111)-like packed rows with Cu(111) terraces underneath. This determines the preference for the facet orientation and the evolution of different phases as a function of coverage. Both Cu and Ag stripe widths can be varied smoothly in the 3-30 nm range by tuning Ag coverage, allowing to test theoretical predictions of elastic theories.

Self-organization on crystal surfaces is a promising alternative for growing uniform nanostructures with regular sizes and spacings Pohl1Plass1Noetzel1Brune1 ; Shchukin1 . In the case of one-dimensional nanostructures one can use the regular step arrays of vicinal surfaces as templates. Ideally atomic steps are preferential adsorption sites, and wires and stripes can be grown by step-flow. However, if the terrace width increases to a few nanometers, single steps tend to meander and their spacing becomes irregular. Yet vicinal surfaces offer another alternative at this mesoscopic scale, i.e. using adsorbates to induce the formation of periodic facets, which are much stiffer than individual steps. In this case the miscut angle as well as the adatom concentration can be tuned to obtain a variety of stripe patterns that are potential templates for further nanostructure growth Seehofer1Aoki1Minoda1 ; Himpsel1 ; Minoda2Horn1 ; Bachmann1 ; Foelsch1 . The fundamental problem is controlling this complex process in order to predict and obtain useful template structures. Therefore it is important to understand the microscopic mechanisms that determine the stripe pattern formation.

In the atomic scale, the driving force of the adsorbate-induced faceting appears to be the lattice matching of the adsorbate layer with a particular substrate orientation. The high step density of the vicinal surface facilitates the step bunching and the formation of facets with better lattice matching in the direction perpendicular to the steps. Thus at a temperature where both adatom and surface step mobilities are high enough the system breaks up into two phases, namely well-matched adsorbate-covered facets and adsorbate-free stepped stripes. A widely studied case is Au on vicinal Si (v-Si). Here a variety of overlayer reconstructions match to different Si crystal planes, leading to a family of striped patterns made of low symmetry facets and wide (111) or (100) terraces with different Au concentrations Seehofer1Aoki1Minoda1 ; Himpsel1 ; Minoda2Horn1 . The final morphology depends on the total amount of Au deposited, the temperature and the vicinal angle of the clean surface. Furthermore, one can tune miscut angle and Au coverage to obtain a well-matched single phase system Himpsel1 . More recently adsorbate-induced faceting and nanostructuration has been observed for Ag and NaCl on v-Cu(111) Bachmann1 ; Foelsch1 . In these cases the lattice matching is good between Ag(111) and Na(100) packed layers and Cu(112) and Cu(223) planes, respectively. This induces stripe patterns of covered and uncovered facets, as well as changes in the facet orientation as a function of coverage.

The important question concerns the mesoscopic scale, i.e. the nanostructure periodicity. Using the elastic theory of continuous media for two-phase systems Marchenko has shown that the difference in stress between the two phases produces long-range interactions, leading to periodic domain patterns Marchenko1 . Although this theory explains well the spontaneous faceting of clean surfaces Rousset1 , or the coverage-dependent periodicity of adsorbate patterns on flat surfaces Vanderbilt1 , it is an open question whether it can be applied for adsorbate-induced faceting. On one hand, because the growth of strained layers might bring in other types of elastic relaxations, thereby suppressing the characteristic wavelength Tersoff1 . On the other hand, the adsorbate-induced faceting is a more complex process that involves huge mass transport and restructuration of two different chemical species, and hence it is more likely to be kinetically restricted. This is the case of Au/v-Si(100) and NaCl/v-Cu(111) Foelsch1 ; Meyer1 . Here we analyze the Ag-induced faceting of v-Cu(111). We show that this self-assembled striped structure can be tuned in the 3-30 nm range by simply controlling Ag coverage Bachman2 . The morphology of the system is stable over a wide temperature range. At the atomic scale we show that faceting is driven by the tendency of Ag to form (111)-like packed layers with minimum stress. At the mesoscopic scale we are able of proving, for the first time in adsorbate-induced faceting, the coverage-dependent periodicity predicted by the continuum elastic theory at thermal equilibrium.

The substrate is v-Cu(111) with $`11.8^{}`$ miscut about the \[$`\overline{1}\overline{1}`$2\] direction. The surface preparation is described elsewhere Bachmann1 . Ag deposition was done with the substrate held at 300K. The regular stripe structure is produced after subsequent annealing to 420K. Further annealing to 700K does not produce any appreciable change in Scanning Tunneling Microcopy (STM) images or Low Electron Diffraction (LEED) patterns, strongly suggesting thermal equilibrium. In order to have a continuous variation of the coverage, Ag is deposited as a wedge. The coverage is calculated from STM images using 3000 Å2 frames and assuming that Ag covered areas consist of packed, 1 monolayer (ML) thick patches. Such assumption is consistent with the structural analysis of Ag-covered areas and also with Auger Electron Spectrocopy measurements Bachmann1 .

In Figures 1 (a) to 1 (e) we study the morphology of the system as a function of Ag coverage by means of STM and LEED. Fig. 2 schematically depicts the side view of the faceted structure deduced from the analysis. The clean surface shows a regular array of monoatomic, {100}-like steps along the \[$`\overline{1}`$10\] direction, giving rise to the characteristic spot splitting in LEED. The average terrace width measured with both LEED and STM is $`d`$=10.2 Å, i.e. the surface mostly shows (223)-like terraces (4$`\frac{2}{3}`$ atomic rows per terrace, $`d`$=10.5 Å) with a few (335)-like terraces (3$`\frac{2}{3}`$ atomic rows, $`d`$=8.4 Å). As shown in Figs. 1(b)-(d), upon Ag adsorption the system undergoes faceting. Below 1 ML we distinguish three different regimes, which differ on the crystallographic orientation of the Ag-covered facet. These are the $`A`$-regime up to $``$0.52 ML, the $`B`$-regime between $``$0.52 ML and $``$0.76 ML, and the $`C`$-regime from $``$0.76 ML to 1 ML.

Figs. 1 (b)-(c) correspond to the $`A`$-regime, which is characterized by a fairly regular hill-and-valley structure of Ag-covered facets oriented in the (112) direction and clean Cu stepped areas. The inset of Fig. 1 (b) shows a closer view of the surface, where we can observe the Moire pattern in the Ag-covered facet, and monoatomic steps at clean Cu stripes with their typical frizzy edges. In the LEED pattern we have both the spot structure from the Ag Moire and the splitting from the step array in the Cu stripe. As observed in Figs. 2 (b)-(c) and depicted in Fig. 2, increasing the Ag coverage in $`A`$-regime results in a higher density of Ag-covered stripes and a reduction of the step density at clean Cu bands. The latter is due to the fact that Cu(112) facets ($`d`$=6.25 Å, 2$`\frac{2}{3}`$ atomic rows per terrace) have higher step density than the average surface, thus their relative growth requires additional steps from Cu stripes. Such effective step removal from Cu stripes is nicely followed in the LEED pattern in Fig. 1 (e). The Cu splitting smoothly shrinks as a function of Ag coverage, and both spots merge into a single one with 0.52 ML, when $`A`$-regime saturates note0 .

The $`B`$-regime from 0.52 ML to 0.76 ML is characterized by a smooth transition in the Ag-covered facet orientation from (112) to (335). The latter is closer to the surface normal, allowing a higher Ag saturation. (335)-covered areas develop directly from the center of preexisting (112) facets or arise by melting two contiguous (112) facets. In both cases, the (335) orientation characterizes the center of the Ag-covered stripes, whereas the boundaries preserve the (112) orientation. For 0.72 ML in Fig. 1 (d) the $`B`$-regime is almost saturated, i.e. we have flat Cu(111) terraces and mostly (335)-oriented Ag stripes, also displaying a Moire structure. During the $`C`$-regime, from 0.76 ML to 1ML, (223) facets develop. With 1 ML we observe a close-packed Ag layer wetting a random distribution of (335)- and (223)-like terraces with the same step density than the starting surface.

In the atomic scale, the driving force for the Ag-induced faceting of the stepped Cu surface is the matching between Ag(111)-like packed rows along the parallel direction of the steps and the Cu facet underneath. Lattice-matching has been also proposed as the microscopic mechanism for the NaCl-induced faceting of Cu(112) Foelsch1 . Our system generally displays a clear tendency to avoid adsorption on wide (111) terraces, which leads to a large mismatch Meunier1 . On the other hand, as deduced from Figs. 1 and 2, there is a preference for (112)-oriented facets at low coverages, in spite that this requires a larger mass transport. This can be naively explained from the width of the Cu terrace required for a good matching in the perpendicular direction. Five Ag packed rows (2.498 Å wide each) fit to two 6.25 Å wide terraces in Cu(112) with only 0.12% mismatch in the direction perpendicular to the steps, whereas Cu(335) requires two 8.37 Å wide terraces to accommodate seven Ag rows with 4.47% mismatch. This preference for (112) facets is also supported by the microscopic analysis of (112) and (335) stripes shown in Fig. 3. The atomic models in the right panels of Fig. 3 reproduce the Moire patterns of the left panels in the simplest way. The open circles represent Ag-close packed layers, the small dots the Cu(112) and Cu(335) unit cell underneath, and the shades the ”on-top” positions. Thus, in this simple approach the step corrugation is being disregarded. We use the smallest compressions and azimuthal rotations of the Ag(111) adlayer to obtain the Moires that fit to the STM observations. For both (112)- and (335)-oriented layers the Ag layer rotation is the same 5.2, allowing the smooth transition from (112) to (335) facets observed beyond 0.52 ML. For (112) facets the Ag adlayer is compressed by 0.5% and 7.8% in the direction perpendicular and parallel to the rows, respectively. For (335) facets the perpendicular and parallel compressions are 6.2% and 7.8%, respectively. Thus, the (335)-oriented layer accumulates much more stress, explaining the preference to form (112)-oriented stripes.

The long range periodicity of the Ag-induced faceting of Fig. 1 is nicely explained within Marchenko’s elastic theory, as shown in Fig. 4. Data points represent the nanostructure period ($`L`$) and the Ag stripe width ($`w_{Ag}`$) measured from the STM images as a function of coverage ($`\theta `$) from 0.1 ML to 0.9 ML. The periodicity is poorly defined above 0.52 ML, and in this case we consider the average value $`L`$, defined as $`L=w_{Ag}+w_{Cu}`$, where $`w_{Ag}`$ and $`w_{Cu}`$ respectively are the average widths of Ag-covered facets and Cu stripes. The lines in Fig. 4 are fit to the data using Marchenko’s model, where the superlattice periodicity $`L`$ and the stripe width $`w_{Ag}`$ are given by Marchenko1 :

$$L(\theta )=\frac{\kappa }{\mathrm{sin}(\pi \theta )};w_{Ag}(\theta )=\frac{\kappa \theta }{\mathrm{sin}(\pi \theta )}$$
(1)
$`\kappa `$ is the only adjustable parameter. The model is strictly defined for a two-phase system, and hence the fit is restricted to the $`A`$-regime. The result is excellent and we obtain $`\kappa _A`$ = 91$`\pm `$6 Å. Data points at regime transitions display a clear discontinuity. Microscopically, this appears related to the way some new facets develop, i.e. by melting two contiguous facets in the previous regime. But, as indicated in Fig. 4, by changing in Eq. 1 to $`\kappa _B`$= 110$`\pm `$8 Å and $`\kappa _C`$= 141 $`\pm `$8 Å both $`w_{Ag}(\theta )`$ and $`L(\theta )`$ in $`B`$\- and $`C`$-regimes are still nicely reproduced.

In Marchenko’s approach $`\kappa `$ is related to elastic constants via $`\kappa =2\pi a\times \mathrm{exp}(1+C_1/C_2)`$, where $`a`$ is a microscopic cutoff and $`C_1`$ is the free energy of the phase boundary per unit length Marchenko1 . If we assume that work function variations between Cu and Ag phases are relatively small Vanderbilt1 , $`C_2`$ has only elastic nature and reflects the stress difference between phases Marchenko1 :

$$C_2=\frac{1\nu ^2}{\pi E}|\stackrel{}{F}|^2$$
(2)
where $`E`$ and $`\nu `$ respectively are the shear modulus and the Poisson ratio for Cu, and $`\stackrel{}{F}=(\stackrel{}{\sigma }_{Ag}\stackrel{}{\sigma }_{Cu})`$ the force exerted at facet edges due to the difference in surface stress between Ag-covered facets and Cu stripes, as defined in the inset of Fig. 4. The general agreement with Marchenko’s theory shown in Fig. 4 indicates that $`\kappa `$ in Eq. 1 does not have a significant coverage-dependence. The $`\theta `$-dependence is expected in the elastic constants due to the changing shape of Cu stripes and Ag-covered facets. Even with coverage-dependent elastic constants Eq. 1 always holds if there is a (relatively) high stress energy in any of the Ag-covered facets, such that $`C_1`$/$`C_20`$. The actual value of $`C_1/C_2`$ can be obtained from the fit in Fig. 4 assuming a cutoff length $`a`$. The cutoff length is defined as the minimum size of the (growing) domain Vanderbilt1 ; Men1 , which in this context is the terrace width. Using $`\kappa _A`$ = 91$`\pm `$2 Å and $`a`$=6.25 Å, we indeed obtain $`C_1/C_2`$= -0.15 $`\pm `$ 0.05. The absolute value of $`C_1/C_2`$ is thus much smaller than any other found in the literature, and can be assumed to be zero note3 . The condition $`C_1`$/$`C_20`$ is expected to hold for $`B`$ and $`C`$ regimes as well. The boundary energy must be very similar, because the atomic structure at facet edges is the same in all regimes, whereas the stress difference is expected to increase, since the preference for (112) orientation indicates that the lowest surface energy occurs for this facet. Thus only a change in $`a`$ would explain the changes in $`\kappa `$ in Fig. 4. The discontinuities in Fig. 4 are in fact consistent with the increasing terrace size (cutoff $`a`$) of the growing facet. The fitting parameter variation from $`\kappa _A`$ = 91 Å to $`\kappa _B`$= 110 Å (21%) and to $`\kappa _C`$= 141 Å (28%) closely correlates with the relative increase in the terrace size from $`d_{112}`$=6.25 Å to $`d_{335}`$=8.4 Å (34%) and to $`d_{223}`$=10.5 Å (25%). Since $`\kappa `$ discontinuities were microscopically related to melting of contiguous facets at regime transitions, this makes the real nature of the cutoff length an interesting issue that encourages a detailed microscopic model.

Note that the $`C_1`$/$`C_20`$ case leads to the shortest possible superlattice period in Ag/v-Cu(111). This contrasts with the long superlattice periods observed at high temperatures in NaCl/v-Cu(111) Foelsch1 . Long periods are only possible if $`C_1/C_21`$, which in the latter case could be due to a more efficient interface relaxation, favored by a large adsorbate-surface charge transfer that weakens Cu-Cu bonds.

A.R.B. and S.S. are supported by the Deutsche Forschungsgemeinschaft (DFG) and the Stichting voor Fundamenteel Onderzoek der Materie (FOM). A.M. and J.E.O. are supported by the Universidad del País Vasco (1/UPV/EHU/00057.240-EA-8078/2000)and the Max Planck Research Award Program. Fruitful discussions with A. Rubio are acknowledged.


format: text

text:
# Hard-sphere limit of soft-sphere model for granular materials: Stiffness dependence of steady granular flow

## I Introduction

The interactions between grains in flowing granular material are roughly classified into two types; the impulsive contact (collision) with the momentum exchange and the sustained contact with the transmission of forces col-fri . The flow in which the impulsive contact is dominant is called collisional flow, while the flow where the sustained contact dominates is called frictional flow. These two types of flow may be found in a simple geometry such as granular flow on a slope. The grains stay at rest when the inclination angle is too small. If the inclination exceeds a critical angle, the material starts flowing frictionally at first. The flow becomes collisional when the inclination is large enough.

As for the collisional flow of granular material, its dynamics has some analogy with molecular fluid, and the kinetic theories based on inelastic binary collisions of particles hold to some extent kinetic . On the other hand, the frictional flow is drastically different from the molecular fluid, and we have little understanding on it. Many models have been proposed for dense granular flows: For example, some models take into account the effect of non-local force transmission which comes from the network of contacting grains densemodel . In the experiments, however, it is difficult to specify the nature of sustained contact in the dense flow.

For the simulations of granular dynamics, the following two models have been commonly used, i.e. the inelastic hard sphere model and the soft sphere model.

In the inelastic hard sphere model, the particles are rigid and the collisions are thus instantaneous, therefore its dynamics can be defined through a few parameters that characterize a binary collision because there are no many-body collisions duran . The model is simple and there are very efficient algorithms to simulate it I99 , but it describes basically only the collisional flow L94  because the sustained contact is not allowed. It is also known that the system often encounters what is called the inelastic collapse K99 ; BM90 ; MY92 ; 2Dcollapse ; the infinite number of collisions take place among a small number of particles in a finite time, thus the dynamics cannot be continued beyond that point without additional assumption.

On the other hand, in the soft sphere model, that is sometimes called the discrete element method (DEM) in the granular community duran ; CS79 , the particles overlap during collision and the dynamics is defined through the forces acting on the colliding particles. Collision takes finite time, and not only binary collision but also many-body collision and sustained contact between particles are possible, therefore, both the collisional and the frictional flows may realize in the model. Many researches have been done on granular flow down a slope using the DEM in both the collisional regime and the frictional regime P93 ; MN01 ; denseDEM .

In actual simulations, however, the stiffness constants used in the soft sphere model are usually much smaller than the one appropriate for real material such as steel or glass ball YASU95  because of numerical difficulty. Therefore, some part of sustained contact in simulations may be decomposed into binary collisions if stiffer particles are used. It is also not clear that the frictional force in the sustained contact may be described by the same forces with the one suitable for the collisional events.

It is, therefore, important to examine how the system behavior may change as the stiffness constant increases in the soft sphere model, or in other words, how the soft sphere model converges to the inelastic hard sphere model in the infinite stiffness limit. In this paper, we present the results of numerical simulations of the granular flow using the soft sphere model, and investigate the system behavior when we change the stiffness constant systematically with keeping the resulting restitution constant unchanged.

In Sec. II, we briefly summarize an inelastic hard sphere model that is used for collisional granular flows. Then we introduce a simple soft sphere model for granular material and discuss how we take the hard sphere limit keeping the restitution coefficient constant. In Sec. III, the simulation results are shown. At first, the stiffness dependence of the steady state of a single particle rolling down a slope is presented to see how the inelastic collapse appears in the hard sphere limit. The collisional flow and the frictional flow are examined. We find that the interactions between particles in the collisional flow smoothly converge to binary collisions of inelastic hard spheres, while those in frictional flow shows non-trivial behavior; the behavior is also different from the “inelastic collapse” in the single particle system. The summary and the discussion are given in Sec. IV.

## II Hard sphere model and soft sphere model

In this section, we introduce the hard sphere model and the soft sphere model for granular material. For simplicity, we consider a two dimensional system and grains are modeled by two dimensional disks.

After showing the correspondence of parameters in the soft sphere model to those of the hard sphere model, we briefly summarize the phenomenon called inelastic collapse, which is the singular behavior in the inelastic hard sphere system. At the last part of this section, we discuss the limiting behavior of the soft sphere model that corresponds to the inelastic collapse of the hard sphere model.

### II.1 Inelastic hard sphere model

In hard sphere models, collisions between particles or between a particle and a wall are considered to happen instantaneously, and its dynamics is defined by the binary collision rule. We consider the collision rule in terms of the normal and tangential restitution coefficients and the sliding friction L94 ; L95 .

Let us consider a collision between the two spheres $`i`$ and $`j`$ of the diameters $`\sigma _i`$ and $`\sigma _j`$ and the masses $`m_i`$ and $`m_j`$ at the contact positions $`𝒓_i`$ and $`𝒓_j`$, respectively. Prior to the collision, the disks have velocities $`𝒄_i`$ and $`𝒄_j`$ and angular velocities $`𝝎_i`$ and $`𝝎_j`$. Then the relative velocity of the point of contact $`𝒗_{ij}`$ is given by

$$𝒗_{ij}=(𝒄_i𝒄_j)+𝒏\times \left(\frac{\sigma _i}{2}𝝎_i+\frac{\sigma _j}{2}𝝎_j\right),$$
(1)
where the normal vector $`𝒏=𝒓_{ij}/|𝒓_{ij}|=(n_x,n_y,0)`$ with $`𝒓_{ij}=𝒓_i𝒓_j`$.

If $`𝒗_{ij}^{}`$ denotes the post-collisional relative velocity, the collision rule for normal direction is

$$(𝒏𝒗_{ij}^{})=e(𝒏𝒗_{ij}),$$
(2)
where $`e`$ is the normal restitution coefficient with $`0e1`$.

In the case without sliding, the collision rule in the tangential direction is given by

$$(𝒏\times 𝒗_{ij}^{})=\beta (𝒏\times 𝒗_{ij}),$$
(3)
where $`\beta `$ is the tangential restitution coefficient with $`1\beta 1`$. The sliding is taken into account so that the tangential component of impulse does not exceed $`\mu |𝒏𝑱|`$ with the Coulomb friction coefficient $`\mu `$, where $`𝑱`$ is the momentum change of the particle $`i`$ through the collision. Namely, when the momentum change of the particle $`i`$ through the collision rules (2) and (3) is $`𝑱^{nos}`$, then $`(𝒏\times 𝒗_{ij}^{})`$ is determined so that

$$|𝒏\times 𝑱|=\text{min}(\mu |𝒏𝑱|,|𝒏\times 𝑱^{nos}|)$$
(4)
is satisfied. The detailed description of the rule is given in Refs. L94 ; L95 .

### II.2 Soft sphere model: DEM

The DEM, or the soft sphere model CS79 ; duran , is often used to simulate the dynamics of granular materials. In the present work, we adopt the two dimensional one with the linear elastic force and dissipation. When the two disks $`i`$ and $`j`$ at positions $`𝒓_i`$ and $`𝒓_j`$ with velocities $`𝒄_i`$ and $`𝒄_j`$ and angular velocities $`𝝎_i`$ and $`𝝎_j`$ are in contact, the force acting on the particle $`i`$ from the particle $`j`$ is calculated as follows: The normal velocity $`v_n`$, the tangential velocity $`v_t`$, and the tangential displacement $`u_t`$ are given by

$`v_n=𝒏𝒗_{ij},`$ $`v_t=𝒕𝒗_{ij},`$ (5)
$`u_t`$ $`=`$ $`{\displaystyle _{t_0}^t}v_t\text{d}t,`$ (6)
where $`𝒗_{ij}`$ is given in Eq. (1), and $`𝒕=(n_y,n_x,0)`$. Here, $`t_0`$ is the time when the two particles start to be in contact. Then the normal force $`F_{ij}^n`$ and the tangential force $`F_{ij}^t`$ acting on the particle $`i`$ from the particle $`j`$ are given by

$`F_{ij}^n`$ $`=`$ $`2Mk_n\left({\displaystyle \frac{\sigma _i+\sigma _j}{2}}|𝒓_{ij}|\right)2M\eta _nv_n,`$ (7)
$`F_{ij}^t`$ $`=`$ $`\mathrm{min}(|h_t|,\mu |F_n|)\text{sign}(h_t)`$ (8)
with

$$h_t=2Mk_tu_t2M\eta _tv_t,$$
(9)
where $`k_n`$ and $`k_t`$ are the elastic constants, $`\eta _n`$ and $`\eta _t`$ are the damping parameters, $`\mu `$ is the Coulomb friction coefficient for sliding, and $`M=m_im_j/(m_i+m_j)`$ is the reduced mass.

### II.3 Hard sphere limit of soft sphere model

When we adopt the linear dependence of the elastic and viscous force on the overlap and the normal relative velocity as in Eq. (7), we can calculate the duration of contact $`\tau _c`$ and the restitution coefficient for a normal binary collision $`e`$ SDW96 ; they are given by

$$\tau _c=\frac{\pi }{\sqrt{2k_n\eta _n^2}}$$
(10)
and

$$e=\mathrm{exp}\left(\frac{\pi \eta _n}{\sqrt{2k_n\eta _n^2}}\right),$$
(11)
respectively.

Neglecting the sliding friction and the variation of $`𝒏`$ during the contact, the half period of the oscillation in the tangential displacement $`u_t`$ is estimated as SDW96

$$\tau _s=\frac{\pi }{\sqrt{6k_t9\eta _t^2}}$$
(12)
for the two dimensional disks with the moment of inertia $`I_i=m_i\sigma _i^2/8`$. Here we choose the parameters $`k_t`$ and $`\eta _t`$ so that the relation $`\tau _s=\tau _c`$, or,

$$\sqrt{6k_t9\eta _t^2}=\sqrt{2k_n\eta _n^2}$$
(13)
is satisfied. Under this condition, the tangential restitution coefficient is given by

$$\beta =\mathrm{exp}(3\eta _t\tau _c).$$
(14)
Equations (10) and (11) can be rewritten as

$`\eta _n`$ $`=`$ $`\left[{\displaystyle \frac{2k_n}{(\pi /\mathrm{ln}e)^2+1}}\right]^{1/2},`$ (15)
$`\tau _c`$ $`=`$ $`\left[{\displaystyle \frac{\pi ^2+(\mathrm{ln}e)^2}{2k_n}}\right]^{1/2}.`$ (16)
Using Eqs. (13)-(16), we can take the inelastic hard sphere limit of this model for given $`e`$ and $`\beta `$ by taking the $`k_n\mathrm{}`$ limit; $`\eta _n`$, $`\eta _t`$, and $`k_t`$ are given by Eqs. (15), (14), and (13), and the duration time of collision $`\tau _c`$ goes to zero as Eq. (16).

### II.4 Inelastic collapse

It is well known that the inelastic hard sphere system can undergo the inelastic collapse, i.e., the phenomenon where infinite collisions take place within a finite period of time K99 . The simplest example is the vertical bouncing motion of a ball under gravity, but it has been shown that the inelastic collapse also occurs in higher dimensional systems without gravity K99 ; BM90 ; MY92 ; 2Dcollapse .

The inelastic collapse never occurs in the soft sphere model because of the finite length of the collision time. However, it is worth to discuss what will happen in the soft sphere model in the simple situation where the inelastic collapse occurs in the hard sphere model. First, let us consider the vertical bouncing motion of a soft ball under gravity. We assume the same force law between the ball and the floor as in Eqs. (7) and (8), except that we replace $`2M`$ by the mass of the ball. While the ball and the floor are in contact, the equation of motion for the overlap of the ball and the floor, $`z`$, is given by

$$\ddot{z}+k_nz+\eta _n\dot{z}=g,$$
(17)
with the acceleration of gravity $`g`$, where the dots represent the time derivatives. One can show that, if the impact velocity $`v_i`$ is below a critical value $`v_c`$, which is $`O(1/\sqrt{k_n})`$ for large enough $`k_n`$, the ball stays in contact with the floor; for $`v_iv_c`$, the restitution coefficient $`e_w`$ can be considered as a constant. Therefore, for a given initial impact velocity $`v_0`$, the number of necessary collisions $`n_c`$ for the ball to stay in contact with the floor is roughly estimated by the condition $`e_w^{n_c}v_0v_c`$, namely, $`n_c\mathrm{ln}\left(v_c/v_0\right)/\mathrm{ln}e_w`$. Because $`v_c1/\sqrt{k_n}`$ for large $`k_n`$, $`n_c`$ behaves as

$$n_c\mathrm{ln}k_n+\text{const.}$$
(18)
in the hard sphere limit. Thus $`n_c`$ diverges logarithmically when $`k_n\mathrm{}`$, which corresponds to the inelastic collapse due to gravity.

In the case without gravity, the inelastic collapse results in the “many-body collision” in the soft sphere model. For example, we consider three soft spheres in one dimension and assume that the binary collision can be approximated by the collision law with a constant restitution coefficient. In the situation where a particle goes back and forth between the two particles approaching each other, the inelastic collapse may take place in the hard sphere model BM90 ; three balls lose relative velocity completely in the limit of infinite collisions. In the soft sphere model, however, when the interval between two collisions becomes smaller than the duration of contact $`\tau _c`$, three balls are in contact at the same time, namely, the three body collision occurs; then they will fly apart. The number of collisions before the three body collision will also diverge logarithmically in the hard sphere limit because $`\tau _c1/\sqrt{k_n}`$.

On the other hand, one should note that a many-body collision in the soft sphere model does not necessarily result in the inelastic collapse in the hard sphere limit. Actually, in most of the cases, a many-body collision will be decomposed into a set of binary collisions in the hard sphere limit.

## III Simulation Results

In this section, we investigate the stiffness dependence of granular material on a slope in the following three situations: (i) a single particle rolling down a slope, (ii) the dilute collisional flow, and (iii) the dense frictional flow. We focus on the steady state in each situation and compare the simulation data with changing $`k_n`$ systematically. The particles are monodisperse in (ii), while they are polydisperse in (iii) in order to avoid crystallization.

In the simulations, the parameters have been chosen to give $`e=\beta =0.7`$, $`\mu =0.5`$ in the hard sphere limit. All values are non-dimensionalized by the length unit $`\sigma `$, the mass unit $`m`$, and the time unit $`\sqrt{\sigma /g}`$. Here, $`\sigma `$ is the diameter of the largest particle in the system and $`m`$ is the mass of that particle. The second order Adams-Bashforth method and the trapezoidal rule are used to integrate the equations for the velocity and the position, respectively integ . Note that the time step for integration, $`dt`$, needs to be adjusted as $`\tau _c`$ becomes smaller. All the data presented in the paper are results with $`dt=\text{min}(\tau _c/100,10^4)`$. We have confirmed that the results do not change for $`dt\tau _c/100`$ by calculating also with $`dt=\tau _c/50`$ and $`dt=\tau _c/200`$ in the case of the single particle.

### III.1 A single particle rolling down a bumpy slope

Let us first consider a single particle rolling down a bumpy slope. It is known that the particle shows the steady motion for a certain range of the inclination angle $`\theta `$ DBW96 . In the simulations, we make the boundary rough by attaching the same particles with the rolling one to the slope with the spacing $`0.002\sigma `$ (see Fig. 1). For the chosen parameters with the normal stiffness $`k_n=2^1\times 10^5`$, the range of $`\theta `$ for which steady motion is realized is $`0.11\mathrm{sin}\theta 0.14`$ MN01 . Here we fix the inclination angle to $`\mathrm{sin}\theta =0.13`$. Figure 2 shows the time evolution of the velocity in the $`y`$ direction, $`v_y`$, with $`k_n=2^1\times 10^5`$ (solid line) and $`k_n=2^{17}\times 10^5`$ (dashed line). It is shown that $`v_y`$ behaves periodically; this period ($`\mathrm{\Delta }t4`$) corresponds to the period for the particle to get past one particle at the floor. The period hardly depends on the stiffness.

In Fig. 3 (a), the $`k_n`$ dependence of the time averaged kinetic energy $`E`$ is shown. It is shown that the energy is an increasing function of $`1/k_n`$ in the softer region ($`1/k_n10^7`$), but no systematic $`1/k_n`$ dependence of $`E`$ exists for $`1/k_n10^7`$. The average collision rate (number of collisions per unit time) between the slope and the particle , $`N_w`$, shows logarithmic dependence on $`1/k_n`$ in the whole region (Fig. 3 (b)). From Fig. 3 (c), we also find that the average contact time fraction with the slope, $`t_w`$, is a decreasing function of $`k_n`$ in the soft region, but it seems to approach a constant value for large enough $`k_n`$.

The logarithmic $`k_n`$ dependences of $`N_w`$ and the constant $`t_w`$ in large $`k_n`$ region agree with our previous analysis of “inelastic collapse under gravity” in the soft sphere model in Section II.4. The motion of the particle in one period is as follows; when the particle comes to a bump (a particle attached to the slope), the particle jumps up, bounces on the bump many times, loses the relative velocity, and finally rolls down keeping in contact with the bump. Therefore, the contact time fraction has a finite value even in the hard sphere limit due to the rolling motion at the last part. $`N_w`$ increases logarithmically in the hard sphere limit as has been expected from Eq. (18).

### III.2 Collisional flow

Next we consider the steady state of the collisional flow. The system considered is shown in Fig. 4. The particles are monodisperse, and the slope is made rough as in the single particle case. The periodic boundary condition is adopted in the flow direction ($`x`$ direction). The length of the slope is $`L=50.1`$ and the number of the particles attached to the slope is $`50`$. The number of flowing particles is also $`50`$, namely, the number of the particles per unit length along the slope is about $`1`$. The inclination angle is set to be $`\mathrm{sin}\theta =0.45`$. The initial configuration of particles is the row of 50 particles at rest with regular spacing in the $`x`$-direction, but each particle is at random height in the $`y`$-direction. After a short initial transient, if the total kinetic energy fluctuates around a certain value, we consider it as the steady flow. All the data were taken in this regime and averaged over the time period of $`1500`$.

As can be seen in the snapshot, Fig. 4, the particles bounce and the number of particles in contact is very small. In Fig. 5 (a), the averaged kinetic energy per one particle, $`E`$, is shown. Though $`E`$ becomes larger as the particles become softer in $`1/k_n10^5`$, the systematic dependence of $`E`$ on $`k_n`$ disappears for large enough $`k_n`$ ($`1/k_n10^6`$) commentyuu . The $`y`$-dependence of the average number density and the flow speed in this region are also shown in Fig. 6 (a) and (b).

In Fig. 5 (b), the average collision rates between particles, $`N_c`$ (filled circles), and between particles and the slope, $`N_w`$ (open circles), per particle are plotted vs $`1/k_n`$. Both of them stay roughly constant in the region where $`E`$ is almost constant, $`1/k_n10^6`$. On the other hand, the average time fractions during which a particle is in contact with other particles, $`t_c`$ (filled circles), and in contact with the slope, $`t_w`$ (open circles), decrease systematically as $`k_n`$ increases as shown in Fig. 5 (c). Both the solid and the dashed lines are proportional to $`1/\sqrt{k_n}`$, namely, $`t_c`$ decreases in the same manner as $`\tau _c`$ in Eq. (16). Actually, the collision time fractions $`t_c`$ and $`t_w`$ converge to $`N_c\tau _c`$ and $`N_w\tau _w`$ ($`\tau _w`$ is the duration of a normal collision of a particle and the floor comment ), respectively, in the hard sphere limit. The differences $`t_cN_c\tau _c`$ and $`t_wN_w\tau _w`$ are plotted in Fig. 5 (d) to show that they go to zero very rapidly upon increasing $`k_n`$. This means that the interactions of the soft sphere model in the collisional flow regime converge to those of the inelastic hard sphere model with binary collisions.

If we look carefully, however, in the large $`k_n`$ region in Fig. 5 (c), we can see slight deviation of $`t_w`$ from the dashed lines; it decreases slower than $`1/\sqrt{k_n}`$. This tendency may indicate that $`t_w`$ remains finite in the hard sphere limit: Actually, in the event-driven simulation of the hard sphere model, we always found the inelastic collapse as long as the restitution constant between a particle and the floor is less than $`1`$. This suggests that $`N_w`$ should diverge and $`t_w`$ should remain finite in the hard sphere limit because of the inelastic collapse due to gravity. The slight deviation of $`t_w`$ from the dashed line in Fig. 5 (c) may be a symptom of it, while we cannot see the logarithmic divergence in $`N_w`$.

### III.3 Frictional flow

Now we study the steady state of the dense frictional flow. We adopted the flat boundary with slope length $`L=10.02`$ and imposed the periodic boundary condition in the flow direction. In order to avoid crystallization, we used the polydisperse particles with the uniform distribution of diameter from $`0.8`$ to $`1.0`$. The number of particles in the system is $`100`$. The inclination angle $`\theta `$ is set to be $`\mathrm{sin}\theta =0.20`$ and the initial condition is given in the similar way with the collisional flow. The ten rows of ten particles with regular spacing in the $`x`$-direction are at rest with large enough spacings between rows in the $`y`$-direction so that particles do not overlap; only the particles in the top row are at random heights in $`y`$-direction.

It turns out that the steady state is not unique and fluctuation is large. Depending on the initial conditions, the particles may flow with a different value of average kinetic energy, or in some cases, the whole system stops. We presume that this is mainly because the system is not large enough: Grains form a layer structure as in Fig. 7, and the flow velocity strongly depends on the configuration of particles in the bottom layer. The time evolutions of the kinetic energy per particle, $`E(t)`$, of three samples with different configurations in the bottom layer are shown in Fig 8. Sample 1 shows stable behavior while sample 2 eventually stops after running with lower energy for some time. Sample 3 shows larger fluctuation with higher energy; the sudden change of $`E(t)`$ in sample 3 results from the change in the configuration of the bottom layer. This large fluctuation should be averaged out if we could simulate large enough system for long enough time, but the amount of computation is too large especially for the system with stiff particles.

In order to make meaningful comparison out of these largely fluctuating data with a variety of behaviors, we select simulation sequences that come from similar flowing behaviors in the following way. First, we define the steady part of the time sequence in each of the samples as the part where the width of the energy fluctuations is smaller than 0.45 over the time period longer than 500. Second, we exclude the data whose averaged energy is out of the range \[0.4,0.8\] select .

Many of the excluded data by this criterion show quite different flowing behaviors. We use only the data selected from this criterion to calculate time averages of physical quantities.

In the snapshot of the frictional flow, Fig. 7, most of the particles seem to be in contact with each other and form a layer structure. The $`y`$-dependence of the average density and the flow speed are shown in Figs. 9 (a) and (b), respectively. We can see that the relative motion between layers is largest at the bottom and very small in the bulk. The stiffness dependence of the average kinetic energy $`E`$ is shown in Fig. 10 (a); the data are scattered due to the non-uniqueness of the steady state.

Nevertheless, the average collision rates show systematic dependence as shown in Fig. 10 (b). Here, the definition of $`N_c`$ ($`t_c`$) is the same as that for the collisional flow, namely, it is the average collision rate (contact time fraction) between particles per particle in the system. On the other hand, $`N_w`$ ($`t_w`$) is defined differently; it is the collision rate (the contact time fraction) between the particles and the slope per particle in the bottom layer, because other particles never touch the slope.

In Fig. 10 (b), we can see that $`N_c`$ (filled circles) and $`N_w`$ (open circles) increase very rapidly as $`k_n`$ becomes larger: they diverge as a power of $`k_n`$. This is quite different from the behavior in the collisional flow in which they are almost constant. Furthermore, this increase is faster than the logarithmic divergence found in the single particle case (see Fig. 3 (b)).

The contact time fractions, $`t_c`$ (filled circles) and $`t_w`$ (open circles), decrease as shown in Fig. 10 (c). The main reason why $`t_c`$ and $`t_w`$ decrease is that $`\tau _c`$ and $`\tau _w`$ decrease faster than $`N_c`$ and $`N_w`$, namely, the contact time fractions estimated from the duration of a binary collision, $`N_c\tau _c`$ and $`N_w\tau _w`$, continue to decrease. Actually, the decrease of the contact time fraction and the increase of the collision rate are natural because a longer multiple contact breaks up into shorter binary contacts as the particles become stiffer.

These contact time fractions, $`t_c`$ and $`t_w`$, however, do not converge to $`N_c\tau _c`$ and $`N_w\tau _w`$, respectively. As shown in Fig. 10 (d), $`t_cN_c\tau _c`$ (filled circles) and $`t_wN_w\tau _w`$ (open circles) are very large as compared to those in the collisional flow, even in the stiffest region. The comparison of this with the rapid convergence in the collisional flow regime (Fig. 5 (d)) indicates that there remains finite multiple contact time in the hard sphere limit and the interaction in the frictional flow can never be considered as many, or even infinite, instantaneous binary collisions. The particles experience the lasting multiple contact even in the hard sphere limit.

## IV Summary and Discussion

The inelastic hard sphere limit of granular flow has been investigated numerically in the steady states of (i) a single particle rolling down the slope, (ii) the dilute collisional flow, and (iii) the dense frictional flow. In (i), it has been found that the “inelastic collapse” between the particle and the slope occurs in the hard sphere limit due to gravity, and the contact time fraction between the particle and the slope remains finite. In (ii), the collision rates $`N_c`$ and $`N_w`$ are almost constant when particles are stiff enough. The contact time fraction between particles $`t_c`$ approaches zero as $`k_n`$ increase in the same manner with the duration of contact for binary collision, $`\tau _c`$. This means that the interaction between particles in the hard sphere limit can be expressed by binary collisions in the inelastic hard sphere model. On the other hand, the decrease in $`t_w`$ is slightly slower than $`1/\sqrt{k_n}`$ in the harder region, which can be a sign of the inelastic collapse between a particle and the slope.

In the case of the frictional flow (iii), the situation is not simple: Although the contact time fractions $`t_c`$ and $`t_w`$ decrease upon increasing $`k_n`$, the collision rates $`N_c`$ and $`N_w`$ increase as a power of $`k_n`$, which is faster than the logarithmic divergence found in the single particle case (i). The origin of this power divergence of collision rate does not seem to be simple because it is a property of the steady state, not a particular dynamical trajectory of the system.

The multiple contact time fraction may be estimated by $`t_cN_c\tau _c`$ and $`t_wN_w\tau _w`$. They were found to be quite large compared to those for the collisional flow, and seems to remain finite even in the infinite stiffness limit: This suggests that the interaction in the frictional flow can never be considered as infinite number of binary collisions. Even in the hard sphere limit, particles experience the lasting multiple contact.

The non-negligible fraction of multiple contact in the hard sphere limit implies the existence of the network of contacting grains even though they are flowing. The models of dense flows should consider the effect of this lasting contacts.

Here we have investigated only the two cases, i.e., the collisional flow and the frictional flow. The system should undergo the transition between the two flows if we change continuously the parameters such as the inclination angle, the roughness of the slope, or the density of particles. It is interesting to investigate how the transition occurs by looking at the quantities measured in this paper, because their stiffness dependences are qualitatively different in the two flows.

###### Acknowledgements.

This research was partially supported by Hosokawa powder technology foundation, the Japan Society for the Promotion of Science, and Grant-in-Aid from the Ministry of Education, Science, Sports and Culture of Japan.


text:
# Jet Physics at Two-Loop Accuracy Presented at RADCOR2002 / Loops and Legs in Quantum Field Theory, September 2002, Kloster Banz, Germany

## 1 Introduction

Jet production observables are among the most sensitive probes of QCD at high energy colliders, where they are used for example to determine the strong coupling constant. At present, the interpretation of jet production data within perturbative QCD is restricted to next-to-leading order (NLO) calculations, with theoretical uncertainties considerably larger than current experimental errors. Going beyond NLO calculations offers moreover a more accurate matching of theoretical and experimental jet definitions and a more detailed modelling of the hadronic final state . The extension of jet calculations to NNLO requires three ingredients: the two-loop corrections to multi-leg amplitudes, the single unresolved limits of one-loop amplitudes and the double unresolved limits of tree amplitudes. Finally, all these contributions have to be combined together into a program for the numerical computation of jet observables from NNLO parton level cross sections. In this talk, I review recent progress made on these subjects as well as the currently open issues.

## 2 Virtual Two-Loop Corrections

Within dimensional regularization, the large number of different integrals appearing in multi-loop calculations can be reduced to a small number of so-called master integrals by using integration-by-parts (IBP) identities . These identities exploit the fact that the integral over the total derivative of any of the loop momenta vanishes in dimensional regularization.

For integrals involving more than two external legs, another class of identities exists due to Lorentz invariance. These Lorentz invariance identities  rely on the fact that an infinitesimal Lorentz transformation commutes with the loop integrations, thus relating different integrals. The common origin of IBP and LI identities is the Poincare invariance of loop integrals within dimensional regularization, as was pointed out by J.J. van der Bij at this conference. Using integration-by-parts and Lorentz invariance identities, all two-loop Feynman amplitudes for $`22`$ scattering or $`13`$ decay processes can be expressed as linear combinations of a small number of master integrals, which have to be computed by some different method. Explicit reduction formulae for on-shell two-loop four-point integrals were derived in . Computer algorithms for the automatic reduction of all two-loop four-point integrals were described in .

The master integrals relevant to $`22`$ scattering or $`13`$ decay processes are massless, scalar two-loop four-point functions with all legs on-shell or a single leg off-shell. Several techniques for the computation of those functions have been proposed in the literature, such as the application of a Mellin-Barnes transformation to all propagators, which was used successfully to compute the on-shell planar double box integral  , the on-shell non-planar double box integral  and two double box integrals with one leg off-shell . Most recently, the same method was used to derive the on-shell planar double box integral with one internal mass scale  as well as the high energy limit of the on-shell planar triple box integral .

A method for the analytic computation of master integrals avoiding the explicit integration over the loop momenta is to derive differential equations in internal propagator masses or in external momenta for the master integral, and to solve these with appropriate boundary conditions. The computation of master integrals from differential equations proceeds as follows . Carrying out the derivative with respect to an external invariant on the master integral of a given topology, one obtains a linear combination of a number of more complicated integrals, which can however be reduced to the master integral itself plus simpler integrals by applying the reduction methods discussed above. As a result, one obtains an inhomogeneous linear first order differential equation in each invariant for the master integral. The inhomogeneous term in these differential equations contains only topologies simpler than the topology under consideration, which are considered to be known if working in a bottom-up approach. The master integral is then obtained by matching the general solution of its differential equation to an appropriate boundary condition.

Using the differential equation technique, one of the on-shell planar double box integrals  as well as the full set of planar and non-planar off-shell double box integrals  were derived. The computer algebra structures applied in the computation of the master integrals from the differential equations are displayed in the right hand column of Figure 1. The differential equation approach has recently been extended to phase space integrals .

A strong check on all these computations of master integrals is given by the completely numerical calculations of , which are based on an iterated sector decomposition to isolate the infrared pole structure. The methods of  were applied to confirm all of the above-mentioned calculations.

A third approach, which avoids the reduction to master integrals, has been presented in . In this approach, all integrals appearing in the two-loop amplitudes are related to higher transcendental functions, which can be expanded in terms of nested harmonic sums.

The two-loop four-point functions with all legs on-shell can be expressed in terms of Nielsen’s polylogarithms . In contrast, the closed analytic expressions for two-loop four-point functions with one leg off-shell contain two new classes of functions: harmonic polylogarithms  and two-dimensional harmonic polylogarithms (2dHPL’s) . Accurate numerical implementations for these functions  are available.

### 2.1 $`22`$ Processes with all legs on-shell

With the explicit solutions of the integration-by-parts and Lorentz-invariance identities for on-shell two-loop four-point functions and the corresponding master integrals, all necessary ingredients for the computation of two-loop corrections to $`22`$ processes with all legs on-shell are now available. The generic structure of such a calculation is outlined in the left hand column of Figure 1. In fact, only half a year elapsed between the completion of the full set of master integrals and the calculation of the two-loop QED corrections to Bhabha-scattering . Subsequently, results were obtained for the two-loop QCD corrections to all parton-parton scattering processes . For gluon-gluon scattering, the two-loop helicity amplitudes have also been derived . Moreover, two-loop corrections were derived to processes involving two partons and two real photons . Since the gluon fusion into photons has a vanishing tree level amplitude, these results form part of the NLO corrections to photon pair production , yielding a sizable correction.

Finally, light-by-light scattering in two-loop QED and QCD was considered in .

The results for the two-loop QED matrix element for Bhabha scattering  were used in  to extract the single logarithmic contributions to the Bhabha scattering cross section.

### 2.2 $`22`$ Processes with one off-shell leg

Using the two-loop master integrals with one off-shell leg , the two-loop QCD matrix element for $`e^+e^{}3`$ jets  and the corresponding helicity amplitudes  were computed following the reduction procedure depicted in Figure 1. The infrared pole structure of these results agrees with the prediction  obtained from an infrared factorization formula.

An independent confirmation of part of these results was performed in , where two of the seven colour factors (corresponding to the terms proportional to $`n_f`$) of the two-loop helicity amplitudes for $`e^+e^{}3`$ jets were derived using the nested sum method of .

Processes related to $`e^+e^{}3`$ jets by crossing symmetry are $`(2+1)`$-jet production in deep inelastic $`ep`$ scattering and vector-boson-plus-jet production at hadron colliders. The analytic continuation of the $`e^+e^{}3`$ jets two-loop helicity amplitudes to the kinematic regions relevant for these scattering processes has been derived in .

## 3 Real Corrections

Besides the two-loop virtual corrections, a NNLO calculation of jet observables has to include the contributions from single unresolved (soft or collinear) real radiation from one-loop processes as well as from double unresolved real radiation at tree level. Only after summing all these contributions (and including terms from the renormalization of parton distributions for processes with partons in the initial state), do the divergent terms cancel among one another. The factorization properties of both the one-loop, one-unresolved-parton contribution  and the tree-level, two-unresolved-parton contributions  have been studied, but a systematic procedure for isolating the infrared singularities has so far been established only for the one-loop, one-unresolved-parton processes. Although this is still an open and highly non-trivial issue, significant progress is anticipated in the near future.

## 4 Numerical Implementation

Prior to their implementation into a numerical program, it is needed to analytically extract the infrared pole terms from the one-loop, one-unresolved-parton and two-unresolved-parton contributions. At NLO, two types of methods have been used very successfully in the past. The phase space slicing method  divides up the final state phase space into resolved and unresolved regions; the infrared subtraction terms are then integrated only over the unresolved regions. This method avoids overcounting of singular contributions, the required integrals over restricted phase space regions can however be very involved beyond NLO. The subtraction method  integrates the subtraction terms over the full phase space. The construction of the subtraction term requires in this case great care to avoid overcounting problems. In general, algorithms based on subtraction are more efficient numerically.

The remaining finite terms must then be combined into a numerical program implementing the experimental definition of jet observables and event-shape variables. The sketch of such a programme to compute $`e^+e^{}3j`$ at NNLO is given in Figure 2.

Programs to compute processes with initial state hadrons involve the additional complication of initial state singularities, which have to be absorbed into the NNLO parton distributions. Lacking the full expressions for the splitting functions at this order, these are not yet available at present, work on them is however well advanced .

A first calculation involving the features of NNLO jet calculations was presented for the case of photon-plus-one-jet final states in electron–positron annihilation in , thus demonstrating the feasibility of this type of calculations. A prerequisite for such a numerical program computing $`n`$ jet final states is a stable and efficient next-to-leading order programme for the processes yielding $`n+1`$ jet final states. For the processes of highest phenomenological interest, these are already available: $`e^+e^{}4j`$ , $`ep(3+1)j`$ , $`pp3j`$ , $`ppV+2j`$ .

## 5 Conclusions and Outlook

Considerable progress has been made in the last two years (in fact since the last “Loops and Legs”- and “RADCOR”-conferences) towards the computation of jet observables at NNLO in QCD. In particular, new methods have been developed for the calculation of two-loop virtual corrections to four-point scattering amplitudes. As a result, the two-loop virtual corrections relevant to all phenomenologically important processes in QCD and QED are now known.

These form however only part of the full calculation required at NNLO accuracy, which also has to take into account contributions from one-loop single-unresolved radiation and tree-level double unresolved radiation processes. While appropriate subtraction terms for these processes have been known for quite some time, their analytic integration (required for the cancellation of infrared poles in physical jet observables) is still an unsolved problem. Once this obstacle has been overcome, the remaining finite parts can be implemented into a numerical programme to compute NNLO jet observables.

## Acknowledgement

I wish to thank Ettore Remiddi, Nigel Glover, Lee Garland and Thanos Koukoutsakis for a pleasant and fruitful collaboration on the topics discussed in this talk.


text:
# A Catalogue of infrared star clusters and stellar groups

## 1 Introduction

The development of infrared detectors in the last decades has unveiled many optically absorbed star clusters which are embedded in dark clouds or in more evolved star forming regions (Hodapp 1994). These clusters which are fully or partially immersed in interestellar gas and dust, are called embedded clusters (Lada & Lada 1991). More generally infrared star clusters can be employed, since some of them can be older and not being embedded in complexes. Embedded clusters are ideal targets to investigate early stages of star formation and evolution. Many studies were made with newly identified clusters (e.g. Deharveng et al. 1999, Carpenter 2000 and Horner et al. 1997). It is important to gather all objects in a catalogue, which in turn will facilitate future developments, and avoid problems such as duplications. It allows the identification of less explored and unexplored areas that could be surveyed. In particular, near infrared imaging surveys such as the Two Micron All Sky Survey (hereafter 2MASS, Skrutskie et al. 1997) and the Deep NIR Southern Sky Survey (DENIS; Epchtein et al. 1997) are providing material for cluster identifications (Soares & Bica 2002, Reylé & Robin 2002).

Dutra & Bica (2000a) identified 58 new rather faint infrared star cluster candidates projected on the central parts of the Galaxy, using the 2MASS survey. The images employed by 2MASS have a relatively low angular resolution and the fields near the Galactic Centre are so crowded that is not possible to construct CMDs and clearly establish the nature of these objects with the available material. Thus they are candidates to infrared clusters to be verified with larger telescopes under good seeing, and we preferred not include them here. In a forthcoming study we will present an infrared survey carried out at a higher angular resolution using the ESO NTT. More recently, Dutra & Bica (2001) have made a search of objects in the 2MASS atlas in different parts of the sky finding 42 star clusters and stellar groups, which are included in the present study.

Mermilliod (1996) compiled data on open clusters which are updated in the Web interface http://obswww.unige.ch/webda. Recently, Mermilliod created a section for embedded clusters in his star cluster database and included a few entries. He points out the need to compile data for these objects and calls for contributions. One of the objectives of the present study is thus to compile a catalogue of infrared star clusters.

Star cluster pairs and multiplets have been studied in detail in the LMC (Dieball et al. 2002), and SMC (de Oliveira et al. 2000). Multiplicity appears to be an important phenomenon which affects the dynamical evolution of young star clusters, e.g. in the SMC (de Oliveira et al. 2000). Some open cluster pairs appear to occur in the Galaxy (Subramaniam et al. 1995). However not much is known about the rate of cluster multiplicity at their birth, and infrared imaging of embedded clusters is now opening the possibility for such studies (e.g. the cluster pair NGC1333, Lada et al. 1996). Another objective of the present work is to verify cluster multiplicity in the infrared cluster catalogue.

In Sect. 2 we present the catalogue and its references. In Sect. 3 we discuss the properties of the sample such as angular and size distribution and cluster multiplicity. Finally, in Sect. 4 concluding remarks are given.

## 2 The catalogue

Infrared detectors and surveys have boosted the studies of infrared clusters, most of them embedded in dark or bright nebulae. It is important to gather these objects into a catalogue. During the search we realized that several objects were not as dense as star clusters, but appear to be physical systems. We illustrate the difference between the two object classes by means of a comparison of richness parameter for the DR22 Cluster and its surrounding DR22 stellar group (Fig. 1 of Dutra & Bica 2001). The DR22 cluster is rather loose and uncrowded, an example of a less dense cluster. For the cluster we obtained 54 stars (arcmin)-2, while for the stellar group 25 stars (arcmin)-2, as measured on the 2MASS atlas. No background or crowding correction was applied. Basically our classification in infrared star cluster or stellar group follows a density estimation by eye on the available material. Some authors also adopted this classification in the object’s original reference. Other objects are concentrations of T Tauri stars (e.g. those in Taurus, Gomez et al. 1993 and Leinert et al. 1993, and Chamaleon, Whittet et al. 1991), also classified as groups. Consequently, the catalogue is formed by two object classes, the standard infrared star clusters and the less dense stellar groups. The present study includes 189 infrared star clusters and 87 infrared stellar groups from 86 references.

The nomenclature of infrared star clusters and groups is not uniform. Some authors adopted IRAS or AFGL related sources, dark nebulae, optical or radio HII regions, running numbers, among others. We included as many as possible designations. We preferred as first designation a related optical HII region, radio source, dark nebula and infrared sources, if multiple possibilities occur.

The catalogue has 276 objects and will be available as Table 1 in electronic form at CDS (Strasbourg) via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5). Each entry occupies 2 catalogue lines. The first line shows by columns: (1) and (2) Galactic coordinates, (3) and (4) equatorial coordinates (J2000.0 epoch), (5) and (6) the major and minor axes dimensions in arcmin and (7) object designations. On line 2 are included a class discriminator (IRC - infrared cluster or IRGr - infrared stellar group), comments and additional information such as distance in kpc (d=), age in Myr (t=), reddening (eFIR=) – see below, linear diameter in pc (LD=) and a reference code (REFnnn). The references of the catalogue will be available as Table 2 in electronic form at CDS (Strasbourg) via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5). It includes in Col. 1 the reference code from Table 1 and in Col. 2 the complete reference. Distances and ages were taken from the individual studies, when available. Angular diameters were measured on the images in the individual studies and on the 2MASS Atlas, when available. Reddening estimates based on CMDs are not systematically available. Consequently, we now prefer to provide the reddening estimates based on the object’s coordinates extracted from the Schlegel et al.’s (1998) all-sky reddening maps, which are based on the 100$`\mu `$m dust emission. For each cluster or stellar group we include the respective E(B-V)FIR reddening value (”eFIR=” in Table 1) associated to its line of sight. We point out that this estimate is an upper limit for the reddening, since it represents the total reddening in the object’s line of sight, and background contributions may occur (Dutra & Bica 2000b). We also point out that it can be overestimated by heated dust associated to star forming regions (Dutra et al. 2002 and references therein).

## 3 Sample Properties

Overall properties of the catalogue can be inferred from the available data. We will analyse the angular distribution, distance and size distributions and multiplicity. Since the vast majority of the objects are embedded, the conclusions basically refer to the initial stages of star cluster and stellar group formation.

### 3.1 Angular distribution

Fig. 1 shows the angular distribution of the objects in galactic coordinates. As expected infrared star clusters are projected close to the Galactic Plane. We note that the Anticentre zone and Cygnus X ($`\mathrm{}=80^{}`$) are well explored. On the other hand, the southern Milky Way ($`240^{}<\mathrm{}<340^{}`$) requires a deep survey. The stellar groups tend to occur away from the Plane since they are mostly related to the nearby complexes (e.g. Orion at $`\mathrm{}`$ 210 and $`b`$ -20, Taurus at $`\mathrm{}`$ 175 and $`b`$ -15, and Chamaleon I at $`\mathrm{}`$ 297 and $`b`$ -15).

### 3.2 Distance and Size distribution

Fig. 2 shows the distance histograms for the infrared cluster and stellar group samples. Most objects in the cluster sample are located up to distances from the Sun of 3 kpc. The present sample is thus biased to relatively nearby zones. This shows the need to explore more distant disk regions. The histogram of stellar groups shows that they are mostly as far as 1 kpc, since they are associated to nearby complexes (e.g. Taurus and Orion).

We have used the distances and angular sizes provided in the catalogue to derive linear diameters for the clusters and stellar groups. Fig. 3 shows the linear diameter histograms for the infrared clusters and stellar groups. The size distribution of clusters shows that the typical one is small. The peak diameter is $``$ 1 pc and the majority has linear diameter smaller than 3 pc. Open clusters older than 15 Myr are well distributed between 1 and 8 pc (Janes et al. 1988). Embedded clusters tend to be smaller than evolved open clusters. This effect appears to be physical since a detection bias among the somewhat older open clusters is not expected for in this case the most prominent IR clusters in a earlier stage would produce a series of prominent compact open cluster population and this is not observed. Nor a star formation bias is expected since the Galaxy has no apparent reason to have currently formed prominent clusters in excess. We are left with an age effect, related to the dynamical evolution of the disk clusters in general. Possibly we are witnessing an evolutionary effect from star cluster formation to those reflecting the internal dynamical evolution and the tidal effects of the Galaxy.

### 3.3 Multiplicity

During the present catalogue construction we realized that often embedded clusters have companions a few diameters away. These systems are indicated as pairs (mP) and triplets (mT) in Table 1. We found 18 pairs and 4 triplets resulting 48 clusters as components. Therefore, 25 % the embedded cluster sample is formed in multiple systems. Not many open cluster pairs are known in the Galaxy (Subramaniam et al. 1995), consequently the dynamical evolution of double and triple embedded clusters probably leads to a high rate of mergers.

We also indicated groups of clusters and stellar groups which can have as many as 9 components (m9 in Table 1). Such ensembles can be related to the same molecular cloud such as the clusters in the Rosette Molecular Cloud (RMC), which has 8 cluster members. The separations in these groups tend to be larger then those in the pair and triplet systems.

The possible link between duplicity effects and age and/or extinction selection effects must await a larger number of objects with CMD’s and homogeneous age and extinction determinations. A comparison sample of non-embedded clusters, with ages of a few 107 yr is also required.

## 4 Conclusions

We built a catalogue of infrared star clusters and stellar groups (less dense objects) in the Galaxy, using 86 references from the literature. A total of 276 entries are included, and we provide coordinates, angular dimensions, designations and additional information when available, such as distances and linear diameters. We analysed overall properties of the sample. The angular distribution of the sample showed a deficiency of objects in the southern Milky Way indicating the need of a systematic search. The distance distribution indicated that the sample is biased to nearby zones ($`d<3`$ kpc). The Milky Way disk beyond this limit also requires more exploration. The linear size distribution revealed that embedded clusters are small, typically 1 pc. We noted that 25 % of the embedded clusters occur in pairs or triplets with separation of a few diameters. Such proximity must affect their dynamical evolution, possibly resulting in mergers. For the first time a catalogue of this nature is presented and can be useful for future searches and developments.

###### Acknowledgements.

We acknowledge support from the Brazilian Institutions CNPq and FAPESP. CMD acknowledges FAPESP for a post-doc fellowship (proc. 00/11864-6).


text:
# 1 fig:tcube

Universal Power-law Decay in Hamiltonian Systems?

The understanding of the asymptotic decay of correlations and of the distribution of Poincaré recurrence times $`P(t)`$ has been a major challenge in the field of Hamiltonian chaos for more than two decades. In a recent Letter, Chirikov and Shepelyansky  claimed the universal decay $`P(t)t^3`$ for Hamiltonian systems. Their reasoning is based on renormalization arguments and numerical findings for the sticking of chaotic trajectories near a critical golden torus in the standard map. We performed extensive numerics and find clear deviations from the predicted asymptotic exponent of the decay of $`P(t)`$. We thereby demonstrate that even in the supposedly simple case, when a critical golden torus is present, the fundamental question of asymptotic statistics in Hamiltonian systems remains unsolved.

As in Ref.  we study the standard map

$$q_{n+1}=q_n+p_n\mathrm{mod}\mathrm{\hspace{0.17em}2}\pi p_{n+1}=p_n+K\mathrm{sin}q_{n+1},$$
(1)
at $`K=K_c=0.97163540631`$, where the golden torus is critical (Fig. 1, inset). We determine the Poincaré recurrence time distributions $`P(t)`$ for trajectories starting below and above the critical golden torus by using the same numerical approach as in Ref. . By considerably increasing the statistics we are able to extend the distribution by almost two orders of magnitude in recurrence times. We verify that our statistical data are not affected by the unavoidable finite numerical precision by comparing data for double ($``$16 significant digits) and quadruple ($``$32 digits) precision. The data for approaching the critical golden torus from above and below are presented in Fig. 1. For times $`t<10^8`$ our data agree with the results presented in Fig. 2 of Ref. . For larger times, however, we find strong deviations from the predicted universal power law $`P(t)t^3`$ (dashed lines in Fig. 1). The deviations might be explained in two ways: The onset of the claimed asymptotic decay might occur for larger times, which is in contradiction to the prefactors determined in Ref. . On the other hand, the long-time trapping of chaotic trajectories might be dominated by islands of stability (non-principal resonances) that are neglected by the renormalization arguments. In fact, the latter possibility is supported by a detailed investigation .

If even in the supposedly simple case of a critical golden torus the decay $`P(t)t^3`$ is not observed, the claim for a universal existence of this decay cannot be maintained. It thus remains a fundamental challenge in the field of Hamiltonian chaos whether the asymptotic behavior of $`P(t)`$ follows a universal power law and what the value of its exponent would be.

We thank Dima Shepelyansky for candid discussions.

M. Weiss1,2, L. Hufnagel1, and R. Ketzmerick1

1 Max-Planck-Institut für Strömungsforschung and Institut für Nichtlineare Dynamik der Universität Göttingen, Germany

2 EMBL, Meyerhofstr. 1, 69117 Heidelberg, Germany

PACS numbers: 05.45.Mt


text:
# 𝑝-adic description of characteristic relaxation in complex systems

## 1 Introduction.

One of the distinctive features of complex systems (such as glasses, macromolecules, and proteins) is an anomalously slow non-exponential relaxation observed in a very wide range of time scales — from that of short-time molecular vibrations to days and months. Relaxation kinetics in complex systems is often fitted into three empirical laws, namely, the stretched exponential curve (the Kohlrausch-Williams-Watts law), $`\mathrm{exp}\left[\left(t/\tau \right)^\alpha \right]`$, $`0<\alpha <1`$, the power decay law, $`\left(t/\tau \right)^\alpha `$, $`0<\alpha `$, and the logarithmic decay law, $`\alpha \left[ln\left(t/\tau \right)\right]^1`$, $`1<\alpha `$ . These types of relaxation are characteristic of complex systems, just as exponential relaxation is typical for gases and liquids.

Although formulation of a microscopic theory of anomalous relaxation in complex systems remains an issue of the day, there exist important concepts regarding the general relaxation behavior of complex systems. The starting point of these concepts is that the states of a micro-component of a system are strongly coupled to configurational rearrangements of the local environment, and it is important that these rearrangements cover a wide range of scales – from microscopic to mesoscopic. In order to describe configurational dynamics of this kind, the concept of energy landscapes is applied -. According to this concept, configurational rearrangements are described by a random walk in the space of configurational states, and this process, in its turn, is determined by the energy landscape of the system.

There is common understanding of the fact that a comprehensive description of energy landscapes of complex systems is hardly possible for several reasons. First of all, due to an uncertainty while selecting local potentials of interaction between the system’s components. Another reason is the large number of metastable configurations, which occur due to multiple quenched constraints and the random nature of the local arrangement. The number of metastable configurations grows exponentially with respect to the size of the system. Thus, energy landscapes of mesoscopic systems (these systems are of the main interest) form multidimensional surfaces with an extremely large number of local minima. It is not quite clear how to describe those surfaces analytically and, accordingly, how to describe the configurational dynamics constrained by such landscapes. In this connection, computer simulation of model systems remains the main field of activity. These studies show that energy landscapes of complex systems may vary significantly. However, the characteristic types of relaxation obviously suggest that such landscapes have some common features. These features might turn out to be as fundamental for the dynamics of complex systems as the space symmetry is for crystal structures.

One of the nontrivial assumptions, which initially appeared in the theory of spin glasses  \- , was that the energy landscapes of complex systems have a hierarchical structure. In recent years, this conjecture was confirmed by the results of numerous computer simulation studies (see, for example, ). It turned out that energy landscapes of complex systems represent hierarchically nested basins of local minima, viz., large basins consist of smaller basins of minima, the latter consist of even smaller ones, etc. The basins of minima are separated from one another by a hierarchy of activation barriers, i.e., the larger basins are separated by higher activation barriers, while the smaller basins within the larger ones are separated by smaller activation barriers, respectively.

The existence of the basins of local minima allows one to describe configurational dynamics in terms of the basin-to-basin kinetics . In this connection, it is assumed that transitions between the basins belonging to the same level of the basin’s hierarchy correspond to configurational rearrangements of a certain scale. Within such basins, a quasi-equilibrium distribution is established during relatively short (for a given scale) period of time, while transitions between the basins determine long-time configurational rearrangements. The hierarchy of nested basins corresponds to the hierarchy of scales of configurational rearrangements.

As indicated in the early works on the spin-glass theory, the hierarchy of nested basins of local minima is similar to the hierarchy of nested ultrametric spheres -. This idea can be found in some publications about anomalous relaxation in glass-like systems -. However, these investigations give no answer to the question whether ultrametricity can be directly used for the description of the dynamics of complex systems.

Recent studies  show that there is an approach which may lead to a noticeable progress in this direction. This approach is based on the $`p`$-adic analysis natural for ultrametric spaces. An introduction to the $`p`$-adic analysis can be found, for instance, in . Some applications of the $`p`$-adic analysis in theoretical physics and theoretical biology are described in -.

The $`p`$-adic master equation describing a random walk constrained by the hierarchical energy landscapes of a certain type (the so-called ultrametric diffusion) was constructed in . In these works, we describe a general method of finding its solutions and examine some Cauchy problems pertaining to conformational dynamics and reactions in proteins, and we also find some solutions of these problems expressed by analytic formulas. It is an interesting fact that the $`p`$-adic pseudo-differential equation introduced earlier as a $`p`$-adic analogue of the diffusion equation (see ) coincides with the particular $`p`$-adic master equation obtained in .

Regarding $`p`$-adic pseudo-differential equations and stochastic processes in non-Archimedean spaces, see , -.

The present paper is aimed at a further development of the approach suggested in  for the description of configurational dynamics of complex systems. We show that the Kohlrausch-Williams-Watts law, the power decay law, and the logarithmic decay law can be described by $`p`$-adic master equations of the same form, and therefore, the relaxation laws characteristic of complex systems actually reflect the same type of ultrametric random processes.

The paper is structured as follows. In Section 2, we introduce a general form of the $`p`$-adic master equation describing random walk constrained by a hierarchical energy landscape, and we also define a $`p`$-adic model of relaxation. In Section 3, we consider three versions of such models and establish their correspondence with three types of relaxation in complex systems.

## 2 Random walk in ultrametric space

The general $`p`$-adic master equation describing a Markovian process of random walk in ultrametric space can be written as follows:

$$\frac{f(x,t)}{t}=__p\left[w(x|y)f(y,t)w(y|x)f(x,t)\right]𝑑\mu (y),$$
(1)
where $`x`$ is a $`p`$-adic number, $`d\mu (x)`$ is the Haar measure on the field of $`p`$-adic numbers $`_p`$, and $`t`$ is time, $`t_+`$. Equation (1) is the usual balance equation for transitions between system states. The only peculiarity is that the space of states is described by the $`p`$-adic numbers, i.e. this space is ultrametric. The function $`f(x,t):_p\times _+_+`$ is a probability density distribution: the integral $`_Bf(x,t)𝑑\mu (x)`$ is the probability of finding the system in a domain $`B_p`$ at the instant $`t`$. The function $`w\left(x|y\right):_p\times _p_+`$ is the probability of the transition from the state $`y`$ to the state $`x`$ per unit time.

The transition from a state $`y`$ to a state $`x`$ can be perceived as overcoming the energy barrier separating these states. The structure of the relation between the transition probability $`w\left(x|y\right)`$ and the parameters of the energy barrier depends, in general, on the mechanism of configurational rearrangements, transition pathways, their distribution, etc. Finding this expression on the basis of the microscopic theory is a nontrivial problem beyond the scope of this article. For the sake of definiteness, we will use the standard Arrhenius relation

$$w\left(x|y\right)A\left(T\right)\mathrm{exp}\left\{\frac{U\left(x|y\right)}{kT}\right\}$$
(2)
where $`U\left(x|y\right)`$ is the height of the (effective) activation barrier for the transition from the state $`y`$ to the state $`x`$, $`k`$ is the Boltzmann constant, and $`T`$ is temperature. Studies of thermal behavior of anomalous relaxation in complex systems show that this expression is quite acceptable for a fairly wide range of temperatures.

Formula (2) establishes a relation between the structure of the energy landscape $`U\left(x|y\right)`$ and the transition function $`w\left(x|y\right)`$. In particular, the condition $`w\left(x|y\right)=w\left(y|x\right)`$ corresponds to a degenerate energy landscape. An important example of a degenerate landscape is a regular hierarchical landscape. In this case, on each hierarchical level, the basins of states split up into the same number of smaller basins; moreover, the activation barriers between the basins of the same hierarchical level have equal values. Graphically, this hierarchical landscape is represented as a regular tree with the same branching rules for all branching nodes. As shown in , the activation barriers $`U\left(x|y\right)`$ separating the states $`x`$ and $`y`$ depend, in this case, only on ultrametric distances between the states, $`\left|xy\right|_p=p^\gamma `$, $`\gamma `$, and therefore, regular energy landscapes are completely described by functions of the form $`U(\left|xy\right|_p)`$. Thus, in the case of random walk constrained by a regular hierarchical landscape, the master equation (1) takes the form

$$\frac{f(x,t)}{t}=__pW(\left|xy\right|_p)\left[f(y,t)f(x,t)\right]𝑑\mu (y)$$
(3)
where

$$W(\left|xy\right|_p)=\frac{A\left(T\right)}{\left|xy\right|_p}\mathrm{exp}\left\{\frac{U\left(\left|xy\right|_p\right)}{kT}\right\}.$$

The $`p`$-adic master equation (3) admits analytic solutions for a fairly wide class of functions $`W(\left|xy\right|_p)`$. The technique of finding its solutions is described, for instance, in .

It may seem that regular hierarchical landscapes represent a very idealistic model of real landscapes and have little relation to the latter. Indeed, regular hierarchical landscapes reflect no details of the energy landscapes of complex systems other than their hierarchical structure. Oddly enough, this is sufficient for the description of the characteristic types of relaxation.

On the basis of general physical considerations, one can identify three essentially different types of regular hierarchical landscapes:

“Logarithmic landscapes” are characterized by activation barriers having a slow (logarithmic) growth with respect to the number $`\gamma `$ of the hierarchical level, $`U(\gamma )\mathrm{ln}\gamma `$. This type of landscapes may be associated with “flexible” systems having a wide range of available scales for configurational rearrangements.

Contrary to logarithmic landscapes, “exponential landscapes” are characterized by a rapid (exponential) growth of activation barriers, $`U(\gamma )e^\gamma `$. Landscapes of this type represent “stiff” systems with a small set of distinct scales available for configurational rearrangements.

Finally, there are “linear landscapes” with linear growth of activation barriers, $`U(\gamma )\gamma `$. These landscapes may be associated with systems of an “intermediate” type.

In all these cases, we consider the Cauchy problem for the respective master equation (3) with the initial conditions:

$$f(x,0)=\mathrm{\Omega }\left(\left|x\right|_p\right)=\{\begin{array}{cc}1\hfill & \left|x\right|_p1\hfill \\ 0\hfill & \left|x\right|_p>1\text{ .}\hfill \end{array}$$
(4)
A relaxation process will be understood as the evolution of population in the domain of the initial distribution $`\left|x\right|_p1`$:

$$S(t)=\underset{\left|x\right|_p1}{}f(x,t)𝑑\mu (x).$$
(5)
In each of the above cases, the analytical estimation of relaxation function $`S(t)`$ is being found. It would enable the relation of $`S(t)`$ to three empirical laws for characteristic relaxation.

## 3 Characteristic relaxation types

The solution of the Cauchy problem for equation (3) with the initial condition (4) has the form

$$f(x,t)=\underset{n=0}{\overset{\mathrm{}}{}}e^{\stackrel{~}{W}(p^n)t}\left[(1p^1)p^n\mathrm{\Omega }(p^n\left|x\right|_p)p^{n1}\delta (\left|x\right|_pp^{1+n})\right]$$
(6)
where

$$\delta (\left|x\right|_pp^\gamma )=\{\begin{array}{cc}1\hfill & \left|x\right|_p=p^\gamma \hfill \\ 0\hfill & \left|x\right|_pp^\gamma \text{ .}\hfill \end{array}$$

Indeed, applying the $`p`$-adic Fourier transformation to (3) and taking into account the initial condition (4), we get

$$f(x,t)=__pe^{\stackrel{~}{W}(\left|k\right|_p)t}\chi (kx)𝑑\mu (k)$$

where

$$\stackrel{~}{W}(\left|k\right|_p)=__pW(\left|x\right|_p)(\chi (kx)1)𝑑\mu (x).$$
(7)
Calculating the integral and using the substitution $`\left|k\right|_p=p^n`$, we obtain

$$\stackrel{~}{W}(p^n)=p^{n+1}\left[(1p^1)\underset{m=0}{\overset{\mathrm{}}{}}p^mW(p^{m+n+1})+p^1W(p^{n+1})\right]$$

and therefore, (6) is the desired solution.

In order to find $`S(t)`$, the solution (6) should be substituted into (5). Calculating the corresponding integral, we find that

$$S(t)=(p1)\underset{n=1}{\overset{\mathrm{}}{}}p^ne^{\mathrm{\Lambda }_nt}.$$
(8)
Here, for the sake of convenience, we have introduced the positive function

$$\mathrm{\Lambda }_n=p^n\left[(1p^1)\underset{m=0}{\overset{\mathrm{}}{}}p^mW(p^{m+n})+p^1W(p^n)\right].$$
(9)
Next, we obtain estimates of $`S(t)`$ for the above three regular hierarchical landscapes. The transition probability $`W`$ is defined here by (2). The frequency factor is taken equal to the unity, and the temperature dependence is characterized by the parameter $`\alpha =T_0/T`$.

### 3.1 Logarithmic landscape.

In view of (2), the transition probability function for the logarithmic landscape can be represented as

$$W(\left|xy\right|_p)=\frac{1}{\left|xy\right|_p}\frac{1}{\mathrm{ln}^\alpha (1+\left|xy\right|_p)},\alpha >1.$$

In this case,

$$\mathrm{\Lambda }_n=(1p^1)\underset{m=0}{\overset{\mathrm{}}{}}\mathrm{ln}^\alpha (1+p^{m+n})+p^1\mathrm{ln}^\alpha (1+p^n).$$
(10)
Then, the following estimates hold:

$$(1p^1)\frac{\mathrm{ln}^\alpha p}{(\alpha 1)}\frac{1}{(n+1)^{\alpha 1}}<\mathrm{\Lambda }_n<\frac{\alpha \mathrm{ln}^\alpha p}{(\alpha 1)}\frac{1}{n^{\alpha 1}}.$$

Using these inequalities, we can estimate the relaxation function $`S(t)`$ as follows:

$`p\mathrm{exp}\left\{{\displaystyle \frac{2\alpha }{\alpha 1}}\left({\displaystyle \frac{t}{\mathrm{ln}p}}\right)^{\frac{1}{\alpha }}\right\}`$
$`<S(t)<p^2\left(\mathrm{exp}\left\{{\displaystyle \frac{1}{\alpha 1}}\left({\displaystyle \frac{t}{\mathrm{ln}p}}\right)^{\frac{1}{\alpha }}\right\}+\mathrm{exp}\left\{\left({\displaystyle \frac{t}{\mathrm{ln}p}}\right)^{\frac{1}{\alpha }}\right\}\right).`$
Thus, the long-time behavior of relaxation kinetics constrained by logarithmic landscapes is fitted into the Kohlrausch-Williams-Watts law.

It is interesting to note, that as $`\alpha 1`$, the upper bound for $`S(t)`$ approaches an exponent. For $`\alpha 1`$, the series (10) is divergent. This fact may be interpreted as a transition to exponential relaxation when a critical temperature $`T_0`$ is attained. This kind of critical behavior of ultrametric diffusion constrained by logarithmic landscape was previously predicted in  on the basis of some qualitative considerations.

### 3.2 Linear landscape.

In the case of a linear landscape, the transition probability function has the form:

$$W(\left|xy\right|_p)=\frac{1}{\left|xy\right|_p^{\alpha +1}}.$$

Estimates for the probability density distribution $`f(x,y)`$ in the case of a linear landscape were obtained in . Calculations for $`S(t)`$ similar to those of  yield the estimate

$$\frac{1}{p}\mathrm{\Gamma }\left(\frac{1}{\alpha }+1\right)\left(\mathrm{\Gamma }_p(\alpha )t\right)^{\frac{1}{\alpha }}<S(t)<\mathrm{\Gamma }\left(\frac{1}{\alpha }+1\right)\left(\mathrm{\Gamma }_p(\alpha )t\right)^{\frac{1}{\alpha }}$$

where $`\mathrm{\Gamma }(\alpha )`$ is the gamma-function and $`\mathrm{\Gamma }_p(\alpha )`$ is the $`p`$-adic gamma-function (see ).

Thus, the long-time behavior of relaxation kinetics constrained by linear landscapes is fitted into the power decay law.

### 3.3 Exponential landscape.

In this case, the transition probability function is

$$W(\left|xy\right|_p)=\frac{1}{\left|xy\right|_p}e^{\alpha \left|xy\right|_p},\alpha >0,$$

and from (9) we have

$$\mathrm{\Lambda }_n=\left(1\frac{1}{p}\right)\underset{m=0}{\overset{\mathrm{}}{}}\mathrm{exp}\left[\alpha p^{n+m}\right]+\frac{1}{p}\mathrm{exp}\left[\alpha p^n\right].$$

The following estimates hold:

$$\mathrm{exp}\left[\alpha p^n\right]<\mathrm{\Lambda }_n<\frac{p^\alpha }{p^\alpha 1}\mathrm{exp}\left[\alpha p^n\right].$$

Hence, for relaxation function $`S(t)`$ we get

$$\frac{p}{2e}\frac{\alpha }{\mathrm{ln}t+\mathrm{ln}\left(p^\alpha /(p^\alpha 1)\right)}<S(t)<(p^2+p+e)\frac{\alpha }{\mathrm{ln}t}.$$

The last estimate shows that the long-time behavior of relaxation kinetics constrained by exponential landscapes is fitted into the logarithmic decay law.

It has been shown that the Kohlrausch-Williams-Watts law, the power decay and the logarithmic decay laws reflect the same type of random processes. Such processes are ultrametric and can be described by the $`p`$-adic master equation formulated in . It seems that ultrametricity is a most general and fundamental feature of the dynamics of complex systems: The configurational space of such systems may be totally disconnected due to numerous quenched constraints, and therefore, it is impossible to accomplish a relatively large configurational rearrangement step-by-step, by way of small configurational changes.

In conclusion, we consider it necessary to make a few remarks.

A direct comparison of the function $`S(t)`$ with experimentally observed relaxation kinetics requires a certain caution. Indeed, the function $`S(t)`$ is the population of a region in the configurational space, whereas in most experiments, the focus of observation is some probe coupled to configurational rearrangements, rather than configurational rearrangements themselves. A specific form of the relation between the probe state and a local configuration may happen to be important for the description of experimental data.

The above three characteristic relaxation laws do not exhaust the diversity of relaxation kinetics in complex systems. There are numerous observations of ”mixed” types of kinetics, with one type of relaxation kinetics observed in some time-window and another type observed outside. For instance, it may happen that at low temperatures the observed kinetics is well fitted into the logarithmic decay law, and at higher temperatures it is described by the power decay law. In other cases, the power decay may transform into the Kohlrausch-Williams-Watts law. It is not difficult to imagine energy landscapes whose activation barriers near the ground state are close to the exponential landscapes, while higher activation barriers correspond to a linear landscape and then transform into a logarithmic landscape. The above approach can be easily used for the description of ultrametric diffusion constrained by such ”mixed” landscapes. Equation (3) can be solved analytically for a wide class of functions $`W(\left|xy\right|_p,\alpha )`$ admitting the $`p`$-adic Fourier transformation.

Finally, our considerations here are restricted to degenerate energy landscapes. These may be directly related only to some specific complex systems, for example, spin glasses. However, energy landscapes of other complex systems are obviously non-degenerate. Local minima in such landscapes are also clustered into a hierarchy of basins, and as the whole, they look like a global potential “hole” (or “funnel”) with extremely rugged “walls”. The description of random walk on such landscapes by $`p`$-adic equations is of great interest, particularly, in connection with the applications of the concept of hierarchical energy landscapes to the protein folding  (see also  and references therein). Some results regarding this problem will be presented in our further publications.

We would like to express our profound gratitude to Academician V.S.Vladimirov and Professor I.V.Volovich for fruitful discussions and valuable remarks. The authors are grateful to Dr. A.P.Zubarev and Dr. S.V.Kozyrev for their active participation in discussions. The given work has been partially supported by RFBR-00-15-97392, RFBR-00-15-96073, and INTAS-9900545 grants.

## References