text:
# Surfactants in semiconductor heteroepitaxy: Thermodynamics and/or kinetics?

## I INTRODUCTION

The presence of surface active species (surfactants) on the surface of growing crystals changes dramatically the mode of growth, the rate of nucleation, the transition from step flow to two-dimensional (2D) nucleation, the structure of steps, the distribution of stress, surface alloying, the concentration of defects, etc. The major effect is the suppression of the three-dimensional (3D) islanding so that smooth films suitable for fabrication of microelectronic devices could be produced. Although the initial attempts to understand the effect of the surfactants on the mode of growth were based on thermodynamic grounds, there is still a debate in the literature between the adherents of the thermodynamic and the kinetic views. That is why we will make an attempt to consider both the thermodynamic and the kinetic aspects of the effect of surfactants on the thin film growth mode. We will do that in the most general way without accounting for the particular properties of the materials. We note only that in the case of semiconductor growth a complete monolayer of the surfactant is needed.

## II THERMODYNAMICS

If we want to change the direction of a process we have to change the sign of the thermodynamic force that drives it. If we want to prevent a particular process to take place we have to suppress the corresponding thermodynamic driving force making it as close to zero as possible. The thermodynamic driving force which determines the occurrence of one or another mechanism of epitaxial growth is the difference $`\mathrm{\Delta }\mu =\mu (n)\mu _{3D}^0`$ of the chemical potential, $`\mu (n)`$, of the overlayer which depends on the film thickness measured in number $`n`$ of monolayers counted from the interface, and the chemical potential, $`\mu _{3D}^0`$, of the bulk 3D crystal. The thickness dependence of the film chemical potential $`\mu (n)`$ originates mostly from the thickness distribution of the strain due to the lattice misfit, but the interaction between the deposit and the substrate, which rapidly decreases with the distance from the interface and could be neglected beyond several monolayers should be also accounted for. If we deposit a crystal $`A`$ on the surface of a crystal $`B`$ in absence of a surfactant the thermodynamic driving force can be written in terms of surface energies $`\sigma _A`$ and $`\sigma _B`$

$$\mathrm{\Delta }\mu =a^2[\sigma _A+\sigma _{AB}(n)\sigma _B]$$
(1)
where $`a^2`$ is the area occupied by an atom, and $`\sigma _{AB}(n)`$ is the interfacial energy which includes in itself the misfit strain energy and the attenuation of the energetic influence of the substrate. In fact, this is the familiar 3-$`\sigma `$ criterion of Bauer.

We can write the above expression in terms of interatomic energies

$$\mu (n)=\mu _{3D}^0+(E_{AA}E_{AB})=\mu _{3D}^0+E_{AA}\Phi $$
(2)
where $`E_{AA}`$ and $`E_{AB}`$ are the energies per atom to disjoin a half-crystal $`A`$ from a like half-crystal $`A`$ and an unlike half-crystal $`B`$. It is in fact the adhesion energy $`E_{AB}`$ which includes in itself the thickness distribution of the strain energy. In the above equation $`\Phi =1E_{AB}/E_{AA}`$ is the so called adhesion parameter which accounts for the wetting of the overgrowth by the substrate. This is the same parameter which enters the work of formation of nuclei on a foreign substrate in the classical nucleation theory.

Equations (1) and (2) are identical. Eq. (1) can be readily obtained from (2) by using the definition of surface free energy and the relation of Dupré. The essential physics, however, is in Eq. (2) as the surface energies are derivatives of the interatomic forces. It shows that the thermodynamic driving force for occurrence of one or another mode of growth is the relative adhesion/cohesion difference, or in other words, the wetting. The incomplete wetting, $`0<\Phi <1`$, is the driving force for 3D nucleation and growth, whereas the complete wetting, $`\Phi 0`$, is the driving force for planar growth by nucleation and growth of 2D monolayer height nuclei.

When $`0<\Phi <1`$ $`\mu (n)>\mu _{3D}^0`$, and 3D islands are formed on top of the substrate surface from the very beginning of deposition. This is the well known Volmer-Weber (VW) growth and $`\mu (n)`$ goes asymptotically to $`\mu _{3D}^0`$ from above (the curve denoted by VW in Fig. 1).

In the other extreme of complete wetting ($`\Phi 0`$) either layer-by-layer (Frank-van der Merwe or FM) growth at negligible lattice misfit, or a Stranski-Krastanov (SK) growth (layer-by-layer growth followed by 3D islands) at a perceptible misfit, are thermodynamically favored. In the case of FM growth $`\mu (n)`$ goes asymptotically to $`\mu _{3D}^0`$ from below (the curve denoted by FM in Fig. 1). The first monolayer has the lowest chemical potential owing to the strongest interaction with the substrate, and a second monolayer can form only after the completion of the first one.

The SK growth mode appears as a result of the interplay of the film - substrate bonding, strain, and surface energies. A wetting layer consisting of an integer number of monolayers is first formed by a FM mode which is driven by a negative thermodynamic driving force $`\mathrm{\Delta }\mu `$ (complete wetting). Further FM growth becomes unfavorable as the strain energy accumulates linearly with film thickness. In addition, the stronger attraction from the substrate, that overcompensates the strain energy, disappears beyond several atomic diameters. On top of the wetting layer 3D islands form and grow under the influence of a positive thermodynamic driving force (incomplete wetting). If misfit dislocations (MDs) are introduced immediately after the completion of the wetting layer unstrained 3D islands are formed and grow on top of it. This is the classical Stranski-Krastanov growth in which the work needed to create new surfaces is overcompensated by a complete strain relaxation. It was found recently that initially coherent (dislocation-free) 3D islands are formed (coherent Stranski-Krastanov growth), in which the formation of a new surface is overcompensated by a gradual strain relaxation. In both cases the 3D islands and the wetting layer represent necessarily different phases separated by an interphase boundary. This boundary is determined by the displacements of the atoms belonging to the first atomic plane of the islands from the bottoms of the corresponding potential troughs of the uniformly strained wetting layer. In the classical SK growth the displacements are described in terms of MDs, whereas in the coherent SK growth the atoms that are closer to the islands edges are displaced owing to the islands finite size. These displacements give rise to weaker adhesion, or in other words, to incomplete wetting. It is namely this incomplete wetting (or weaker adhesion) which appears as the thermodynamic driving force for the 3D islanding in the SK mode of growth. In other words, we can treat the SK mode as a FM mode driven by complete wetting ($`\mathrm{\Delta }\mu <0`$), followed by VW mode driven by incomplete wetting ($`\mathrm{\Delta }\mu >0`$) (the curve denoted by SK in Fig. 1).

When considering the effect of the surfactant on the mode of growth we note that the chemical potential $`\mu _{3D}^0`$ of the infinitely large crystal does not depend on whether there are impurities adsorbed on its surface or not. This, however, is not true for sufficiently thin films. If a monolayer thick film of $`A`$ is covered by a monolayer of a surfactant (S) atoms the growth mode criterion (1) turns into

$$\mathrm{\Delta }\mu _S=a^2(\sigma _{AB}+\sigma _{SA}\sigma _{SB})$$
(3)
where the $`\sigma `$’s are the corresponding interfacial energies. This expression is equivalent to

$$\mu _S(n)=\mu _{3D}^0+(E_{AA}E_{AB})(E_{SA}E_{SB})$$
(4)
where $`\mu _S(n)`$ is the chemical potential of the film covered by S atoms. It is immediately seen that if 3D islanding is thermodynamically favored in absence of a surfactant, i.e. $`E_{AA}>E_{AB}`$, the inequality (4) may have the opposite sign if the adhesion of the surfactant is stronger to the overlayer rather than to the substrate and thus the third term in (4) overcompensates the second one. Then the surfactant is ”good” from thermodynamic point of view if it adheres more strongly to the overlayer, thus changing the sign of the thermodynamic driving force.

This is valid, however, only for the transition from VW to FM growth as is the case of the As mediated growth of Si on Ge(001). Even in this case it is valid only at zero misfit. After several monolayers of $`A`$ the energetic contact between the surfactant and the substrate $`B`$ is completely lost ($`E_{SB}E_{SA}`$), and if the misfit is large 3D islands should appear resulting in SK mode. In systems that follow the SK growth mode where, beyond the wetting layer, $`A`$ is deposited on strained $`A`$, the difference $`E_{SA}E_{SB}`$ will be nearly equal to zero. Many papers both theoretical and experimental have been devoted to the problem considering mainly the surfactant growth of Ge on Si(001). However, in most of the papers quoted above the mechanism of the exchange of the surfactant (As or Sb) atoms with Ge atoms belonging to the first monolayer on the Si(001) surface have been studied. It is obvious, that such studies cannot explain the effect of the surfactant to suppress the 3D islanding as the first monolayer belongs to the wetting layer and it will grow in a FM mode pseudomorphous with the substrate irrespective of there is a surfactant on top or not.

Obviously, in the SK mode the surfactant should make the thermodynamic driving force for 3D islanding on top of the wetting layer equal to zero for considerable film thickness. As mentioned above in the SK mode the atoms which are in contact with the wetting layer should be displaced from their ”bulk-like” positions either to make MDs or near to the islands edges. These displacements make the wetting incomplete, or which is the same, make the adhesion weaker than the cohesion ($`E_{AB}<E_{AA}`$), which is the thermodynamic driving force for 3D islanding. As found in a series of papers the surfactant group-V atoms form dimers on Si(Ge)(001) which reside exactly where the Si(Ge) dimers should be located. After the surfactant and overlayer dimers exchange places the latter occupy the epitaxial sites that are provided by the outgoing surfactant. The surfactant dimers on top do not permit large displacements of the overlayer atoms underneath, and in turn incomplete wetting. In other words, the surfactant atoms press the overlayer atoms into bulk-like positions and do not allow the formation of misfit dislocations. Thus, the surfactants suppress the thermodynamic driving force for 3D islanding by not allowing large atomic displacements. Once the thermodynamics do not require 3D islanding a change of the growth mode with the temperature should not be observed.

## III KINETICS

In absence of a surfactant the atoms arriving at the crystal surface diffuse on it, join pre-existing steps or islands, or give rise to new islands. At sufficiently high temperatures the diffusivity of the atoms on the terraces is high and they reach the pre-existing steps before meeting with each other. This results in a step-flow growth. At low temperatures and in turn low diffusivity the atoms meet with each other before reaching the steps and give rise to 2D nuclei on the terraces. In presence of a complete monolayer of S atoms two new phenomena take place. First, the surfactant changes the energetics of incorporation of atoms into ascending and descending steps, and second, the overlayer atoms have to exchange places with the S atoms on the terraces or at the steps in order to join the crystal lattice. Both phenomena lead to a change of the diffusivity of the adatoms across the steps and on terraces, and in turn exert an effect on the direction of transport (upward or downward) of atoms. Moreover, the change of diffusivity on terraces affects strongly the kinetics of nucleation and growth of the islands thus stimulating either the step flow or 2D nucleation. We consider in more detail the effect of the surfactant on the kinetics of attachment and detachment of atoms to and from the steps, and the exchange-deexchange kinetics on terraces.

Whereas an atom approaching an ascending step joins it upon striking, an atom approaching a descending step has to overcome an additional Ehrlich-Schwoebel (ES) barrier. If the ES barrier is low the overlayer grows more or less in a layer-by-layer mode. Otherwise, instabilities of different kind appear. Zhang and Lagally pointed out the possibility of a reduction of the ES barrier when overlayer atoms exchange sites with S atoms that decorate the steps, in such a way that the last barrier before descending is shifted downwards below the level of the surface diffusion barrier. The latter is in good agreement with first principle molecular dynamics calculations of displacement of Sb dimer at the step edge by a Si dimer through a push-out mechanism. In addition, Markov has suggested that atoms that approach ascending steps have to overcome an extra energetic barrier owing to the necessity to displace the S atoms already adsorbed at the step edges. These effects which are of purely kinetic origin can suppress the thermodynamic driving force for 3D islanding in heteroepitaxial growth.

When thermodynamics predict 3D islanding, ($`\mathrm{\Delta }\mu >0`$) the atoms are more strongly bound to the upper layers rather than to the lower layers. In other words the upper layers have lower chemical potentials (Fig. 1). This results in a gradient of the chemical potential that drives the atoms upwards and is the driving force for the 2D-3D transformation. The appearance of an additional barrier at the ascending steps, and the reduction of the ES barrier at the descending steps owing to the decoration of the steps by S atoms, reverse the asymmetry of the attachment - detachment kinetics and thus lead to a suppression of the thermodynamic driving force for 2D-3D transition. A reverse gradient of adatoms of a kinetic origin appears that drives the atoms downwards. The latter results in a planar rather than a 3D growth. Thus the attachment - detachment kinetics suppress the thermodynamic driving force for 3D islanding. At high temperature the thermodynamics prevail and SK growth takes place irrespective of the presence of the surfactant. Surfactant induced FM growth occurs at lower temperatures. Note that this kinetic effect does not require a complete monolayer of the surfactant but a small amount which is sufficient to decorate the steps.

In the presence of a complete monolayer of the surfactant a new kinetic effect appears owing to the necessity of the overlayer atoms to join the crystal lattice and of the S atoms to float over the surface. These are the phenomena of exchange and deexchange of overlayer and S atoms. Zhang and Lagally first assumed that an exchange process between overlayer and S atoms on the terraces should take place on the S precovered surfaces. Kandel and Kaxiras assumed later that a deexchange process is also possible and can play a significant role in epitaxial growth. They calculated the values of 0.8 and 1.6 eV for the activation energies, $`E_{ex}`$ and $`E_{dex}`$, for exchange and deexchange processes, respectively. Ko, Chang and Yi computed through first-principles pseudopotential total-energy calculations the values of the barriers for surface diffusion, exchange and deexchange for Ge on Ga, As and Sb precovered Si(111) surface. They found that in all cases $`E_{dex}>E_{ex}`$ but the differences are not as large as the ones calculated by Kandel and Kaxiras.

We introduce two time constants which characterize the exchange and the deexchange processes both being normalized to the time, $`t_1=1/F`$, necessary for deposition of a complete monolayer, $`F`$ being the the atom arrival frequency. First, this is the mean residence time of the atoms on top of the S layer before exchange $`\tau _{ex}=(F/\nu )\mathrm{exp}(E_{ex}/kT)`$ where $`\nu `$ is the attempt frequency. Second, this is the mean residence time before deexchange of an atom embedded in (or under) the S layer $`\tau _{dex}=(F/\nu )\mathrm{exp}(E_{dex}/kT)`$. Obviously, if a particular time constant is greater than unity the corresponding process will not occur. As $`\tau _{dex}\tau _{ex}`$ three possibilities exist. The first is $`\tau _{dex}1`$. This extreme describes the reversible exchange as the adatoms have sufficient time to go back on top of the S layer through a deexchange process. A dynamic exchange - deexchange equilibrium is established. With the value of $`E_{dex}=1.6`$ eV calculated by Kandel and Kaxiras we find, at 600K, $`\tau _{dex}0.002`$. Then the exchange - deexchange equilibrium is established in the very beginning of the deposition process. The nucleation process takes place in a more or less disordered (because of the different size of the atoms involved) 2D phase - a monolayer consisting of mixed S and overlayer atoms. A considerable fraction of the incoming atoms remains on top of the S layer and diffuse fast to the pre-existing steps. The formation, growth and decay of nuclei take place through exchange and deexchange processes. The second extreme is $`\tau _{ex}1<\tau _{dex}`$. This is the case of irreversible exchange. The incoming atoms rapidly exchange places with S atoms and remain buried under the S layer. All processes of diffusion, nucleation and incorporation into islands and steps occur under (or between) the S atoms. The case $`1<\tau _{ex}\tau _{dex}`$ should be excluded as it means that the S atoms will remain buried under the arriving overlayer atoms.

We calculate further the nucleus density, $`N_S`$, in both cases of reversible and irreversible exchange. This is extremely important because of two reasons. First, the large number of small islands will promote layer-by-layer growth owing to the existence of a critical island size for second layer nucleation. Second, higher nucleation rate leads to growth by 2D nuclei, whereas in the opposite case step-flow growth will take place. In all cases the expression for $`N_S`$ can be written in the form

$$N_S=N_0\mathrm{exp}\left(\frac{\chi }{i}E_S/kT\right)$$
(5)
where $`E_S`$ combines all energy contributions which depend on the presence of the surfactant, and

$$N_0=\left(\frac{\nu }{F}\right)^\chi \mathrm{exp}\left(\frac{\chi }{i}E_0/kT\right)$$
(6)
is the nucleus density in absence of a surfactant ($`E_S=0`$), Thus, the island density scales with the ratio $`\nu /F`$, the scaling exponent $`\chi `$ being a function of the number $`i`$ of atoms in the critical nucleus. In the case of low barriers for attachment of atoms to 2D nuclei, the latter ”feel” each other through the diffusion fields around them, and

$$\chi =\frac{i}{i+2}.$$
(7)
Kandel has recently shown by using a rate equation approach that in presence of surfactants which give rise to barriers for incorporation of atoms to the critical nuclei

$$\chi =\frac{2i}{i+3}.$$
(8)
The scaling exponent (8) has been independently derived by Markov by using a different approach. He showed that it is characteristic for kinetic regime of growth where the adatom concentration between the islands is practically constant and the islands ”do not feel” the presence of each other.

As has been noted in Ref. () the scaling exponent (7) varies with $`i`$ from 1/3 to 1, whereas (8) has values larger than unity already at $`i>2`$. Thus, one could distinguish between diffusion and kinetic regimes of growth if $`\chi `$ is smaller or greater than unity. Typical examples of the scaling exponent (8) are the homoepitaxy of Si on Sn precovered surface of Si(111), and of Ge on Pb precovered surface of Si(111).

It is immediately seen that the exponential multiplying $`N_0`$ in Eq. (5) can be smaller or greater than unity depending on the interplay of the energies involved in $`E_S`$. If $`E_S<0`$ the surfactant will stimulate the nucleus formation ($`N_SN_0`$). The latter means that a greater density of smaller islands will be formed, and the film will grow by 2D nucleation rather than by step-flow. What is more important is that the greater density of smaller nuclei leads to layer-by-layer growth rather than to the mound formation. In the opposite case the surfactant will inhibit the nucleation and will promote the step-flow growth. In the case of reversible exchange $`E_S`$ contains the term $`i[(E_{dex}E_{ex})\mathrm{\Delta }E_{sd}]`$ where $`\mathrm{\Delta }E_{sd}=E_{sd}^{}E_{sd}`$ is the the difference between the barriers for diffusion on the clean substrate, $`E_{sd}^{}`$, and on the surfactant passivated surface, $`E_{sd}`$, ($`E_{sd}^{}>E_{sd}`$). If $`\mathrm{\Delta }E_{sd}`$ overcompensates the difference $`E_{dex}E_{ex}`$, we could expect that $`E_S>0`$, the nucleus formation will be inhibited, ($`N_SN_0`$), and surfactant-induced step-flow growth will take place.

In the case of irreversible exchange Eq. (5) is still valid but the barriers $`E_{dex}`$ and $`E_{ex}`$ do not enter the energy term $`E_S`$. The largest term in $`E_S`$ is the increment of the surface diffusion barrier $`i(E_{sd}E_{sd}^{})`$ where now $`E_{sd}>E_{sd}^{}`$ owing to the fact that the process of diffusion under the S layer is inhibited compared with that on the clean surface. It is unlikely this term to be overcompensated by the smaller positive terms, and hence, we can always expect that the surfactant will stimulate the 2D nucleation rather than the step-flow growth.

Finally, we discuss briefly the reasons that lead to a change of the island density. Usually the latter is attributed to the change of atoms diffusivity. However, the change the island density might be due to other reasons. The nucleation rate depends on the Gibbs free energy for nucleus formation which is in turn a function of the nucleus size. The presence of the surfactant could lead to a change of the number of atoms in the critical nucleus. As a result the values of $`E_S`$ and $`E_0`$ in Eqs. (5) and (6), as well as the scaling exponent $`\chi `$ will change in one and the same temperature interval. In addition, the regime of growth can change from diffusion to kinetic, so that the scaling exponent (7) will be replaced by (8). In this case, the island density will increase by orders of magnitude even if $`i`$ remains constant.

## IV GROWTH OF Ge ON Si

We discuss in more detail a ”model” system such as Ge/Si. As pointed out by Kaxiras, one and the same surfactant can have qualitatively different effect depending on the crystallographic orientation of the substrate. On Si or Ge(001) the surfactants form dimers which do not break during the exchange process. On (111) surface three different geometries (substitutional, trimers, and zig-zag chains) of the group-V adsorbates (P, As, Sb) are possible. In the first $`1\times 1`$ structure one monolayer of S atoms is incorporated into the upper part of the (111) bilayer, whereas in the last two structures a monolayer of S atoms is bonded on top of a complete (111) bilayer. In addition, the growth of the two surfaces of Si or Ge differs drastically. Whereas the (001) surface grows by single monolayers, the (111) surfaces grow by bilayers. Thus, we have completely different mechanisms of exchange of S and overlayer atoms. In the case of (001) surfaces the growth unit is a dimer so that an exchange of S and overlayer dimers takes place through a single monolayer. On (111) surfaces the mechanism of exchange depends on the surface geometry of the surfactant. For example, if the surfactant geometry is substitutional we could suppose that single overlayer atoms should displace single S atoms from their positions in the upper monolayer that constitutes the bilayer. After the complete building of the upper monolayer atoms belonging to the next lower monolayer should be incorporated first in order to serve as a template for the S atoms in the next upper monolayer. Thus one surfactant monolayer should be displaced by a complete bilayer of the overgrowth the process taking place by exchange of single atoms. Comparing both surfaces we could speculate that on (111) surface the surfactant could not suppress successfully the thermodynamic driving force for 3D islanding by not allowing large atomic displacements (formation of misfit dislocations).

We consider first the growth of Ge on Si(111). Hwang, Chang and Tsong studied this system by using Pb as a surfactant. They established that the Pb atoms make a substitutional $`1\times 1`$ structure on top of the Ge surface. It was found that the saturation nucleus density in the submonolayer regime scales with the deposition flux with $`\chi =1.76`$. This clearly shows that in the Pb mediated epitaxy of Ge on Si(111) the scaling exponent is given by Eq. (8) and the growth proceeds in a kinetic regime. In other words, significant energy barriers appear for attachment of atoms to the critical nuclei owing to the necessity of exchange of overlayer and S atoms. The same scaling exponent ($`\chi =1.8`$) has been obtained in the homoepitaxial growth of Si on Sn precovered surface of Si(111), after comparing the experimental data for the critical terrace width for step flow growth with the corresponding theory. Under clean conditions in the same system Voigtländer et al. obtained a scaling exponent smaller than unity ($`\chi =0.85`$).

The observations of Hwang et al. show unambiguously that Pb has a prominent kinetic effect on the submonolayer growth of Ge on Si(111). In a further study the same authors found that 3D islanding is suppressed below 470K. At temperatures higher than 470K 3D islanding takes place. This observation could be explained by the suppression of the thermodynamic driving force for 3D islanding by the attachment - detachment kinetics mentioned above. The upward surface flux is hindered whereas the downward flux is enhanced. We could conclude that the suppression of 3D islanding in Ge growth on Pb precovered Si(111) is due more to kinetic rather than to thermodynamic reasons.

Voigtländer and Zinner studied the growth of Ge on Si(111) mediated by Sb. They established the same substitutional structure of the Sb atoms on the surface of Ge. They found that a transition from surfactant mediated layer-by-layer growth to the equilibrium SK mode takes place above twice higher temperatures (900K) as compared with the Pb mediated growth. The higher temperatures could be related to stronger chemical bonding but qualitatively the surfactant effect is the same.

The system Ge/Si(001) is the most studied one. Ide has observed that after annealing at about 800K for 10 min the 2D Ge islands formed in submonolayer regime on As precovered Si(001) did not change whereas the islands deposited under clean conditions disappeared. This observation shows that 2D islands formed under the surfactant layer are more stable than those without S atoms on top of them. Antimony mediated deposition and annealing experiments carried out by Osten et al. showed that Sb prevents 3D islanding, but also can smooth an already islanded film. Katayama et al. found that as the coverage of the surfactant increases, the intermixing of Ge and Si, as well as the nucleation and the growth of macroscopic Ge islands, are suppressed. Both observations could be explained with the forced occupation of epitaxial sites by the Ge atoms under the influence of the S dimers on top of them. On the other hand, Ide also established that at 800K the Ge atoms reach and join the preexisting steps of Si (step flow growth) on the clean surface, whereas 2D islands are formed in the presence of As. The latter unambiguously demonstrates the kinetic effect of the surfactant on the mode of growth by decreasing the diffusivity of the incoming atoms.

Most of the studies of surfactant mediated growth of Ge on Si(001) surface are carried out at or near to 500C. To the author’s knowledge higher temperature deposition has been carried out by Tromp and Reuter (630C), and by Horn von Hoegen et al. at 700C. No 3D islanding except for mound formation has been noticed on top of the wetting layer at the highest temperature used. A change of the growth mode with temperature as in the case of Ge/Si(111) has not been reported so far. This can be interpreted as follows. The surfactant suppresses successfully the thermodynamic driving force for 3D islanding by not allowing large atomic displacements and the growth continues in a layer-by-layer mode (with the unavoidable roughness to relieve the strain) untill MDs are introduced at the interface. Once the thermodynamics do not require 3D islanding a change of the growth mode with the temperature should not be observed.

Jenkins and Srivastava carried out first principle density functional theory (DFT) calculations of the mode of growth of Ge on Sb precovered Si(001). They studied the behavior of the chemical potential of the Ge film relative to the bulk chemical potential, $`\mu _{3D}^0`$, in absence and presence of Sb as a surfactant. They found that in absence of Sb the chemical potentials of the first three monolayers of Ge, which belong to the wetting layer, are smaller (more negative) than $`\mu _{3D}^0`$ (complete wetting). The chemical potential of the forth monolayer becomes more positive than $`\mu _{3D}^0`$ (incomplete wetting) which is an indication of 3D islanding after the first three monolayers. In presence of Sb the chemical potential of the first monolayer is slightly more negative that $`\mu _{3D}^0`$ but the chemical potentials of the remaining three monolayers are equal to $`\mu _{3D}^0`$ within the accuracy of the calculations. An error is introduced in the calculations by not allowing the Ge overlayers to relax laterally. Thus after the exhange of Ge and Sb dimers the Ge atoms occupy the bulk-like positions left by the outgoing Sb dimers and the surfactant successfully suppressed the thermodynamic driving force for 3D islanding. In fact, the surfactant replaced the curve SK by the curve FM in Fig. (1).

Similar first principle DFT calculations have been carried out by Gonzá- les-Méndez and Takeuchi for the growth of Si on Ge(001) by using As as a surfactant. In absence of a surfactant Si grows in VW mode as 3D islands directly on the Ge surface. In presence of As the overlayer prefers to grow in a layer-by-layer mode due to thermodynamic reasons, the energy of the 3D islands being always lower than that of the complete monolayers. The deposit atoms are almost in bulk-like positions. Equally important, the authors established that As greatly reduces the Si and Ge intermixing. This is another indication of the thermodynamic influence of the surfactant as the alloying also requires large atomic displacements.

## V CONCLUSION

The presence of a third element in the system unavoidably changes the kinetics of growth. Hence, we have to decide whether the kinetics is solely responsible for suppression of the 3D islanding, or the thermodynamics plays a decisive role, with the kinetics exerting an additional effect. We argue that this question could be answered by studying the temperature dependence of the capability of the surfactant to suppress 3D islanding. If the kinetics play a decisive role we should expect a change of the mode of growth from layer-by-layer at lower temperature to 3D (SK or VW) islands at higher temperature. The physical reason is that at high temperature the system is closer to equilibrium whereas at lower temperature the kinetic effects suppress the thermodynamic driving force for 3D islanding. This is most probably the case of the system Ge/Si(111) where one monolayer of the surfactant should be displaced by two monolayers of the overgrowth and the exchange process takes randomly place by single atoms. The S atoms do not press the overlayer atoms into bulk-like positions and formation of misfit dislocations or displacements of the edge atoms is not inhibited. We have seen that the growth is accompanied by typical kinetic phenomena as the kinetic regime of growth ($`\chi >1`$). On the contrary, the system Ge/Si(001) demonstrates the decisive role of the thermodynamics. The kinetics play an auxiliary role in the same direction. The physical reason for suppression of the 3D islanding is that a monolayer of the surfactant is displaced by a monolayer of the film during growth, and the exchange process occurs by dimers rather than by single atoms. As a result the overlayer atoms occupy nearly bulk-like positions provided by the outgoing surfactant dimers. Large atomic displacements are forcibly prohibited which in turn leads to complete wetting and a zero thermodynamic driving force for 3D islanding. It is thus clear that the third element in the system does not act as a surfactant at all.

We conclude that the relative weight of the effect of the third element on the thermodynamics and kinetics in semiconductor heteroepitaxy depends mostly on the mechanism of exchange of overlayer and surfactant atoms. Undoubtedly the third element changes the surface energies but it is the change of the relative adhesion of the overlayer to the substrate that is crucial from thermodynamic point of view.

###### Acknowledgements.

The author is greatly indebted to M. Paunov for the fruitfull discussion and the critical reading of the manuscript.


format: text

text:
# SYNCHRONIZATION OF KAUFFMAN NETWORKS

## I Introduction

Synchronization of coupled elements is a form of collective evolution present in a variety of complex real systems and mathematical models. This class of emergent behavior has been observed in biological populations , chemical reactions , neural networks , and human social phenomena , among other instances. Models that account for synchronization consider, for example, globally coupled logistic maps , chaotic oscillators , Hamiltonian systems , and formal neural networks .

In the usual formulation, two identical dynamical systems whose individual dynamics is governed by the equation $`\dot{𝐰}=𝐅(𝐰)`$ are coupled to each other in the form

$$\dot{𝐰}_{1,2}=𝐅(𝐰_{1,2})+ϵ\left(𝐰_{2,1}𝐰_{1,2}\right),$$
(1)
where $`ϵ`$ is the coupling parameter. Full synchronization takes place when both systems converge asymptotically to a common trajectory, $`𝐰_1(t)=𝐰_2(t)`$. When the individual dynamics is chaotic—a particularly relevant case in connection with the description of real systems—full synchronization occurs above a critical value $`ϵ_c`$ of the coupling intensity. This critical point is determined by the competition between chaos and coupling, and can be calculated in terms of the Lyapunov exponent of the individual dynamics .

While globally coupled chaotic elements with a few internal variables have been extensively studied, synchronization of spatially extended systems remains quite unexplored. Recently, synchronization has been reported for a system consisting of two coupled complex Ginzburg-Landau equations . Globally coupled neural networks , stochastically coupled cellular automata , and nonidentical complex Ginzburg-Landau systems , are other examples of spatially extended systems that present a critical transition to synchronization.

In this paper we study the synchronization dynamics of two coupled identical Kauffman networks, which are discrete extended dynamical systems with quenched disorder. With respect to previous work on coupled extended systems, the interest of Kauffman networks resides in the fact that the transition to synchronization admits an analytical description which results to be in excellent agreement with numerical simulations. In Sect. II we briefly review the definition of Kauffman networks and the annealed model for the calculation of their overlaps. Next, in Sect. III, we introduce a stochastic coupling mechanism for Kauffman networks and propose an analytical approach in the framework of the annealed model, which identifies the transition to synchronization in our system as a transcritical bifurcation. Section IV, where we report our numerical results, is the core of the present paper. There, we study the effects of spurious synchronization in finite-size networks, consider the application of noise to the system to eliminate such effects, and compare the results with the analytical description. Remanent finite-size effects are numerically quantified and their analytical treatment, which requires a formulation beyond the annealed model, is outlined. In Sect. V we discuss the synchronization transition in some subclasses of Kauffman networks, which may be thought of as interpolations between generic Kauffman networks and ordered cellular automata. Finally, in Sect. VI, our results are summarized and discussed.

## II Kauffman networks

Kauffman networks, also known as random Boolean networks, were introduced as a model for the problem of cell differentiation . Since then, they have been the object of many studies concerning their properties .

A Kauffman network (KN) is a disordered deterministic dynamical system. It consists of an $`N`$-site network, where each site is connected to $`K`$ randomly chosen sites. The parameter $`K`$ is known as the connectivity of the network. We refer to the set of $`K`$ sites connected to a given site as its neighborhood. The state of each site is given by a Boolean variable $`\sigma _i\{0,1\}`$, and evolves according to the inputs coming from the neighbor sites. The evolution rule is chosen independently and randomly for each site. To each possible configuration of the neighborhood—there are $`2^K`$ such configurations—an output is assigned, namely, $`1`$ with probability $`p`$, or $`0`$ with probability $`1p`$. The parameter $`p[0,1]`$ is known as the bias of the rule. Then, for each variable $`\sigma _i`$ we have Boolean functions $`f_i`$ such that $`\sigma _i(t+1)=f_i\left[𝝂_i(t)\right]`$, where $`𝝂_i=(\sigma _{i_1},\mathrm{},\sigma _{i_K})`$ is the set of inputs of site $`i`$. The state of all sites is updated simultaneously according to the corresponding functions. We can write an evolution equation for the state vector of the network $`𝝈=(\sigma _1,\sigma _2,\mathrm{},\sigma _N)`$, as

$$𝝈(t+1)=𝐟\left[𝝈(t)\right],$$
(2)
with $`𝐟\left[𝝈(t)\right]=(f_1\left[𝝂_1(t)\right],f_2\left[𝝂_2(t)\right],\mathrm{},f_N\left[𝝂_N(t)\right])`$.

The $`K`$ connections and the evolution rule of each site are chosen at the beginning and kept fixed during the evolution. Thus, the disorder is quenched and the dynamics is deterministic. For a finite number of sites $`N`$, the number of states in phase space is also finite—it equals $`2^N`$. Then, for any initial condition, the system will eventually fall into a cycle.

In the $`(p,K)`$ parameter space, Kauffman networks present phases of frozen and chaotic evolution, separated by a critical line. The transition between these phases has been extensively studied, and characterized by means of several order parameters, such as the Hamming distance  and the stable core size . In most of this work we will deal with the case $`p=1/2`$ and $`K=3`$, which lies within the chaotic phase.

The annealed model (AM) was introduced to study the evolution of overlaps between states in KNs . In this model, the $`K`$ connections $`\{i_1,\mathrm{},i_K\}`$ of each site as well as the Boolean functions $`f_i`$ are randomly changed at each time step. This means that an entirely different realization of the network is used at each step. Note that, while ordinary KNs are deterministic, the annealed model works as a probabilistic automaton. The asymptotic periodic behavior of KNs is absent in the annealed model. The advantage of this model is that it allows for analytical calculations, and it has been shown that its predictions are in good agreement with the behavior of KNs in the limit $`N\mathrm{}`$ .

Suppose that we have two identical KNs with the same connections and rules. We feed them with different initial conditions, and let them evolve in time. We define the overlap $`a(t)`$ between the networks as the fraction of homologous sites that are in the same state at time $`t`$. In the AM, it is possible to calculate the time evolution of the overlap. At time $`t+1`$ the connections and the rules $`f_i`$ are reassigned, but the same changes are applied to both networks, keeping them identical. The probability for a site having all its inputs coming from sites in the same state in both networks is $`a(t)^K`$. At the next time step, consequently, such site will be in the same state in both networks, no matter the evolution rule chosen for it. Thus, there is a fraction $`a(t)^K`$ of homologous sites whose state will coincide at $`t+1`$. The remaining $`1a(t)^K`$ sites still have a probability of overlapping. Even if the state of the neighborhoods of a given site are different in the two networks, it may happen that the evolution rule assigns the same output to them. The probabilities for $`f_i\left(𝝂_i^1\right)=f_i\left(𝝂_i^2\right)=0`$ and for $`f_i\left(𝝂_i^1\right)=f_i\left(𝝂_i^2\right)=1`$ are, respectively, $`(1p)^2`$ and $`p^2`$. The overlap at time $`t+1`$ is then

$$a(t+1)=a(t)^K+\left[1a(t)^K\right]\left[p^2+(1p)^2\right].$$
(3)
An alternative way to characterize the difference between two networks is the difference automaton, defined by

$$d_i(t)=\sigma _i^1(t)\sigma _i^2(t),$$
(4)
where $``$ denotes Boolean addition. The density of this automaton is given by

$$D(t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}d_i(t),$$
(5)
and coincides with the Hamming distance between the networks. Note that $`D(t)=1a(t)`$ so that, from Eq. (3), we get

$$D(t+1)=2p(1p)\left(1[1D(t)]^K\right).$$
(6)
The Hamming distance has proven to be a suitable order parameter in the study of the synchronization transition in coupled elementary cellular automata , where the analysis of overlaps between states is a basic tool to define the effects of coupling. In the next section, we adapt the annealed model to the description of coupled KNs.

## III Coupled Kauffman networks

### A Stochastic coupling

In order to establish a coupling mechanism between two KNs, we first observe that, due to the discrete nature of KNs, the usual deterministic coupling used for maps  cannot be applied here. Consequently, we introduce a form of stochastic coupling between networks as previously done for cellular automata , where the continuous parameter $`q`$ that controls the strength of the coupling is a probability, as explained in the following.

The evolution of the coupled system is implemented by the successive application of two operators. First, the evolution function $`𝐟`$ is applied to both networks as if they were not coupled \[see Eq. (2)\], yielding $`𝐟\left(𝝈^1\right)`$ and $`𝐟\left(𝝈^2\right)`$. Next, the stochastic coupling operator $`𝒮`$ is applied:

$$\{𝝈^1(t+1),𝝈^2(t+1)\}=𝒮(𝐟\left[𝝈^1(t)\right],𝐟\left[𝝈^2(t)\right]).$$
(7)
The operator $`𝒮`$ compares the states of the networks site by site. If $`\sigma _i^1(t)=\sigma _i^2(t)`$ the state of the site is not modified. If, on the other hand, $`\sigma _i^1(t)\sigma _i^2(t)`$, with probability $`q`$ the states of the homologous sites in both networks are set to the same value. This value is chosen among $`\sigma _i^1(t)`$ and $`\sigma _i^2(t)`$, with the same probability $`1/2`$ for each instance. With probability $`1q`$, even if $`\sigma _i^1(t)\sigma _i^2(t)`$, the coupling does not act, leaving the state of that site unchanged in both networks. We call $`q`$ the coupling probability. The whole evolution can be formally expressed as

$$\{\sigma _i^1,\sigma _i^2\}\{\begin{array}{cc}\{f_i\left(𝝂_i^1\right),f_i\left(𝝂_i^2\right)\}\hfill & \text{with probability }1q,\hfill \\ \{f_i\left(𝝂_i^1\right),f_i\left(𝝂_i^1\right)\}\hfill & \text{with probability }q/2,\hfill \\ \{f_i\left(𝝂_i^2\right),f_i\left(𝝂_i^2\right)\}\hfill & \text{with probability }q/2.\hfill \end{array}$$
(8)
We stress that we are dealing with two identical extended systems, and that the coupling mechanism connects homologous elements of these two systems, namely, the $`i`$th site of network 1 is connected by coupling with the $`i`$th site of network 2, as schematically illustrated in Fig. 1. The coupling mechanism defined above is symmetric, since each network may influence the other. It could also be possible to consider a biased, non symmetric coupling, in which one network drives the other .

For $`q=0`$ the two networks are uncoupled, and evolve independently from each other. For $`q=1`$, in contrast, the networks synchronize completely at the first time step. From then on, they follow a common trajectory in phase space without further intervention of the coupling mechanism. Our aim in the following is to study the behavior of coupled KNs for intermediate values of the coupling probability, $`q(0,1)`$.

In the frozen phase, were no damage spreading takes place, an arbitrary small coupling intensity $`q>0`$ leads eventually to synchronization. The situation is different in the chaotic phase. There, we find two competing driving forces acting on the coupled system, namely, the chaotic dynamics which induces the separation between two trajectories to grow , and the coupling, by which the Hamming distance between the networks decreases. In this paper, we focus attention on the chaotic phase.

### B Annealed model for coupled networks

The annealed model can be used to predict the behavior of the pair of coupled KNs. We recall that the time evolution equation for the Hamming distance in the case of two free networks is given by Eq. (6). Now suppose that the Hamming distance of two coupled KNs at time $`t`$ is $`D(t)`$. The first substep in the dynamics of this system consists of the free evolution of both networks. The Hamming distance after this substep, $`D(t+\delta t)`$, is therefore given by Eq. (6). At the second substep, coupling acts on the system, and a fraction $`q`$ of the homologous sites that were in different states are assigned the same value. This leaves a fraction $`(1q)D(t+\delta t)`$ of sites with different states in the two networks. Thus, the evolution of the Hamming distance for coupled networks is given by the map

$$D(t+1)=F\left[D(t)\right]=\phi (p,q)\left(1\left[1D(t)\right]^K\right),$$
(9)
with $`\phi (p,q)=2p(1p)(1q)`$.

It can be shown that the map (9) has a stable fixed point $`D^{}>0`$ for $`q<q_c`$, with

$$q_c=1\left[2p(1p)K\right]^1.$$
(10)
At $`q=q_c`$ the system undergoes a transcritical bifurcation , and $`D^{}=0`$ becomes a stable fixed point for $`q>q_c`$. Thus, within the AM approximation, $`q_c`$ stands for the critical coupling at which synchronization sets on. In this paper we deal mostly with the case $`K=3`$, for which the stable equilibrium $`D^{}`$ can be given analytically as a function of the coupling probability $`q`$. In this case, in fact, the map is defined by the cubic function $`F(x)=\phi (p,q)x\left(33x+x^2\right)`$. The stable Hamming distance is

$$D^{}(q)=\{\begin{array}{cc}\frac{3}{2}\frac{1}{2}\left[3+\frac{4}{\phi (p,q)}\right]^{1/2}\hfill & \text{for }q<q_c,\hfill \\ & \\ 0\hfill & \text{for }qq_c.\hfill \end{array}$$
(11)
Note that near the critical point, $`qq_c`$, this Hamming distance is approximately given by

$$D^{}(q)=6p(1p)|qq_c|.$$
(12)
Therefore, the corresponding critical exponent is equal to unity.

For $`qq_c`$, the Hamming distance approaches $`D^{}(q)`$ exponentially in time. For $`q=q_c`$, on the other hand, Eq. (9) can be approximately written, for $`D(t)0`$, as $`D(t+1)=D(t)(K1)D(t)^2/2`$. This implies a power-law decay for long times, $`D(t)t^1`$. In the next section, we compare these analytical results with those of extensive numerical simulations.

## IV Numerical results

We have performed numerical simulations of pairs of KNs coupled under the scheme presented above. The results reported in this section correspond to the case of $`p=1/2`$ and $`K=3`$. We have recorded the time evolution of the Hamming distance, performing averages over $`r`$ realizations of $`N`$-site networks. The number of realizations is chosen in such a way that, for different values of $`N`$, $`rN10^6`$. In a typical realization, we start with two identical networks with different random initial conditions. For each realization, new connections and local functions are chosen. The networks are allowed to evolve freely, without coupling, for a transient time $`\tau `$ of, typically, $`10^3`$ steps. This is done for the networks to reach their asymptotic dynamics before coupling is allowed to act. After this transient, we turn the coupling on, reset the time to zero, and start measuring $`D(t)`$.

In Fig. 2 we show the time evolution of the Hamming distance $`D(t)`$ for different values of the coupling parameter $`q`$. The values of $`q`$ have been chosen to display the three typical behaviors, namely, synchronization for $`q>q_c`$, critical decay for $`qq_c`$, and convergence to a finite distance for $`q<q_c`$.

For the present values of $`p`$ and $`K`$, the annealed model predicts a critical coupling probability $`q_c=1/3`$ \[cf. Eq. (10)\]. It is however clear from Fig. 2 that the power-law decay in $`D(t)`$ is observed for a lower coupling, $`q=0.29`$. Simulations of the same networks with $`q=1/3`$, on the other hand, always lead to synchronization. In fact, the annealed model is expected to provide a good approximation to our system in the limit $`N\mathrm{}`$ . Figure 3 shows the time evolution of the Hamming distance for a fixed coupling probability, $`q=0.29`$, and several values of $`N`$. The AM prediction from Eq. (9) is also shown. The strong dependence with the size of the network is apparent. In particular, we find that for this coupling strength $`10^3`$-site networks synchronize whereas $`10^4`$-site networks do not. The AM result gives a good description for the case of $`N=10^4`$.

The dependence of $`D(t)`$ on the size of the system can be partially ascribed to the effect of spurious synchronization. Because of the discrete nature of KNs, two finite-size networks can be brought to the same state by a fluctuation caused by the stochastic coupling. In such case, the two networks will remain synchronized from then on. This event is more frequent for small networks, where the relative size of fluctuations increases. Near the critical point, furthermore, where the Hamming distance vanishes asymptotically, the effect of fluctuations is strongly enhanced. The net result of spurious synchronization is that the effective critical coupling for finite-size networks shifts to lower values as $`N`$ decreases. As a consequence, the average Hamming distance in our simulations may vanish even for coupling probabilities below $`q_c`$, as illustrated in Figs. 2 and 3.

Spurious synchronization can be avoided by adding noise to the system. Such strategy has already been adopted in this field, specifically, in the study of globally coupled chaotic maps , to prevent synchronization due to round-off errors in computer simulations. We implement the addition of noise as a new substep in the dynamics of our system. After the evolution and coupling substeps, we flip the state of each site in one of the networks with a small probability $`\eta `$. Figure 4 illustrates the effect of noise in the evolution of $`D(t)`$ for $`q=0.29`$ and $`N=10^3`$. For this coupling intensity, where the consequences of spurious synchronization are vast, the behavior of the Hamming distance with and without noise changes drastically.

Note that noise eliminates spurious synchronization for $`q<q_c`$, but also prevents the KNs to exactly synchronize even for $`q>q_c`$. Therefore, the critical behavior that characterizes the synchronization transition in the absence of noise disappears as noise is added, and is recovered only for $`\eta 0`$ (but $`\eta 0`$). The effect of noise can be straightforwardly incorporated to the AM approximation.

The map that gives the time evolution of the Hamming distance is now

$$D(t+1)=(1\eta )F\left[D(t)\right]+\eta \left(1F\left[D(t)\right]\right),$$
(13)
with $`F\left[D(t)\right]`$ given by Eq. (9). As for the model without noise, for $`K=3`$ it is possible to analytically find the asymptotic distance $`D_\eta ^{}`$ predicted by Eq. (13). Figure 5 shows $`D_\eta ^{}`$ as a function of the coupling probability $`q`$ for various noise intensities. Note the approximation to the critical behavior as $`\eta 0`$. The insert displays the difference between $`D_\eta ^{}`$ and the asymptotic distance in the absence of noise as a function of $`q`$.

In Fig. 6, we compare the prediction of Eq. (13) with numerical results for the Hamming distance of two coupled $`10^4`$-site KNs with noise intensity $`\eta =10^4`$, for three values of the coupling intensity. The agreement is excellent during the transients, but some noticeable discrepancies persist in the asymptotic value, especially, for $`qq_c`$.

To study such discrepancies in detail, and thus test the AM results, we compute from our numerical simulations the average asymptotic value of $`D(t)`$, defined as

$$D=\frac{1}{T}\underset{t=t_0}{\overset{t_0+T}{}}D(t).$$
(14)
Averages are performed over a time span $`T`$ of $`5\times 10^3`$ time steps, when the asymptotic regime of the coupled system has been reached, i.e. for sufficiently large values of $`t_0`$. As above, we choose $`\eta =10^4`$, and determine $`D`$ for several values of $`N`$. Results for $`N=2^7`$, $`2^{10}`$, and $`2^{14}`$ are presented in Fig. 7. The AM prediction for this value of $`\eta `$ is also shown, as a curve. We see that the AM results systematically overestimate the values of $`D`$, and that the agreement improves for larger values of $`N`$. Therefore, even when noise has been added to avoid spurious synchronization, finite-size effects persist.

These remanent finite-size effects are measured by the difference $`\delta D^{}=D_\eta ^{}D`$ between the AM prediction and the numerical average defined in Eq. (14). In Fig. 8 we plot $`\delta D^{}`$ as a function of $`N`$ for different coupling intensities. The insert shows $`D_\eta ^{}`$ and $`D`$ as a function of $`N`$ for the same values of $`q`$. There are two well-differentiated regimes in the size dependence of $`\delta D^{}`$. For small $`N`$, the deviation between the AM estimate and numerical results is practically constant. For large values of $`N`$ the deviation decreases, seemingly as a power-law, $`\delta D^{}N^z`$. Least-square fits for $`N>10^2`$ yield $`z=1.1\pm 0.2`$ for the exponent.

To find an explanation for these finite size effects, it is necessary to go beyond the annealed approximation. In the following, we outline an approach to the calculation of the Hamming distance based on statistical averages over states of the KNs. The density of a state $`𝝈=(\sigma _1,\mathrm{},\sigma _N)`$ of a single KN is defined as

$$\rho \left(𝝈\right)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\sigma _i$$
(15)
\[cf. Eq. (5)\], whereas the distance between two states $`𝝈`$ and $`𝝈^{}`$ reads

$$D(𝝈,𝝈^{})=\frac{1}{N}\underset{i=1}{\overset{N}{}}(\sigma _i\sigma _i^{})^2.$$
(16)
For fixed $`p`$, $`K`$ and $`N`$, a realization $`R`$ of the network is defined by the connections and the local rules. We define $`\mathrm{\Omega }_R`$ as the set of all the states visited by the KN for this realization, at asymptotically large times and from all the possible initial conditions. In other words, the set $`\mathrm{\Omega }_R`$ contains all the states that belong to the limit cycles of the dynamics. It is possible to introduce a probability distribution $`𝒫_R\left(𝝈\right)`$ over $`\mathrm{\Omega }_R`$, given by the frequency with which a given state $`𝝈`$ is visited at asymptotically large times averaged over all initial conditions. Averages over $`\mathrm{\Omega }_R`$ will be computed with this distribution. For instance, the average of the density $`\rho `$ is

$$d_R=\rho (𝝈)_{\mathrm{\Omega }_R}=\underset{𝝈\mathrm{\Omega }_R}{}𝒫_R(𝝈)\rho (𝝈).$$
(17)
The average of the distance between two states, Eq. (16), can be written as

$$D(𝝈,𝝈^{})_{\mathrm{\Omega }_R}=2d_R(1d_R)\frac{2}{N}\underset{i=1}{\overset{N}{}}\xi _i^2,$$
(18)
where $`\xi _i=\sigma _i_{\mathrm{\Omega }_R}d_R=\sigma _i^{}_{\mathrm{\Omega }_R}d_R`$. Here, we have assumed that the occurrence of states $`𝝈`$ and $`𝝈^{}`$ are probabilistically uncorrelated.

The quantities $`\xi _i`$ measure, for each realization of the network, the average deviation of the state of each site from the average density $`d_R`$ over the whole network. Let us introduce the distribution $`\mathrm{\Gamma }_R(\xi )`$ as the fraction of sites in the KN with deviation $`\xi `$ for a specific realization $`R`$. Unfortunately, the explicit form of $`\mathrm{\Gamma }_R(\xi )`$ is not known. It is however known that this is a nontrivial distribution, in particular, due to existence of the so-called stable core . The stable core is a set of sites that have always the same asymptotic, fixed states, irrespectively of the initial condition. For these sites, $`\sigma _i_{\mathrm{\Omega }_R}=0`$ or $`1`$, so that the deviations $`\xi _i`$ adopt their extremal values, $`\xi _i=d_R`$ or $`1d_R`$, respectively. Using the distribution $`\mathrm{\Gamma }_R(\xi )`$ to replace the sum in Eq. (18) by an integral, we have

$$D(𝝈,𝝈^{})_{\mathrm{\Omega }_R}=2d_R(1d_R)2_{d_R}^{1d_R}\mathrm{\Gamma }_R(\xi )\xi ^2𝑑\xi .$$
(19)
Now, analytical results on the size of the stable core  suggest that, as $`N\mathrm{}`$, $`\mathrm{\Gamma }_R(\xi )`$ approaches an asymptotic profile $`\mathrm{\Gamma }(\xi )`$ which depends on $`p`$ and $`K`$, but becomes independent of the specific realization of the network. For large sizes, we may assume a dependence of the form

$$\mathrm{\Gamma }_R(\xi )\mathrm{\Gamma }(\xi )N^1\mathrm{\Gamma }_R^{}(\xi ),$$
(20)
where $`\mathrm{\Gamma }_R^{}`$ is the first (analytical) correction due to finite sizes. Within these assumptions, Eq. (19) takes the form

$$D(𝝈,𝝈^{})_{\mathrm{\Omega }_R}=D_0\frac{2}{N}_{d_R}^{1d_R}\mathrm{\Gamma }_R^{}(\xi )\xi ^2𝑑\xi ,$$
(21)
where

$$D_0=2d_R(1d_R)2_{d_R}^{1d_R}\mathrm{\Gamma }(\xi )\xi ^2𝑑\xi $$
(22)
is the asymptotic distance for $`N\mathrm{}`$.

We now associate $`D(𝝈,𝝈^{})_{\mathrm{\Omega }_R}`$ with the distance between the states of two coupled KNs at a given (long) time. Indeed, in our system both networks have the same connections and rules, and correspond therefore to the same realization $`R`$ of the network. According to the definition (4) and (5), after an average over realizations of the network for fixed $`p`$, $`K`$, and $`N`$ is performed, the distance $`D(𝝈,𝝈^{})_{\mathrm{\Omega }_R}`$ coincides with the Hamming distance $`D`$, Eq. (14), considered above. In Eq. (21), thus, $`D_0`$ should correspond to the Hamming distance predicted by the AM approximation, valid for $`N\mathrm{}`$, and the first correction to the AM estimate is given by the additional term. As found from our simulations, this term depends on the network size as $`N^1`$. The ansatz (20) is therefore supported by numerical results. Note that these conclusions are independent of the strength of coupling, measured by the probability $`q`$, since the only effect of the coupling dynamics in connection with the above analysis is to change the set of asymptotic states $`\mathrm{\Omega }_R`$ and the profile of the distribution $`\mathrm{\Gamma }_R(\xi )`$.

## V Disordered cellular automata and partially ordered networks

According to the results reported in the previous sections, the synchronization transition of stochastically coupled KNs on the one hand, and of stochastically coupled elementary cellular automata (ECA, ) on the other, are qualitatively different. Namely, they belong to different universality classes. Whereas we have found that synchronization in KNs appears through a transcritical bifurcation—with a critical exponent equal to unity for the asymptotic Hamming distance—the corresponding transition in ECA has been shown to exhibit nontrivial exponents , which seem to suggest that it belongs to the universality class of directed percolation . It is therefore relevant to study a third class of systems, intermediate between ECA and generic KNs.

In a generic KN there are two sources of disorder. We have the network topology, determined by the random choice of connections, and the local evolution rules, which are also chosen at random. On the other hand, in cellular automata both the topology and the dynamical rules are fully homogeneous. ECA can indeed be interpreted as a very special subclass of KNs with $`K=3`$, where the choice of dynamical rules and connections is deterministic. In order to distinguish between the effects of disorder in the topology and in the dynamics on the synchronization transition, we consider now the subclass of KNs where the connections are still chosen at random but the local rule is the same for all sites. We refer to these networks as disordered cellular automata (DCA).

We focus here on the elementary rule $`22`$, which exhibits chaotic evolution . This rule is defined by the Boolean function $`f(\{0,0,1\})=f(\{0,1,0\})=f(\{1,0,0\})=1`$ and $`f=0`$ for the remaining five possible neighborhoods (see Table I). Thus, the bias for rule $`22`$ is $`p=3/8`$. This DCA, however, cannot be thought of as a KN with $`p=3/8`$. In a generic KN with this bias, in fact, the local evolution could be substantially different from the behavior of rule $`22`$. In particular, the dynamics at some sites may be governed by nonchaotic rules, giving rise to sensible differences in the global behavior. This is clearly illustrated, for instance, by a measurement of the asymptotic density of a rule-$`22`$ DCA, which yields $`d0.423`$, instead of $`d=p=0.375`$.

The formulation of an annealed model of DCA requires a different approach, in order to account for the homogeneity of the dynamical rules. We define the AM by reassigning all the connections at each time step, but keeping the functions $`f_i=f`$ fixed. Suppose that we have two networks with overlap $`a(t)`$. The probability for a site to have exactly $`Kl`$ of its reassigned inputs coming from homologous sites in the same state is

$$p_l(t)=\left(\genfrac{}{}{0pt}{}{K}{l}\right)a(t)^{Kl}\left[1a(t)\right]^l.$$
(23)
We introduce the quantity $`A_l`$ as the probability for two homologous sites having all but $`l`$ inputs coming from homologous sites in the same state to give the same output. The overlap at time $`t+1`$ is then given by

$$a(t+1)=\underset{l=0}{\overset{K}{}}A_lp_l(t).$$
(24)
The quantities $`A_l`$ depend on the evolution rule, and their value is fixed. They can be evaluated within some approximations, as shown in the Appendix. Note that $`A_0=1`$ because, no matter the rule, if the inputs are all equal the outputs will coincide. In the annealed model for KNs we had $`A_0=1`$ and $`A_l=p^2+(1p)^2`$ for $`l=1,\mathrm{},K`$. For $`K=3`$ the map for the Hamming distance can be casted into the form

$$D(t+1)=B_1D(t)+B_2D(t)^2+B_3D(t)^3,$$
(25)
with $`B_1=3\left(1A_1\right)`$, $`B_2=3\left(2A_1A_21\right)`$, $`B_3=3A_1+3A_2A_3+1`$. Coupling enters then the formulation exactly as in Eq. (9), as an additional factor $`1q`$ in the evolution of $`D(t)`$.

In Fig. 9, the curve stands for the asymptotic value $`D^{}`$ as a function of the coupling probability predicted from Eq. (25) for rule $`22`$. The prediction is qualitatively similar to that for KNs, in particular, in the region close to the synchronization transition. Numerical results on DCA with rule $`22`$ for $`N=2^{11}`$ are also shown in Fig. 9. To avoid spurious synchronization, a small amount of noise, $`\eta =10^4`$, has been added. The agreement with the AM is reasonably good, though some systematic deviations are clearly visible in the zone of the transition. As before, these deviations may be attributed to finite-size effects.

Finally, as an interpolation between DCA and ECA we consider partially ordered topologies (POCA), constructed in the following way. We start with an $`N`$-site ECA, which consists of a one-dimensional array where each site is connected to itself and to its two nearest neighbors. Then, for each site, each neighbor is redrawn at random from the whole network with probability $`\theta `$. With the complementary probability, $`1\theta `$, the neighbor is preserved. The resulting topological structure of connections is analogous to that of small-world networks . For $`\theta =1`$, we recover the DCA networks discussed above. Results of numerical simulations on POCA with rule $`22`$ are shown in Fig. 9 for $`\theta =0.1`$ and $`\theta =0.3`$. Data for ECA with rule $`22`$ ($`\theta =0`$) are also shown. We see that, in spite of the relatively low values of $`\theta `$, the asymptotic Hamming distance of POCA depends on the coupling probability $`q`$ in a way qualitatively similar to that of KNs. In particular, it does not exhibit the abrupt dependence on $`q`$ observed in ECA near the critical point. This would be in agreement with the crossover scenario in small-world networks , where it is known that even small amounts of disorder induce behaviors which already resemble that of fully random structures.

## VI Summary and discussion

We have studied the behavior of two Kauffman networks interacting through a form of symmetric stochastic coupling. As for many other localized or extended, random or deterministic, coupled dynamical systems , we have found that coupled Kauffman networks can synchronize their evolution if coupling is strong enough. In our case, there is a critical value of the coupling probability $`q`$ beyond which the two networks converge to the same trajectory as time elapses.

In contrast with the situation encountered for other extended systems , however, for Kauffman networks it has been possible to give an analytical description of the synchronization transition, in excellent agreement with numerical results for large-size systems. This formulation is provided by an extension of the so-called annealed model  to the system of coupled networks. The model gives the evolution of the overlap between the two networks—or, in other words, of their Hamming distance—and makes it possible to evaluate its asymptotic value. The asymptotic Hamming distance $`D^{}`$ is used as an order parameter for the synchronization transition. For coupled Kauffman networks in the chaotic phase, the annealed model predicts the existence of a critical coupling probability $`q_c`$, such that $`D^{}`$ is finite for $`q<q_c`$ and vanishes for $`q>q_c`$. At the critical point, $`D^{}|qq_c|`$, and the transition has the character of a transcritical bifurcation.

The transition predicted by the annealed model is qualitatively different from that observed in the deterministic version of this kind of systems, more specifically, in cellular automata. For elementary cellular automata, in fact, the critical behavior of the Hamming distance exhibits nontrivial exponents . It has been suggested that, at least for some evolution rules, the synchronization transition in cellular automata belongs to the class of directed percolation . On the other hand, by analogy with the problem of damage spreading, synchronization in Kauffman networks could belong to the class of directed percolation in disordered systems .

Results from extensive numerical simulations of relatively large networks ($`10^4`$ sites), performed for Kauffman networks in the chaotic regime ($`p=1/2`$, $`K=3`$), are in good agreement with the predictions of the annealed model both in the temporal evolution and in the asymptotic behavior of the system. However, for small networks some systematic departures from the analytical results are apparent. These deviations can be partially explained taking into account the occasional events of spurious synchronization for $`q<q_c`$, due to the effect of suitably large fluctuations on our discrete finite-size system. The effect of these fluctuations becomes more important as the coupling probability approaches its critical value. In any case, spurious synchronization can be successfully eliminated by adding a small amount of noise to the dynamics. Moreover, noise can also be encompassed into the annealed model. We have introduced noise in the simulations and observed that remanent finite-size effects persist. These must now be ascribed to the failure of the annealed model in describing finite networks. By means of a more detailed statistical description of the overlaps between networks, in fact, we have been able to account for such remanent effects, also explaining the dependence of the deviation from the annealed model on the network size.

Finally, we have presented preliminary numerical results on disordered and partially ordered stochastically coupled cellular automata. These systems can be seen as providing a kind of interpolation between Kauffman networks, with their completely random connections and dynamical rules, and cellular automata, which are fully ordered. In the case of disordered networks with the same evolution rule on all the sites, it is possible to extend the annealed model, which predicts the same class of synchronization transition as for Kauffman networks. Numerical simulations agree with these predictions. Increasing the order in the connections by means of a scheme analogous to the construction of a small-world network , we have also considered partially disordered structures. Even for small amounts of structural disorder, the behavior associated with the synchronization transition seems to resemble that of Kauffman networks more than that of cellular automata. This leads us to conjecture that the synchronization transition of partially disordered automata is in the same universality class as for coupled Kauffman networks. Nevertheless, further extensive simulations and a careful determination of the critical exponents are necessary to draw any conclusions on this point.

## Acknowledgment

This work has been partially carried out at the Abdus Salam Centre for Theoretical Physics (Trieste, Italy). The authors thank the Centre for hospitality.

## Coefficients for the annealed model

In this Appendix we illustrate the computation of the coefficients $`A_l`$ in Eq. (24) with an explicit example. We recall that $`A_l`$ is defined as the probability for two homologous sites having all but $`l`$ inputs coming from homologous sites in the same state to give the same output.

Let us assume that the frequency with which a given neighborhood appears in a state of the whole network is fully determined by the density $`d`$ of the state. Namely, we neglect the correlations between neighborhood frequencies, associated with the spatial patterns generated by the dynamics. For $`K=3`$, there are eight possible neighborhoods, which we label from $`0`$ to $`7`$ as shown in Table I. Within the above assumption, the frequency $`p_i`$ of each neighborhood $`i`$ can be estimated in terms of the density $`d`$ as

$$\begin{array}{c}p_0=(1d)^3,\hfill \\ p_1=p_2=p_4=d(1d)^2,\hfill \\ p_3=p_5=p_6=d^2(1d),\hfill \\ p_7=d^3.\hfill \end{array}$$
(26)
Moreover, we note that the density $`d`$ can in turn be written in terms of the frequencies $`p_i`$ of neighborhoods with nonzero output. For rule $`22`$ these neighborhoods are $`\{1,0,0\}`$, $`\{0,1,0\}`$ and $`\{0,0,1\}`$ (see Table I), so that we have

$$d=p_1+p_2+p_4.$$
(27)
Combining Eqs. (26) and (27) yields $`d=1\sqrt{3}/30.423`$, which agrees with the numerical measurement reported in Sect. V. The corresponding values of the frequencies $`p_i`$ are shown in Table I. They are in very good agreement with the values obtained from numerical simulations, also shown in the table. This supports our above assumption of uncorrelated neighborhood frequencies.

Once the frequencies $`p_i`$ are known, we are able to calculate the coefficients $`A_l`$. As a specific example, we discuss $`A_3`$. By definition, this is the probability for two homologous sites having all but $`3`$ inputs coming from homologous sites in the same state to give the same output. Thus, we are in the case where the two neighborhoods have all the homologous sites in different states. The pairs of neighborhoods that satisfy this condition are $`\{0,7\}`$, $`\{1,6\}`$, $`\{2,5\}`$, and $`\{3,4\}`$. Among them, however, only the pair $`\{0,7\}`$ has the same output for both neighborhoods (see Table I). The coefficient $`A_3`$ is therefore given by

$$A_3=\frac{p_0p_7}{p_0p_7+p_1p_6+p_2p_5+p_3p_4}=\frac{1}{4}.$$
(28)
The computation of the other coefficients is accomplished in a similar way, yielding $`A_1=4(198\sqrt{3})/169`$ and $`A_2=(9+\sqrt{3})/13`$. Moreover, $`A_0=1`$.


text:
# Recent structure function results from neutrino scattering at Fermilab To be published in proceedings of the XXXth INternational Conferencee on High Energy Physics, Osaka, Japan, July, 2000

## Abstract

We report on the extraction of the structure functions $`F_2`$ and $`\mathrm{\Delta }xF_3=xF_3^\nu xF_3^{\overline{\nu }}`$ from CCFR $`\nu _\mu `$-Fe and $`\overline{\nu }_\mu `$-Fe differential cross sections. The extraction is performed in a physics model independent (PMI) way. This first measurement of $`\mathrm{\Delta }xF_3`$, which is useful in testing models of heavy charm production, is higher than current theoretical predictions. The ratio of the $`F_2`$ (PMI) values measured in $`\nu _\mu `$ and $`\mu `$ scattering is in agreement (within 5%) with the NLO predictions using massive charm production schemes, thus resolving the long-standing discrepancy between the two sets of data. In addition, measurements of $`F_L`$ (or, equivalently, $`R`$) and $`2xF_1`$ are reported in the kinematic region where anomalous nuclear effects in $`R`$ are observed at HERMES. \[Preprint UR-1614, ER/40685/952\]

Deep inelastic lepton-nucleon scattering experiments have been used to determine the quark distributions in the nucleon. However, the quark distributions determined from $`\mu `$ and $`\nu `$ experiments were found to be different at small values of $`x`$, because of a disagreement in the extracted structure functions. Here, we find that the neutrino-muon difference is resolved by extracting the $`\nu _\mu `$ structure functions from CCFR neutrino data in a physics model independent (PMI) way. In addition, measurements of $`\mathrm{\Delta }xF_3`$, $`F_L`$, and $`2xF_1`$ are presented.

The sum of $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ differential cross sections for charged current interactions on an isoscalar target is related to the structure functions as follows:

| $`F(ϵ)`$ | $`\left[\frac{d^2\sigma ^\nu }{dxdy}+\frac{d^2\sigma ^{\overline{\nu }}}{dxdy}\right]\frac{(1ϵ)\pi }{y^2G_F^2ME_\nu }`$ |
| --- | --- |
|  | $`=2xF_1[1+ϵR]+\frac{y(1y/2)}{1+(1y)^2}\mathrm{\Delta }xF_3.`$ (1) |

Here $`G_F`$ is the Fermi weak coupling constant, $`M`$ is the nucleon mass, $`E_\nu `$ is the incident energy, the scaling variable $`y=E_h/E_\nu `$ is the fractional energy transferred to the hadronic vertex, $`E_h`$ is the final state hadronic energy, and $`ϵ2(1y)/(1+(1y)^2)`$ is the polarization of the virtual $`W`$ boson. The structure function $`2xF_1`$ is expressed in terms of $`F_2`$ by $`2xF_1(x,Q^2)=F_2(x,Q^2)\times \frac{1+4M^2x^2/Q^2}{1+R(x,Q^2)}`$, where $`Q^2`$ is the square of the four-momentum transfer to the nucleon, $`x=Q^2/2ME_h`$ is the fractional momentum carried by the struck quark, and $`R=\frac{\sigma _L}{\sigma _T}`$ is the ratio of the cross-sections of longitudinally- to transversely-polarized $`W`$ bosons. The $`\mathrm{\Delta }xF_3`$ term, which in leading order $`4x(sc)`$, is not present in the $`\mu `$-scattering case. In addition, there is a threshold suppression originating from the production of heavy $`c`$ quarks in a $`\nu _\mu `$ charged current interaction with $`s`$ quarks. For $`\mu `$-scattering, there is no suppression for scattering from $`s`$ quarks, but more suppression when scattering from $`c`$ quarks.

In previous analyses of $`\nu _\mu `$ data, structure functions were extracted by applying a slow rescaling correction to correct for the charm mass suppression in the final state. In addition, the $`\mathrm{\Delta }xF_3`$ term from a leading order charm production model was used as input in the extraction. These resulted in physics model dependent (PMD) structure functions. In the new analysis reported here, slow rescaling corrections are not applied. $`\mathrm{\Delta }xF_3`$ and $`F_2`$ are extracted from two-parameter fits to the $`F(ϵ)`$ distributions according to Eq. (1). However, in the $`x>0.1`$ region, we extract values of $`F_2`$ with $`\mathrm{\Delta }xF_3`$ constrained to the NLO TR-VFS(MRST) predictions. Since $`\mathrm{\Delta }xF_3`$ for $`x>0.1`$ is small, the extracted values of $`F_2`$ are insensitive to $`\mathrm{\Delta }xF_3`$.

Fig. 1(left) shows the extracted values of $`\mathrm{\Delta }xF_3`$ as a function of $`x`$ (above $`Q^2=1`$ ), including both statistical and systematic errors, compared to various theoretical methods for modeling heavy charm productions within a QCD framework. Fig. 1(right) shows the sensitivity to the choice of scale. With reasonable choices of scale, all the theoretical models yield similar results. However, at low $`Q^2`$, our $`\mathrm{\Delta }xF_3`$ data are higher than all theoretical models.

Our $`F_2`$ (PMI) measurements divided by the NLO TR-VFS(MRST) predictions are shown in Fig. 2(left). Also shown are $`F_2^\mu `$ and $`F_2^e`$ divided by the theory predictions. Nuclear effects, target mass, and higher twist corrections are included in the calculation. As shown in Fig. 2, within 5% both the neutrino and muon structure functions are in agreement with the NLO TR-VFS(MRST) predictions, and therefore in agreement with each other, thus resolving the long-standing discrepancy between the two sets of data. A comparison using the NLO ACOT-VFS(CTEQ4HQ) predictions yields similar results. Note that previously there was up to a 20% difference between the CCFR $`F_2`$ (PMD) and NMC data at $`x=0.015`$, as shown in Fig. 2(right).

Recently, there has been a renewed interest in $`R`$ at small $`x`$ and $`Q^2<1`$, because of the large anomalous nuclear effect that has been reported by the HERMES experiment. Their measurement implies a large enhancement in $`F_L`$ but suppression in $`2xF_1`$ in heavy nuclear targets. It is expected that any nuclear effect in $`R`$ would be enhanced in the CCFR iron target with respect to the nitrogen target in HERMES, unless the origin of this effect depends on the incident probe (electron versus neutrino).

Values of $`F_L`$ and $`2xF_1`$ are extracted from the sums of the corrected $`\nu _\mu `$-Fe and $`\overline{\nu }_\mu `$-Fe differential cross sections in different energy bins according to Eq. (1). An extraction of $`F_L`$ requires knowledge of $`\mathrm{\Delta }xF_3`$. which we obtain from the NLO TR-VFS(MRST) calculation. Because of the large uncertainty in $`\mathrm{\Delta }xF_3`$ at low $`Q^2`$ region, an extrapolation of the curve which describes the measured CCFR $`\mathrm{\Delta }xF_3`$ data above $`Q^2=1`$ is used for the systematic error. Here we are interested in the relative $`Q^2`$ dependence of $`F_L`$ and $`2xF_1`$.

Fig. 3 shows the preliminary values of $`F_L`$ and $`2xF_1`$ as a function of $`Q^2`$ for $`x<0.1`$. The inner errors include both statistical and experimental systematic errors. The outer errors represent the $`\mathrm{\Delta }xF_3`$ model errors added in quadrature. The curves are the predictions from a QCD-inspired leading order fit to the CCFR differential cross section data with $`R`$= $`R_{world}`$ (for neutrino scattering) which does not include the HERMES effect. Large anomalous deviations from the fit (e.g. 200 - 300%) are not seen in the CCFR data.

More details on this work can be found in reference 6 and 7.


text:
# The detection of spiral arm modulation in the mass distribution of an optically flocculent galaxy

## References

Bertin, G., Lin, C.C., Lowe, S.A., Thurstans, R.P., 1989, ApJ 338, 78

Block, D.L., Puerari, I., 1999, A&A 342, 627

Block, D.L., Wainscoat, R.J., 1991, Nature 353, 48

Considère, S., Athanassoula, E., 1982, A&A 111, 28

Elmegreen, B.G., 1995, in Physics of the Interstellar Medium and Intergalactic Medium, eds. A. Ferrara, C.F. McKee, C. Helles and P.R. Shapiro, ASP Conference Series, 80, 218

Elmegreen, B.G., Elmegreen, D.M., Seiden, P.E., 1989, ApJ 343, 602

Lin, C.C., Shu, F.H., 1964, ApJ 140, 646

Puerari, I., 1993, PASP 105, 1290

Puerari, I., Dottori, H.A., 1992, A&AS 93, 469

Puerari, I., Block, D.L., Elmegreen, B.G., Frogel, J.A., Eskridge, P.B., 2000, A&A, 359, 932


text:
# Implications of the ALEPH 𝜏-Lepton Decay Data for Perturbative and Non-Perturbative QCD

\[

## Abstract

We use ALEPH data on hadronic $`\tau `$ decays in order to calculate Euclidean coordinate space correlation functions in the vector and axial-vector channels. The linear combination $`VA`$ receives no perturbative contribution and is quantitatively reproduced by the instanton liquid model. In the case of $`V+A`$ the instanton calculation is in good agreement with the data once perturbative corrections are included. These corrections clearly show the evolution of $`\alpha _s`$. We also analyze the range of validity of the Operator Product Expansion (OPE). In the $`VA`$ channel we find a dimension $`d=6`$ contribution which is comparable to the original SVZ estimate, but the instanton model provides a different non-singular term of the same magnitude. In the $`V+A`$ case both the OPE and the instanton model predict the same $`d=4`$ power correction induced by the gluon condensate, but it is masked by much larger perturbative contributions. We conclude that the range of validity of the OPE is limited to $`x<0.3`$ fm, whereas the instanton model describes the data over the entire range.

\]

1. Quantitative understanding of the interface between perturbative and non-perturbative effects is the central problem in QCD dynamics. Historically, QCD sum rules based on the Operator Product Expansion (OPE)  constituted the first serious attempt to describe non-perturbative phenomena in QCD. The initial application of QCD sum rules to vector and axial vector meson lead to very promising results. It was soon discovered, however, that not all hadrons are alike . Phenomenology demands that non-perturbative effects in scalar and pseudo-scalar channels, both meson and glueball, are much bigger than in the vector channels. This fact is not reproduced by the OPE but it was realized that direct instanton effects appear in exactly those channels in which non-perturbative effects are large. This observation gave rise to the instanton liquid model .

The available information on hadronic correlation functions, both from experimental data, the OPE and other exact results was reviewed in . Since then, the high statistics measurement of hadronic $`\tau `$ decays $`\tau \nu _\tau +\mathrm{hadrons}`$ by the ALEPH experiment at CERN  has significantly improved the experimental situation in the vector and axial-vector channel. The purpose of this paper is to compare these results with theoretical predictions, both from the OPE and instanton models. In particular, we would like to assess the range of applicability of the two approaches and put improved constraints on the parameters that enter into the calculations. Translating the spectral functions measured by the ALEPH collaboration into Euclidean coordinate space correlation functions will also allow precise comparison of the experimental data with improved lattice calculations along the lines of .

2. In the following, we shall consider the vector and axial-vector correlation functions $`\mathrm{\Pi }_V(x)=j_\mu ^a(x)j_\mu ^a(0)`$ and $`\mathrm{\Pi }_A(x)=j_\mu ^{5a}(x)j_\mu ^{5a}(0)`$. Here, $`j_\mu ^a(x)=\overline{q}\gamma _\mu \frac{\tau ^a}{2}q`$, $`j_\mu ^{5a}(x)=\overline{q}\gamma _\mu \gamma _5\frac{\tau ^a}{2}q`$ are the isotriplet vector and axial-vector currents. The correlation functions have the spectral representation

$$\mathrm{\Pi }_{V,A}(x)=𝑑s\rho _{V,A}(s)D(\sqrt{s},x),$$
(1)
where $`D(m,x)=m/(4\pi ^2x)K_1(mx)`$ is the Euclidean coordinate space propagator of a scalar particle with mass $`m`$. We shall focus on the linear combinations $`\mathrm{\Pi }_V+\mathrm{\Pi }_A`$ and $`\mathrm{\Pi }_V\mathrm{\Pi }_A`$. These combinations allow for a clearer separation of different non-perturbative effects. The corresponding spectral functions $`\rho _V\pm \rho _A`$ which are measured by the ALEPH collaboration are shown in Fig. 1.

In QCD, the vector and axial-vector spectral functions have to satisfy chiral sum rules. If we assume that $`\rho _V(s)\rho _A(s)=0`$ above the maximum invariant mass $`s=m_\tau ^2`$ for which the spectral functions can be measured, then we find that the ALEPH data satisfy all chiral sum rules within the experimental uncertainty. However, the central values of the sum rules differ significantly from the chiral predictions . In general, both $`\rho _V`$ and $`\rho _A`$ are expected to show oscillations of decreasing amplitude . If we set $`\rho _V(s)\rho _A(s)=0`$ above an arbitrarily chosen invariant mass $`s_0`$ this will lead to the appearance of spurious dimension $`d=2,4`$ operators in the correlation functions at small $`x`$. For this reason we have decided to use $`s_0=2.5\mathrm{GeV}^2`$, which is slightly below the tau mass but allows all chiral sum rules to be satisfied. The reader should be aware of the fact that we have, in effect, slightly moved the data points in the small $`x`$ region within the error bars reported by the ALEPH collaboration. Finally we add the pion pole contribution, which is not shown in Fig. 1, to the axial vector spectral function. This corresponds to an extra term $`\mathrm{\Pi }_A^\pi (x)=f_\pi ^2m_\pi ^2D(m_\pi ,x)`$. The resulting correlation functions $`\mathrm{\Pi }_V(x)\pm \mathrm{\Pi }_A(x)`$ are shown in Fig. 2.

3. We begin our analysis with the combination $`\mathrm{\Pi }_V\mathrm{\Pi }_A`$. This combination is sensitive to chiral symmetry breaking, while perturbative diagrams, as well as gluonic operators cancel.

In Fig. 2 we compare the measured correlation functions with predictions from the instanton model. These predictions are described in great detail in  and the review . The main assumption is that the QCD vacuum is dominated by strong non-perturbative field configurations, instantons. In the simplest model, the random instanton liquid (RILM), the instanton positions and color orientations are distributed randomly. The ensemble is characterized by two numbers, the instanton-anti-instanton density $`(N/V)=1\mathrm{fm}^4`$ and the average instanton size $`\rho =1/3`$ fm. These parameters were fixed a long time ago using the requirement that they must reproduce the phenomenological values of the quark and gluon condensates . The agreement of the instanton prediction with the measured $`VA`$ correlation is impressive and extends all the way from short to large distances. At distances $`x>1.25`$ fm both combinations are dominated by the pion contribution while at intermediate $`x`$ the $`\rho ,\rho ^{}`$ and $`a_1`$ resonances contribute.

In order to study the validity of the operator product expansion we have to study the short distance region in more detail. The OPE predicts that the $`VA`$ correlation function starts with the the following quark-anti-quark operators of dimension $`d=4`$ and $`d=6`$

$`{\displaystyle \frac{\mathrm{\Pi }_V(x)\mathrm{\Pi }_A(x)}{2\mathrm{\Pi }_0(x)}}`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{4}}m\overline{q}qx^4`$ (3)
$`+{\displaystyle \frac{\pi ^3}{9}}\alpha _s(x)\overline{q}q^2\mathrm{log}(x^2)x^6+\mathrm{}.`$
The value of the dimension $`d=4`$ operator is determined by the Gell-Mann-Oakes-Renner (GMOR) relation to be $`(x/1.66\mathrm{fm})^4`$. Using $`\overline{q}q=(230\mathrm{MeV})^3`$ and the one-loop running coupling constant we also estimate the $`d=6`$ operator as $`(x/0.66\mathrm{fm})^6`$. This implies that the $`d=6`$ operator totally dominates over the $`d=4`$ operator.

This estimate can be checked by considering the value of the $`d=6`$ operator as a free parameter and trying to extract it from the measured data. Because the $`d=4`$ operator is so small we use the GMOR value. A similar determination of power corrections was already done by the ALEPH collaboration using moments sum rules . Nevertheless, fitting the OPE coefficients in coordinate space provides important additional insight. The results depend on the coordinate range $`[0,x_m]`$ used in the fit, but for $`x_m<0.3`$ fm this dependence is weak. The result for $`x_m=0.3`$ fm, also shown in Figs. 2,3, is

$$\frac{\mathrm{\Pi }_V(x)\mathrm{\Pi }_A(x)}{2\mathrm{\Pi }_0(x)}=\left(\frac{x}{1.66\mathrm{fm}}\right)^4+\left(\frac{x}{0.66\mathrm{fm}}\right)^6+\mathrm{}.$$
(4)
We find that the size of the dimension $`d=6`$ term agrees with the SVZ prediction. However, the range of convergence of the OPE is only $`x<0.3`$ fm. We also note that the accuracy of the data in this regime is very poor.

We can also check the short distance behavior of the correlation function in the instanton liquid. Instantons generate the same $`d=4`$ operator in the OPE but the nature of the $`d=6`$ operator is different. To leading order in the semi-classical expansion there is no radiatively generated $`\alpha _s\overline{q}q^2\mathrm{log}(x^2)x^6`$ operator, but instead there is a non-singular $`\overline{q}q^2x^6`$ term. Such terms are dropped in standard OPE, but they are present in the correlation functions. The numerical value of this term is $`(x/0.64\mathrm{fm})^2`$, close to the data and the OPE term.

4. We shall now focus our attention on the $`V+A`$ correlation function. The unique feature of this function is that the full correlator is close to the free field result for distances as large as 1 fm. This phenomenon was referred to as “super-duality” in .

The instanton model reproduces this feature of the $`V+A`$ correlator. We also notice that for small $`x`$ the deviation of the correlator in the instanton model from free field behavior is small compared to the perturbative $`O(\alpha _s/\pi )`$ correction. This opens the possibility of precision studies of the pQCD contribution. But before we do so, let us compare the correlation functions to the OPE prediction

$`{\displaystyle \frac{\mathrm{\Pi }_V(x)+\mathrm{\Pi }_A(x)}{2\mathrm{\Pi }_0(x)}}`$ $`=`$ $`1+{\displaystyle \frac{\alpha _s}{\pi }}{\displaystyle \frac{1}{384}}g^2(G_{\mu \nu }^a)^2x^4`$ (6)
$`{\displaystyle \frac{2\pi ^3}{81}}\alpha _s(x)\overline{q}q\mathrm{log}(x^2)x^6+\mathrm{}`$
Note that the perturbative correction is attractive, while the power corrections of dimension $`d=4`$ and $`d=6`$ are repulsive. Direct instantons also induce an $`O(x^4)`$ correction $`1\frac{\pi ^2}{12}\left(\frac{N}{V}\right)x^4+\mathrm{}`$ , which is consistent with the OPE because in a dilute instanton liquid we have $`g^2G^2=32\pi ^2(N/V)`$. This term can indeed be seen in the instanton calculation and causes the correlator to drop below 1 at small $`x`$.

It is possible to extract the value of $`\mathrm{\Lambda }_{QCD}`$ together with the power corrections from the data. Because the dimension 6 operator is relatively small we have fixed it from a joint fit with the $`VA`$ correlator. We find

$`{\displaystyle \frac{\mathrm{\Pi }_V(x)+\mathrm{\Pi }_A(x)}{2\mathrm{\Pi }_0(x)}}`$ $`=`$ $`1+{\displaystyle \frac{\alpha _s(x)}{\pi }}\left({\displaystyle \frac{x}{1.52\mathrm{fm}}}\right)^4`$ (8)
$`\left({\displaystyle \frac{x}{0.85\mathrm{fm}}}\right)^6+\mathrm{}.`$
The value of $`\alpha _s(m_\tau )0.35`$  is consistent with other determinations , but the value of the gluon condensate term is smaller than the standard SVZ value . In fact, the data do not show any kind of “dip” and as soon as the $`d=4`$ power correction becomes comparable to the perturbative correction it is in dramatic disagreement with the data. Unfortunately, it will be hard to improve on this situation even if high precision data that cover a larger range of invariant masses in the vector channel become available. Within the range of validity of the OPE in the $`V+A`$ channel, $`x<0.3fm`$, the power corrections are always small as compared to perturbative corrections. This makes it doubtful that one will ever be able to extract the value of the gluon condensate.

5. Finally, we address the purely perturbative contribution to the $`V+A`$ correlation function, using the instanton calculation as a representation of the non-perturbative part of the correlation function. This is supported by the fact that instantons provide a very accurate description of the $`VA`$ correlator which is free of perturbative contributions. The difference between the full correlation function and the instanton calculation is shown by the squares in Fig. 4. For comparison, we also show the full correlation function with only the free field behavior subtracted. At short distance, there is no difference between the two curves, and both follow the first order perturbative result $`\alpha _s(x)/\pi `$. At larger distances $`\mathrm{\Pi }_V(x)+\mathrm{\Pi }_A(x)2\mathrm{\Pi }^0(x)`$ starts to drop, but the non-perturbatively subtracted correlator continues to grow. This behavior nicely shows the running of $`\alpha _s`$ even at moderately large $`x`$ . At $`x<0.3`$ fm the agreement becomes even better if the two-loop contribution is added, but in this case the Landau pole is reached earlier. For this reason, the good agreement of the data with the one-loop result even for large $`x>0.3`$ fm may be somewhat coincidental. The reason one is able to follow the pQCD behavior well outside the usual perturbative domain is the remarkable degree of cancellation among non-perturbative effects. Further high statistics studies of this issue using lattice simulations would be very interesting.

6. In summary, we have used the high statistics ALEPH data on hadronic $`\tau `$ decays to calculate Euclidean space correlation functions in the vector and axial-vector channel. We focussed our discussion on the linear combinations $`\mathrm{\Pi }_V(x)\pm \mathrm{\Pi }_A(x)`$. The combination $`VA`$ receives no contribution from perturbation theory and provides a clean probe for chiral symmetry breaking and the quark condensate. $`V+A`$, on the other hand, allows for a study of perturbative QCD and gluonic operators.

We have compared the two correlation functions with the predictions of the random instanton liquid and the OPE. The instanton model provides a very accurate description of the $`VA`$ correlation function for all distances. In the $`V+A`$ channel the instanton model, supplemented by pQCD corrections with a running coupling constant, also provides excellent description of the data for $`x<1`$ fm.

The remarkable degree of cancellation of non-perturbative effects in the $`V+A`$ channel provides a unique opportunity to access perturbative corrections well beyond the usual pQCD domain. In the $`VA`$ channel, on the other hand, there is an opportunity to extract the dimension $`d=6`$ $`\overline{q}q^2`$ operator from the data. The result agrees with the SVZ prediction, but the accuracy is limited by the largest invariant mass accessible in $`\tau `$ decays. In addition to that, the instanton model suggests the presence of a non-singular $`d=6`$ contribution of the same magnitude. Attempts to extract the $`d=4`$ gluon condensate operator from the $`V+A`$ channel fail because in the range of validity of the OPE the $`d=4`$ power correction remains a small correction to the pQCD contribution.

We conclude that the range of validity of the OPE in the vector channels is quite small, $`x<0.3`$ fm. This means that there is essentially no “window” in which both the OPE is accurate and the correlation function is dominated by the ground state. Instantons, on the other hand, provide a quantitative tool at all distances. This is true even though the vector channels, because of the smallness of direct instanton effects, are generally considered to be the best system to study the OPE.

Acknowledgments: We thank M. Shifman for valuable comments. The work is partly supported by US DOE grant No. DE-FG02-88ER40388.