text:
# Emergence of Hierarchy on a Network of Complementary Agents

## 1 Introduction

In the last years great attention has been paid to the study of different complex systems which share the common feature of organizing themselves in networks. Examples include social networks , the Internet , food webs , metabolic networks , immune networks , economic networks  and ecological networks . Despite their differences, these systems may be cast in a common framework regarding the structure of their connectivity, being mostly classified as either regular, scale-free , small-world  or random networks. Scale-free networks are characterized by power law distributions, indicating there is no typical number of connections per site, whereas the small-world phenomenon is associated with high clustering coefficient and distances between sites that increase only logarithmically with system size , often presenting exponential connectivity distributions .

So far, the models proposed to study the formation and dynamics of networks are mainly based on two approaches: dynamical processes taking place in static network structures and the study of the dynamics of the network structures themselves. Most of the characterization of the network structures is focused on the dynamics of aggregation of new nodes and links to a growing network. A key point in the study of such networks is that their dynamics are based on nodes which lack identity. Or, to put it correctly, their identity and therefore their mutual affinity is uniquely determined by their connectivity. For instance, scale-free networks are obtained by linking newly introduced nodes to old ones with a probability proportional to the number of links of the latter . In other words, it is based on a continuous growth of nodes with preferential attachment. When these complex networks are used as a substrate to study a given dynamical process, the network itself usually remains unchanged while the states associated to the nodes evolve according to pre-determined rules . A mixed case has also been considered in the framework of game theory , where not only the state of the nodes evolve in time, but also their connectivity. In the last case the topological changes are biased by some state-dependent rules.

Although it has been pointed out that many of the network structures mentioned above are scale-free, the assumed rich-gets-richer dynamics that leads to this kind of structure may be appropriate only for some cases . It is not clear, for instance, whether biological systems typically belong to this class of networks, since the addition of new nodes does not necessarily occur continuously, and certainly not without a strong dependence on the nodes identity. Many biological systems organize themselves through pattern recognition based on molecular interactions. Molecules, natural or synthetic, are able to replicate or interact when their shapes and chemical features are complementary. Depending on the way atoms or group of atoms are distributed spatially in a molecule, it may fit into chemical “nooks and crannies” of another . Thus complementarity is essential for pattern recognition interactions underlying biological and biochemical processes. It is the basis of the DNA replication mechanisms and the generation of the immune responses, among many other important biological processes. In the case of the immune system, the control of foreign and harmful elements (antigens) to the organism is based on the recognition between molecular receptors and antigenic determinants through lock-and-key interactions. Another example comes up when we think of symbiotic species in ecological networks as complementary nodes, which generate a mutual reinforcement when they interact. In this sense the role of complementarity can be extended to social and economical systems, where agents need to cooperate with other agents in order to prevail in a competing environment.

Inspired on the complementary interactions observed in the immune system and other biological and biochemical systems , the model we present here may be considered as a third approach to describe the formation of networks. The effective connectivity of the system in this case emerges from the dynamical interaction among the nodes, which is in turn inherently tied to their identity. Using the parlance of socio-economical models, each node is hereafter referred to as an agent. In order for each agent to guarantee its participation in the network, it needs not only to make use of its own capabilities, but it also depends on the support of its complementary peers. The maximum number of peers each agent may have is given by a neighborhood which includes all possible complementary agents within some window of tolerance. Thus, the maximum connectivity of each agent is fixed, although large in most cases. In principle all agents could be active but in effect only a fraction of them remains activated due to intrinsic regulation mechanisms. The intensity of the interactions depends on the current state (hereafter referred to as “activity”) of the elements of the interacting pair, which will turn out to dynamically mold their connectivity distribution. When a pair matches due to the complementarity of their attributes, their activity will increase, thus increasing the possibility of new interactions. Note that this amounts to a mutually reinforcing mechanism which operates on the activities of the nodes, not on their connectivity. The appearance of new agents is introduced by diffusion of activity in agent space, strengthening agents belonging to the similarity neighborhood of the interacting agents. These new agents are chosen by generating random mutations in one of the two interacting agents. The newly added agents, however, will not survive unless they manage to maintain themselves in the network, so that chance in this case is subjected to regulation. As we shall see, the rate at which the activity diffuses in the agent space is essential for the maintenance and stability of the network. Self-regulating mechanisms also take place, preventing agents from increasing their activity too much.

Different from the approach used recently by Castellano et al. , where the similarities between agents play a key role in a social influence model, here we assume that complementarity is the main mechanism in establishing the interactions among agents. The question we address here is which kind of organization is generated by individuals interacting locally by complementarity, having a limited (although large) number of interacting neighbors and whose connections are dynamically established depending on their activity?

In the studies of social behavior it is very difficult to test hypotheses concerning the relationship of individual behaviors and macroscopic regularities observed. So far the main concern in social sciences, especially in game theory and general equilibrium theory , has been focused on the stationary states rather than the dynamical evolution of the system. Our aim in this work is to study the essential micro or local mechanisms that lead to the dynamical emergence of global regularities and properties in a social community where individuals interact due to the complementarity of their attributes.

Surprisingly, the fact that a given environment can support only finite population, referred on the literature as the ecological principle of carrying capacity, emerges naturally in our model. We also obtain highly skewed distributions of activity. If we interpret this activity as power, influence or wealth, the inequalities observed in these distributions are characteristic of the actual human society . Starting with agents having the same kind of preferences and possibilities, but a unique neighborhood, the system will end up with a multi-modal distribution of activities, since the differences produced in their states after each interaction may be amplified with time, generating the so called horizontal inequality phenomenon .

## 2 The model

The identity of each agent is defined by a bit string with length $`b`$, each bit representing for instance the presence or absence of relevant attributes. To each site $`i`$ ($`i=1,\mathrm{},2^b`$) of this $`b`$-dimensional hyper-cubic unit cell we associate a discrete variable $`N_i\{0,1,\mathrm{},N_{max}\}`$ which represents the activity of agent $`i`$. If $`N_i=0`$, we say that the $`i`$-th agent does not take part in the network, or is inactive.

At each time step a parallel update of the individual activities of the entire system is performed. The individual update of a given activity $`N_i`$ follows two steps. In the first one, the activity of agent $`i`$ is decreased with probability $`p_i\text{min}(1,N_i/N_{max})`$. This factor accounts for suppression mechanisms generated by e.g. competition for leadership or resources, restricting the growth of an agent’s activity, which is limited to a maximal value of $`N_{max}`$:

$$P\left(N_i^{}|N_i\right)=\left(\genfrac{}{}{0pt}{}{N_i}{N_i^{}}\right)p_i^{N_iN_i^{}}(1p_i)^{N_i^{}},N_i^{}N_i$$
(1)
where $`N_i`$ and $`N_i^{}`$ are the current and updated states of agent $`i`$, respectively, and $`P\left(N_i^{}|N_i\right)`$ is the probability with which the activity of a site $`i`$ changes from $`N_i`$ to $`N_i^{}`$ in a time step. If the agents were allowed to be indefinitely activated, they would simply accumulate activity and the system might never attain a stationary state. In this sense the limitation on each agents activity effectively amounts to the possibility of new agents being incorporated into the network, according to the rules below. In the absence of interactions the agents activity will thus monotonically decrease, on average, like $`N_i^{}=N_i(1N_i/N_{max})`$. If $`N_i`$ eventually surpasses $`N_{max}`$, the activity of the $`i`$-th agent is reset to zero.

The second step on the individual activity update accounts for the mutually reinforcing interaction among agents, through which activity increases . We say that $`j`$ and $`i`$ are perfect complementary matches if the Hamming distance $`d(i,j)`$ between them (the number of different bits) is $`b`$ (in this case we write $`j=\overline{i}`$). However, an agent $`i`$ can interact not only with its perfect complementary match $`\overline{i}`$, but also with agents that are partially complementary. We define a parameter which measures this mismatch window, so that an agent may interact with agents that have at most a number Maxfit of equal attributes, or similarities.

If agents $`i`$ and $`j`$ are coupled, they will increase the total activity by, at the most, $`n_{ij}`$ units, with a probability $`p_{ij}`$ per unit , where

$`n_{ij}`$ $`=`$ $`N_iN_j`$ (2)
$`p_{ij}`$ $`=`$ $`{\displaystyle \frac{N_iN_j}{\left(\genfrac{}{}{0pt}{}{N_{total}}{2}\right)}}`$ (3)
and $`N_{total}(t)=_{i=0}^{2^b}N_i(t)`$. The agent to which a unit of activity is added will be similar to either $`i`$ or $`j`$: after having chosen one of them randomly, each of its bits may be flipped, with probability $`M`$ per bit (eq. 4). This provision allows for a diffusion of the activity generated in agent space, which will turn out to be essential for an agent’s maintenance in the network for two main reasons: 1) by diffusing their activity to a similarity neighborhood, a given pair of agents strengthen other (similar) agents with which they can further interact, thus creating a supportive environment; 2) diffusion also helps protect an agent from the suppression factor, preventing its activity from rising too fast and reaching $`N_{max}`$. $`M`$ will thus be hereafter referred to as the diffusion rate.

Denoting by $`l_k`$ the agent to which the $`k`$-th unit in the $`(i,j)`$ interaction will be potentially added, the dynamics can be described by the iteration of equations 4-5 below, for $`k=1,\mathrm{},n_{ij}`$ and $`i<j`$ such that $`d(i,j)\text{Maxfit}`$:

$`P(l_k)`$ $`=`$ $`{\displaystyle \frac{1}{2}}M^{d(l_k,i)}(1M)^{bd(l_k,i)}+{\displaystyle \frac{1}{2}}M^{d(l_k,j)}(1M)^{bd(l_k,j)}`$ (4)
$`P\left(N_{l_k}^{\prime \prime }|N_{l_k}^{}\right)`$ $`=`$ $`\delta (N_{l_k}^{\prime \prime },N_{l_k}^{}+1)p_{ij}+\delta (N_{l_k}^{\prime \prime },N_{l_k}^{})\left[1p_{ij}\right],`$ (5)
where $`\delta `$ is the Kronecker delta. Note that the activity $`i`$ and $`j`$ created remains within the pair with probability $`(1M)^b+M^b`$.

Several quantities of interest are kept track of. As one would expect, the stability of the activity of a given agent ultimately depends on the activity of its complementary agents, including perfect and slightly defective mismatches. To capture the instantaneous average correlation between agents and their complements, we therefore introduce a variable which measures the symmetry of the activation in agent space:

$$s(t)=\frac{2^b_iN_i(t)N_{\overline{i}}(t)\left(_iN_i(t)\right)^2}{2^b_iN_i^2(t)\left(_iN_i(t)\right)^2}.$$
(6)
Note that $`s=1`$ when all $`N_i=N_{\overline{i}}`$. We also compute the fraction of agents with nonzero activity $`f(t)=2^b_{\{i|N_i0\}}^{2^b}1`$ and the total activity $`N_{total}(t)`$ (see above) of the network. The latter is very useful from the dynamical point of view to help detect when the system reaches the stationary state.

It is instructive to observe how the populations of agents with low ($`0<N_i/N_{max}<10^2`$), intermediate ($`10^2N_i/N_{max}<10^1`$) and large activity ($`10^1N_i/N_{max}1`$) change as the dynamics evolve, which gives a first coarse-grained description of how activity is allocated in the network. The fractions of active agents ($`N_i0`$) falling into the three logarithmic bins above will be respectively called $`C_0(t)`$, $`C_1(t)`$ and $`C_2(t)`$. Once the stationary state is reached, we study in more detail the activity and connectivity distributions of the system. We define an agent’s local connectivity $`k_i`$ as the number of its interacting peers which are active (i.e. we count the number of active agents within the Maxfit-delimited region). Note that the maximal connectivity of an agent

$$k_{max}=\underset{n=0}{\overset{\text{Maxfit}}{}}\left(\genfrac{}{}{0pt}{}{b}{n}\right)$$
(7)
is finite, but increases rapidly as the system size grows (e.g. $`b=16`$ yields $`k_{max}=697`$ for Maxfit=3 and $`k_{max}=2517`$ for Maxfit=4). Therefore we have an interesting situation where the number of potential connections is very large, yet much smaller than the total number of agents (thus we are far away from mean field interactions). Out of these $`k_{max}`$ potential individuals, an agent will dynamically choose some, depending on how much activity the group manages to gather.

In order to further characterize our analysis on the inequality of the activity distribution in a given time step, we also calculate the Gini coefficient, which is a well known index used in socio-economical studies . The Gini coefficient measures the inequality of the distributions by weighing the total activity (or income, or power, etc) of the activated population by the activity of the poorest elements. Compared to other measures of inequality, it has the advantage of being applicable independently of the functional form of the activity distribution (as opposed to the Pareto index, for instance ). Agents are sorted by their activity and the fraction $`x`$ of the least excited agents (i.e. the cumulative population) is computed. Denoting by $`y`$ the cumulative activity associated to those agents, the function $`y(x)`$ is called the Lorenz curve. For perfectly “egalitarian” distributions ($`N_i=N`$, $`i`$), $`y_{eg}(x)=x`$. The Gini index is twice the area between $`y(x)`$ and $`y_{eq}(x)`$, attaining zero for egalitarian distributions and one if a single agent has all the activity in the system.

## 3 Results

For the sake of clarity in all the cases discussed below we have kept Maxfit = 2 and $`N_{max}=10^3`$. In figures 1 and 2 we show the typical behavior observed for single runs for $`b=8`$ and $`b=12`$, respectively. We randomly select an initial fraction $`f(0)`$ of the total number of possible agents and set their activities to $`1\%`$ of the maximal activity.

The dynamics exhibited for $`b=8`$, although interesting in its own and revealing the appearance of emergent properties, is slightly different from the dynamics observed for larger system sizes (see the results for $`b=12`$ below). Figure 1(a) shows the behavior of the correlation $`s(t)`$ among complementary agents for different values of the diffusion rate $`M`$. For M = 0 there is a very low correlation among complementary agents, since $`s`$ is essentially established by the initial configuration. The activation will basically last among these initially activated agents, due to the lack of diffusion of activity in agent space. Increasing the diffusion rate $`M`$ from $`10^4`$ to $`10^1`$, we observe that the correlation initially increases, reaching a maximum value around one for $`M=10^2`$, in which case we interpret that a network is mounted. For $`M=10^1`$ the correlation oscillates non-periodically, indicating that a large rate of diffusion does not guarantee the stability of the network, which mounts for a while but does not last long. This last kind of behavior is reminiscent of the “global dynamical cascades” observed in a model of heterogeneous agents for economics .

Figures 1(b)-(d) show the result for the fraction of active agents with low, intermediate and large activity for $`M=0,10^4`$ and $`10^1`$ respectively. In the first case \[Fig. 1(b)\], without diffusion of activity in agent space, the system does indeed reach a stationary state, with the overwhelming majority of the active agents with large activity. The fact that no intermediate- or low-activity agents are present turns out to be an indication of failure ($`s0`$) of the network to self-organize hierarchically (note that agents with $`N_i=0`$ do not enter the statistics of the $`\{C_i\}`$). A mild amount of diffusion ($`M=10^4`$) is enough to help mounting a hierarchical structure, as shown in Fig.1(c). The behavior is now completely different from the previous case and we do observe the formation of a network. This is not only characterized by the large value of $`s`$, but also by the stratified concentrations of active agents, which stabilize in the long run. The stationary values suggest a rough balance between low, intermediate and large activity in order to sustain the structure of the network. This pattern persists up to values $`M=𝒪(10^2)`$. The results for $`M=10^2`$ (not shown) indicate the formation of a stronger and highly activated network, with the majority ($`80\%`$) of the agents with intermediate activity, and $`20\%`$ with large activity, showing some hierarchy, where stronger agents are kept by a large number of intermediate-activity agents and essentially none with low activity.

For larger values of the diffusion rate, the formation of the network is destabilized, as shown for $`M=10^1`$ in Figure 1(d). We observe oscillations in the values of $`C_0(t)`$ and $`C_1(t)`$, suggesting that the network sustains itself for a brief period. However, this balance is quickly destroyed by the rate at which activity disperses in the network, which prevents the stronger agents and their complementary partners to satisfactorily keep their dominance. The peak in $`C_1`$ \[$`t1100`$ in Fig. 1(d)\] indicates that the activity which the interacting agents (and their neighbors) create spreads too much and too fast in agent space, eventually destroying their prevalence and returning the network to a non-hierarchical situation, where no particular pair or group of agents dominates the game. In the following it may happen that a new dominant group emerges out of the inherent fluctuations of the model, only to be destabilized later, in a quasi-periodic fashion. This behavior is very robust and has been exhaustively checked over many long runs for different initial configurations. For higher diffusion rates ($`M0.4`$), the peaks described in the previous case disappear and $`s`$ fluctuates around zero, indicating the complete absence of correlation among the agents. In this case $`C_0`$ is essentially one, showing that the system consists of a very large number of weakly activated (and disconnected) agents. It is interesting to point out that in the whole range $`10^4<M<10^1`$ (where the network is mounted), the fraction of active agents $`f`$ remains above $`80\%`$ in the stationary state, indicating that the majority of agents is active.

Results for $`b=12`$ are shown in Figure 2 for two representative values of the diffusion rate. The first interesting feature we observe when $`M=10^3`$ \[Fig. 2(a)\] is that, unlike the 8-bit case, there is a long transient period before the fraction of active agents $`f`$ and the correlation $`s`$ stabilize. After that, the emergence and establishment of the network then occur in a relatively short time. Moreover, $`f`$ remains around $`0.2`$ in the plateau, signaling that the majority of the agents have zero activity and do not take part in the network. In the evolution of concentrations $`\{C_i\}`$ shown in Fig 2(c), most of the active agents have low activity, followed by agents with intermediate and large activity. The same behavior is observed for different samples, although the stationary concentrations may fluctuate around slightly different mean values. The larger 12-bit system is therefore much less excited and sparser than the smaller 8-bit one: not only is the fraction of active agents smaller, but also the activity of those agents is typically less than for the smaller system. Fig. 2(b) shows the results for a higher diffusion rate ($`M=10^1`$) where again we observe a long transient, after which a peak in $`s`$ and in $`f`$ occurs, corresponding to an attempt to mount the network. The attempt fails due to the instability generated by the high rate at which activity diffuses. The system will keep spiking indefinitely, similarly to what has been observed for 8 bits in the same parameter region. Finally, a complementary view of these phenomena can be composed by the results in Fig. 2(d), which shows the evolution of the total activity for both diffusion rates. For $`M=10^3`$ the stabilization of $`s`$ and $`f`$ is associated with the stabilization of $`N_{total}`$ at a large value, while for $`M=10^1`$, the peaks in $`s`$, $`f`$ and $`N_{total}`$ occur simultaneously, the latter attaining a much lower level (note the different vertical scales) than in the $`M=10^3`$ case.

Finally, it is interesting to observe the behavior of the Gini coefficient for both diffusion rates. Figures. 2(a)-(b) show that it departs from zero (since the initial agents have the same activity) and immediately grows, indicating the immediate creation of inequalities due to the fluctuations. Note however that during this growth $`s`$ remains near zero, showing that the initial fluctuations do not immediately lead to a structuring in agent space. The inequality decreases a little and then increases, now lifting with it the correlation. The Gini coefficient then stabilizes around $`0.80`$, indicating a highly uneven distribution.

The diffusion rate thus plays an extremely important role in connecting the agents, its effect being: too low or too high diffusion rates deteriorate the stability of the network, while intermediate values increase the diversity of active agents in the right amount to allow new pairs to match, thus providing the mechanism for the system to adapt, incorporating new agents. Note that in the current model this is the only source providing new active agents. The parameter Maxfit plays a complementary role to that of the diffusion rate on the regulation of the network. A large diffusion rate may eventually be compensated by a smaller Maxfit, which renders the agents more specific in the interaction with their mirrors (in the sense that it reduces the number of interacting peers). However, Maxfit cannot be too low lest there be not enough agents to create a supportive neighborhood, specially for low diffusion rates.

The dependence of the formation of the network for $`b=12`$ on the diffusion rate is summarized in Fig. 3, which shows the stationary value of $`s`$, averaged over 500 time steps and 50 runs, as a function of $`M`$. In order to obtain this plot we have adopted the following protocol: we consider that the network is mounted if $`s`$ stays larger than a threshold value $`s_{min}=0.8`$ (estimated from many individual runs for the parameters under consideration) during a period of 1000 time steps. This time interval is long enough for an eventual peak not to be erroneously detected as a successful attempt. The “mount probability” is defined as the fraction of times that this criterion is satisfied. We observe then two transitions, one at $`M10^4`$ and another at $`M0.03`$. The large error bars in $`s`$ around the transition points are the signature of the peaks observed in Fig. 2. The first transition is mainly provoked by the suppression factor: since the rate at which activity disperses is very low, the initially correlated agents reinforce themselves, leading to a fast increase in their activity which strongly submits them to the suppression factor. In this case no correlation is maintained long enough to mount the network. The second non-equilibrium phase transition is due to the dynamics of the system as discussed below.

Another interesting feature that emerges from this model is the weak dependence on the initial conditions. For several values of $`M`$ we have tested the evolution of the system for different initial fractions of nonzero agents. What we have observed, for both 8 and 12 bits, is that the results are qualitatively the same for whatever initial conditions we choose, except of course for the limiting case $`f(0)0`$, in which the activities eventually decay to zero due to a lack of matching pairs. This result is clearly observed in the inset of Fig. 3, which shows the stationary value of $`s`$ as a function of $`f(0)`$ for 12 bits. The mechanisms behind this independence can be understood on the basis of the transients of both $`N_{total}`$ and $`s`$ (Fig. 2). In order for the network to stabilize, the unpaired initial agents will first be “cleaned up”, leaving mostly pairs that match within the margins set by Maxfit. The dynamics automatically takes care of this initial dismantling, since all those agents which do not have their matching pairs will gradually decrease their activity \[see Fig. 2(d)\]. This in turn increases the probability that the remaining matching agents interact, and stationarity is then achieved when activity increases and losses are balanced. This explains why the model is not very sensitive to the initial conditions. Larger initial values of $`f`$ only increase the time spent in the transient, therefore we have typically used $`f(0)=0.05`$ for 12 bits. In the 8-bit case, the initial transient is much shorter and of a different nature because the network somehow manages to accommodate a much larger fraction of agents in highly excited states, showing a different dynamics. Due to the small system size all potential pairs will be effectively connected, which is not the case for larger system sizes. This particular behavior might be related to what is observed in population models and commonly referred to as finite population effect .

It is also important to understand how the behavior of the system changes when the size of the space gets larger. Representative simulations for $`b=16`$ (using $`f(0)=0.015`$) show that the fraction of active agents $`f`$ in the stationary state is even smaller than for 12 bits, but the behavior of the concentrations $`\{C_i\}`$ remains qualitatively the same, indicating that those are the relevant quantities for larger system sizes. The dynamical behavior for 16 bits exhibits the same features as for 12 bits, the corresponding plots (not shown) being very similar to Fig. 2. Even though a more detailed analysis of the parameter space is currently under investigation, it is clear that the effects of the diffusion rate $`M`$ follow the same trends both for larger systems and for smaller systems: if it is either too large or too small, no balance can be achieved and the agents will fail to establish a connected network. It is interesting to note that for Maxfit=2 the absolute number of active agents in the stationary state remains $`𝒪(10^3)`$ when one goes from 12 to 16 bits, explaining the decrease in $`f`$. This number grows when Maxfit=3, but still corresponds to a small fraction of the $`2^b`$ potential agents, for large systems, and it is the diffusion rate that guarantees its stability. Without any direct suppression mechanism acting on the number of agents belonging to the network, the local rules of the model therefore lead to an effective “carrying capacity” of the system . The proper scaling of Maxfit with the growth of the system is yet to be fully understood, i.e., how does the stationary value of $`f`$ depend on the relation between Maxfit and b? We also conjecture whether a larger agent space would allow the system to adapt more easily to perturbations, which normally would destabilize the network for small systems. These and other questions concerning the robustness of the network on the parameter space are currently under investigation and will be published elsewhere.

We finally turn to the characterization of the ordered phase, focusing at intermediate values of $`M`$. Our analysis are based on the emergent distribution of activity and connectivity in the network, which turns out to be highly non-trivial. $`P(N)`$ measures how activity is allocated, whereas $`P(k)`$ is the connectivity distribution. These distributions are calculated in the stationary state for $`b=8`$ ($`M=10^2`$) and for $`b=12`$ and $`b=16`$ ($`M=10^3`$) over 30 runs. Results are shown in Fig. 4. For the three cases \[Figs. 4(a)-(c)\], we obtain a multi-modal distribution of activity. Therefore, the system does not have a single typical activity, nor does it possess a scale-free activity distribution: as opposed to that, we observe that there are four well-defined “typical” values of activity (other than zero) for 8 and 12 bits, and three for 16 bits, defining a hierarchy in agent space. We therefore justify in hindsight the arbitrary logarithmic binning as defined by the variables $`C_0`$, $`C_1`$ and $`C_2`$, which are actually reasonable coarse grained descriptions of the peaks of Fig. 4. It is also worth noting that, for 12 and 16 bits, there is a gap in the activity distribution, a feature that has also been observed in larger system sizes.

The distribution presented in Fig. 4(d) exhibits an even richer behavior for 12 bits. Not only do we observe multi-modality again, but also a non-monotonic behavior of the peaks. In other words, the majority of the agents are “connected” to $`40\%`$ or $`50\%`$ of their potential collaborators, while two smaller groups either have $`20\%`$ or $`90\%`$ of their possible connections actually established. The distribution for 16 bits has similar features, while for 8 bits $`P(k)`$ in its most ordered regime is simply given by $`P(k)=\delta (k,k_{max})`$ (where $`k_{max}`$ is defined in eq. 7), reflecting the fact that $`f`$ saturates at one. The increase of the variance of the activity distribution corresponds to an increase of the horizontal inequality. It depends on the rate at which new agents are activated (diffusion of activity in agent space), since new agents take some time until their activity is brought into line with those already prevailing in the system.

## 4 Conclusions

In this paper we proposed a network model of mutually reinforcing agents by using a bit-string approach. The model is inspired in the complementary interactions observed in biological and biochemical systems. The interaction among the agents is based on their identity. All agents have potentially the same connectivity although the actual individual connectivities are established by dynamical rules based on interactions that allow perfect and slightly defective mismatches between agents. A suppression mechanism acting upon agents limits the maximum activity that may be associated to each agent. New agents may be added to the system through the interactions among active agents, which allow the diffusion of activity in agent space. The rich dynamical behavior exhibited by this kind of network depends on the diffusion rate, the size of the mismatch window and the maximum activity allowed to each agent. The agents prospering in activity may activate new agents. For very low or very large diffusion rates there is a lack of correlation between complementary agents and the network is not formed, while for intermediate values of this parameter we observe the formation of a stable network.

For small systems (8 bits) we observe a highly excited network involving most of the available agents, while for larger systems this network is much sparser. However, regardless of the size of the system, we observe a hierarchical organization reflected by the multi-modal distribution of activities and connectivities. When we increase the size of the system from $`12`$ to $`16`$ bits, we observe that the absolute number of agents participating in the network has the same order of magnitude. Therefore the number of competent inactive agents increases as the system size grows, a feature that might favor the adaptation of large networks to perturbations.

The skewed distributions obtained correspond to the emergence of structures or self-organization in the network of activated agents, induced by local interactions between the agents. The system organizes itself in a hierarchical structure where we find grosso modo three classes: highly, intermediate an low activated agents. Since we considered the network of active agents, we do not include among the classes the non-active one. For the parameters considered in this paper the Gini coefficient stays always between $`0.700.80`$, indicating a rather unequal distribution of activities.

In the stationary state we do observe the organization of the system with a finite number of active agents, without any need of ab initio considerations about finite population. This means that at this point there is an equilibrium between agents being deactivated and those being activated. The model may explain the hierarchy observed in business networks and social networks. It may be also a good model to describe the formation of corruption networks. For corruption to take place, the interaction among the agents in this case cannot be based on similarities or acquaintances, but rather on complementarity that generates a dependence among the participants. Diffusion can then be understood as a way of sharing power in order to protect leading agents from betrayal. As another application we may consider the introduction of different types of interactions (other than mutually reinforcing), which might give rise to multi-modal distributions which could then be interpreted as the emergence of trophic levels in ecology.

### Acknowledgments

This work was partially supported by CNPq, CAPES and FAPERJ. We would like to thank Nestor Caticha and Américo T. Bernardes for their invaluable suggestions and criticisms.


format: text

text:
# Generalized ensemble algorithm for U(1) gauge theory

## 1 INTRODUCTION

The four dimensional pure compact U(1) gauge theory is known to posses a phase transition (PT) separating a confined phase and a Coulomb one and the determination of the order of PT is of great importance. In the pure compact U(1) gauge theory, however, the determination of the order of PT is very difficult since on large lattices the standard updating algorithms like Metropolis, heat-bath, hybrid Monte Carlo (HMC), fail to generate enough tunneling between metastable states. Many efforts have been done to clarify the order of the PT. In spite of such efforts the order of the PT is still controversial.

A promising algorithm to overcome this difficulty is the multicanonical algorithm which uses a multicanonical weight in stead of the Boltzmann one and can enhance the tunneling rate between metastable states. Although the multicanonical algorithm works effectively, one disadvantage of the multicanonical algorithm might be that a multicanonical weight used in the simulations is not known a priori and it must be determined before the simulations. Usually the weight is estimated from a short run of the standard updating algorithm. However this estimation is more difficult as a lattice size becomes larger.

Recently in several fields where the usual Monte Carlo technique does not work effectively, simulations were performed with the Tsallis weight, first introduced by Tsallis. The Tsallis weight controlled by a parameter $`q`$ can be easily defined, and the usual Boltzmann weight is given by taking the limit $`q1`$. The results have showed that the use of the Tsallis weight might be more advantageous than that of the Boltzmann one.

Inspired by such studies, here we apply the Tsallis weight for the hybrid Monte Carlo simulation of the pure compact U(1) lattice gauge theory and investigate whether the Tsallis weight solves the difficulty in determining the order of the PT.

## 2 GENERALIZED ENSEMBLE AND TSALLIS WEIGHT

The partition function of a system with the action $`S`$ is given by

$$Z=\mathrm{exp}(\beta S[U])𝑑U.$$
(1)
Introducing a function $`G`$, the partition function can be rewritten as

$`Z`$ $`=`$ $`{\displaystyle \mathrm{exp}(\beta S[U])/GG𝑑U}`$ (2)
$`=`$ $`<{\displaystyle \frac{\mathrm{exp}(\beta S)}{G}}>_GZ_G,`$
where $`Z_G`$ is a generalized partition function:

$$Z_G=G𝑑U,$$
(3)
and $`<O>_G`$ stands for the expectation value of an observable $`O`$ in the ensemble generated with $`Z_G`$. The expectation value of $`O`$ in the canonical ensemble generated with eq.(1) can be obtained in terms of the generalized ensemble as

$`<O>`$ $`=`$ $`{\displaystyle O\mathrm{exp}(\beta S[U])𝑑U/Z}`$ (4)
$`=`$ $`<O{\displaystyle \frac{e^{\beta S}}{G}}>_G/<{\displaystyle \frac{e^{\beta S}}{G}}>_G.`$
Previously it was suggested that in some cases Tsallis weight works better than the Boltzmann one. In the Tsallis formulation,

$$G=[1(1q)\beta \overline{S}]^{\frac{1}{(1q)}},$$
(5)
or in the familiar exponential form,

$$G=\mathrm{exp}(\beta S_T)W_T,$$
(6)
with

$$S_T=\frac{1}{\beta (q1)}\mathrm{ln}[1(1q)\beta \overline{S}],$$
(7)
where $`\overline{S}=SS_0`$. $`S_0`$ is a certain constant shifting the origin of the action. Here we call $`W_T`$ ”Tsallis weight”. The Boltzmann weight can be obtained by taking the limit $`q1`$ in the Tsallis weight $`W_T`$.

## 3 Hybrid Monte Carlo with Tsallis weight

The action of pure compact U(1) lattice gauge theory is given by

$$S=\underset{x,\mu >\nu }{}[1\mathrm{cos}(\theta _{\mu \nu }(x))],$$
(8)
where $`\theta _{\mu \nu }(x)`$ is the sum of link angles contributing a plaquette at a site $`x`$ on the $`\mu \nu `$ plane. The Hamiltonian with the Tsallis weight, used for HMC simulations, is defined by

$$H_T=\frac{1}{2}p_i^2+\beta S_T.$$
(9)
Here note that $`S_T`$ is not always well-defined. Namely $`1(1q)\beta \overline{S}`$ can be negative. However, provided that $`q`$ is close to one, $`1(1q)\beta \overline{S}`$ is made positive.

## 4 RESULTS

As an exploratory study we tested the above algorithm varying $`q`$ on an $`8^4`$ lattice in the vicinity of the PT ( at $`\beta =1.00737`$ ). We performed HMC simulations with eq.(9). $`S_0`$ ( = 0.37 ) was taken to be the average value of the action $`S`$. In order to make comparisons, we also performed HMC simulations with the Boltzmann weight. Typically we simulated $`10^6`$ trajectories for each $`q`$.

We observed more tunneling for $`q<1`$ than for $`q=1`$ (the Boltzmann weight). Fig.1 shows representative history of actions of q=1 and 0.999997 and we see more tunneling for the case of q=0.999997.

Fig.2 shows action densities of q=1, 0.999997, 0.999995 and 0.99999. The height of the valley between two peaks increases as $`q`$ decreases which indicates that the tunneling rate increases as $`q`$ decreases. For too small $`q`$, however, values of the action will locate near $`S_0`$ and in this case there may be a difficulty in reweighting. Fig.3 and 4 show that the action densities reweighted to q=1 and compare those with the Boltzmann action density ( $`q=1`$ ). The reweight densities reproduce well the Boltzmann action density. However, for $`q=0.99999`$ ( the smallest $`q`$ here ) the reweighted action density has large fluctuations at $`S`$ far from $`S_0`$, which may indicate that one needs caution when using the Tsallis weight.

It is also possible to take $`q>1`$ in the simulations. For $`q>1`$, however, the desired behavior, i.e. more tunneling, was never observed for our case.

## 5 DISCUSSION

The exploratory study on an $`8^4`$ lattice showed that the Tsallis weight enhances the tunneling rate. The lattice size used here is small. For precise determination of the order of the PT, larger lattices are needed. Therefore it is important to confirm that this preferred feature of the Tsallis weight is preserved for larger lattices, possibly up to $`18^4`$ lattice size or so. Simulations on larger lattices are in progress.

## ACKNOWLEDGMENTS

This work was partially supported by the Ministry of Education, Science, Sports and Culture, Grant-in-Aid, No.13740164.


text:
# MPI-PhT-2001-20 PM/01-34 Hints on the power corrections from current correlators in 𝑥-space

## 1 Introduction

We shall be concerned here with the current-current correlators in the coordinate space:

$$\mathrm{\Pi }(x)=0|J(x)J^{}(0)|0,$$
(1)
in case of the $`(V\pm A)`$ and (pseudo)scalar currents:

$`J_\mu ^{(V\pm A)}`$ $`=`$ $`\overline{q}_i\gamma _\mu (1\pm \gamma _5)q_j,`$
$`J^{S\pm P}`$ $`=`$ $`[(m_im_j)\pm (m_i+m_j)]\overline{q}_i(1\pm \gamma _5)q_j,`$ (2)
where $`q_{i,j}`$ and $`m_{i,j}`$ are the quark fields and masses. The two-point functions (1) obey a dispersion representation:

$$\mathrm{\Pi }(x)=\frac{1}{4\pi ^2}_0^{\mathrm{}}𝑑t\frac{\sqrt{t}}{x}K_1(x\sqrt{t})\mathrm{Im}\mathrm{\Pi }(t),$$
(3)
where Im$`\mathrm{\Pi }(t)`$ is related to the current induced cross section and $`K_1(z)`$ is the modified Bessel function, which behaves for small $`z`$ as:

$$K(z0)\frac{1}{z}+\frac{z}{2}\mathrm{ln}z.$$
(4)
In the limit $`x0`$, $`\mathrm{\Pi }(x)`$ coincides with the free-field correlator and the main theoretical issue is how the asymptotic freedom gets violated at intermediate $`x`$.
From pure theoretical point of view, the use of the $`x`$-space is no better than the use of the momentum space, which is the traditional tool of the QCD sum rules . Each representation has its own advantages and inconveniences (for a recent discussion see ). The $`x`$-space approach is motivated and described in detail in Ref. . In particular, the current correlators (1) are measured in the most direct way on the lattice. The importance of the lattice measurements  is that they allow to study the correlation functions for currents with various quantum numbers, while direct experimental information is confined to only vector and axial-vector currents . The well-known $`\tau `$-decay data were widely used for theoretical analyses both in the $`Q`$\- and $`x`$-spaces ( see, e.g.,  ). Most recently, new lattice data on the $`S,P`$ channels were obtained . The most interesting observation is that in the $`S+P`$ channel there are noticeable deviations from the instanton liquid model  while in the $`V\pm A`$ channels the agreement of the existing data with this model is quite good .
Such deviations were in fact predicted in Ref.  where unconventional quadratic corrections, $`1/Q^2`$ were introduced. The primary aim of the present note is to perform a more detailed comparison of the lattice data with the model of Ref. . We, indeed, find further support for the quadratic corrections. However, the overall picture is far from being complete and we are trying to analyze the data in a more generic way. The central assumption is that the violations of the parton model for the correlators at moderate $`x`$ are due to power-like corrections.

## 2 Current-current correlators

For the sake of completeness, we begin with a summary of theoretical expressions for the current correlators, both in the $`Q`$ and $`x`$spaces. We will focus on the $`(V\pm A)`$ and $`(S\pm P)`$ channels since the recent lattice data  refer to these channels. In case of $`(V\pm A)`$ currents the correlator is defined as:

$$\mathrm{\Pi }_{\mu \nu }(q)=id^4xe^{iqx}TJ_\mu (x)J_\nu (0)^{}=(q_\mu q_\nu g_{\mu \nu }q^2)\mathrm{\Pi }(q^2),$$
(5)
where $`q^2Q^2>0`$ in the Euclidean space-time. For the sake of definiteness we fix the flavor structure of the light-quark current $`J_\mu `$ as:

$$J_\mu ^{V\pm A}=\overline{u}\gamma _\mu (1\pm \gamma _5)d.$$
(6)
In the chiral limit one has in the $`(V+A)`$ case (see, e.g., ):

$$\mathrm{\Pi }^{V+A}(Q^2)=\frac{1}{2\pi ^2}\left\{\left(1+\frac{\alpha _s}{\pi }\right)\mathrm{ln}\frac{Q^2}{\nu ^2}\frac{\alpha _s}{\pi }\frac{\lambda ^2}{Q^2}+\frac{\pi }{3}\frac{\alpha _s(G_{\mu \nu }^a)^2}{Q^4}+\frac{256\pi ^3}{81}\frac{\alpha _s\overline{q}q^2}{Q^6}\right\}.$$
(7)
The corresponding relation for the $`(VA)`$ case reads as:

$$\mathrm{\Pi }^{VA}(Q^2)=\frac{64\pi }{9}\frac{\alpha _s\overline{q}q^2}{Q^6},$$
(8)
In the x-space the same correlators, upon dividing by $`\mathrm{\Pi }_{pert}^{V+A}`$ where $`\mathrm{\Pi }_{pert}^{V+A}`$ stands for the perturbative correlator, are obtained by applying the equations collected for convenience in the Table 1.

In particular,

$$\frac{\mathrm{\Pi }^{V+A}}{\mathrm{\Pi }_{pert}^{V+A}}1\frac{\alpha _s}{4\pi }\lambda ^2x^2\frac{\pi }{48}\alpha _s(G_{\mu \nu }^a)^2x^4\mathrm{ln}x^2+\frac{2\pi ^3}{81}\alpha _s\overline{q}q^2x^6\mathrm{ln}x^2.$$
(9)
Note that $`\mathrm{ln}x^2`$ is negative since we start from small $`x`$. An important technical point is that on the lattice one measures the trace over the Lorentz indices $`\mu ,\nu `$, see Eq (5). In the $`Q`$-space this is equivalent to considering $`Q^2\mathrm{\Pi }(Q^2)`$ instead of $`\mathrm{\Pi }(Q^2)`$. The $`x`$-transform of the $`Q^2\mathrm{\Pi }(Q^2)`$ is given by:

$$\frac{Q^2\mathrm{\Pi }^{V+A}}{Q^2\mathrm{\Pi }_{pert}^{V+A}}1\frac{\pi }{96}\alpha _s(G_{\mu \nu }^a)^2x^4+\frac{2\pi ^3}{81}\alpha _s\overline{q}q^2x^6\mathrm{ln}x^2.$$
(10)
Next, we will concentrate on the currents having the quantum numbers of the pion and of $`a_0(980)`$-meson. The correlator of two pseudoscalar currents is defined as

$$\mathrm{\Pi }^P(Q^2)id^4xe^{iqx}T\{J^\pi (x)J^\pi (0)\},$$
(11)
where

$$J^P=i(m_u+m_d)\overline{u}\gamma _5d,$$
(12)
In the momentum space, to represent the result in terms of the running coupling, and masses it is more convenient to consider the second derivative in $`Q^2`$ of $`\mathrm{\Pi }^P(Q^2)`$ defined in Eq. (11), which obeys an homogeneous RGE:

$$\frac{^2\mathrm{\Pi }^P}{(Q^2)^2}=\frac{3}{8\pi ^2}\frac{(\overline{m}_u+\overline{m}_d)^2}{Q^2}\left\{1+\frac{11}{3}\frac{\overline{\alpha }_s}{\pi }\frac{4\alpha _s}{\pi }\frac{\lambda ^2}{Q^2}+2\frac{\pi }{3}\frac{\alpha _s(G_{\mu \nu }^a)^2}{Q^4}+23\frac{896\pi ^3}{81}\frac{\overline{\alpha }_s\overline{\overline{q}q}^2}{Q^6}\right\}.$$
(13)
Here, the standard OPE terms can be found in  while the gluon-mass correction was introduced first in  111We assume that $`\alpha _s\lambda ^2`$ does not run like $`\alpha _s(G_{\mu \nu }^a)^2`$..

In what follows, we shall work with the appropriate ratio where the pure perturbative corrections are absorbed into the overall normalization and concentrate on the power corrections assuming that these corrections are responsible for the observed rather sharp variations of the correlation functions. Thus, in the x-space we have for the pion channel:

$$\frac{\mathrm{\Pi }^P}{\mathrm{\Pi }_{pert}^P}1\frac{\alpha _s}{2\pi }\lambda ^2x^2+\frac{\pi }{96}\alpha _s(G_{\mu \nu }^a)^2x^4\frac{7\pi ^3}{81}\alpha _s\overline{q}q^2x^6\mathrm{ln}x^2.$$
(14)
Note that the coefficient in front of the last term in Eq. (14) differs both in the absolute value and sign from the corresponding expression in . The channel which is crucial for our analysis is the $`(S+P)`$. In this channel:

$`R_{P+S}{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mathrm{\Pi }^P}{\mathrm{\Pi }_{pert}^P}}+{\displaystyle \frac{\mathrm{\Pi }^S}{\mathrm{\Pi }_{pert}^S}}\right)1{\displaystyle \frac{\alpha _s}{2\pi }}\lambda ^2x^2+{\displaystyle \frac{\pi }{96}}\alpha _s(G_{\mu \nu }^a)^2x^4+{\displaystyle \frac{4\pi ^3}{81}}\alpha _s\overline{q}q^2x^6\mathrm{ln}x^2.`$ (15)
This expression concludes the summary of the power corrections to the current correlators.

## 3 Quadratic power corrections.

One of the central points of the present note is that there are no $`\lambda ^2`$ corrections to the $`V\pm A`$ correlators as can be seen in Eq (10). On the other hand, these terms are present in the case of the $`(S\pm P)`$ channels as can be seen in Eq (15). There is no such asymmetry in the $`Q`$-space, see Eqs (7), (8) and (11). Thus, the $`(V\pm A)`$ correlator, as measured on the lattice, are $`\lambda ^2`$-term blind 222Note that a reversed case is well known. Namely, there are certain terms which are seen in the small-$`x`$ expansion and are not seen in the large-$`Q`$ expansion . ! Thus, we are coming to a kind of a theorem. Namely:

If one assumes:
(1) that the $`V\pm A`$ channels are described by the instanton liquid model while in the $`S+P`$ channel there are considerable deviations from this model (as the lattice data seem to strongly indicate )
and
(2) that this difference is due to some power corrections,
then
the power corrections can be uniquely identified as the gluon-mass corrections (see terms proportional to $`\lambda ^2`$).

Note that the $`\lambda ^2`$ corrections are singled out for two reasons: First, since taking the trace over the Lorentz indices $`\mu ,\nu `$ corresponds to multiplying by $`Q^2`$ as can be seen in the discussion of Eq. (10) and it is only a $`1/Q^2`$ correction which can become a polynomial as a result of multiplying $`\mathrm{\Pi }(Q^2)`$ by $`Q^2`$. Second, in the $`Q`$-space there should be no log factor in front of $`1/Q^2`$, $`1/Q^2\mathrm{ln}Q^2`$. These two conditions are satisfied in case of $`(V\pm A)`$ currents and are not fulfilled in the $`(S\pm P)`$ channels. The difference between the channels is that in the latter case the $`\lambda ^2`$ correction is present in the imaginary part of the $`\mathrm{\Pi }^{S,P}(Q^2)`$ .
Thus, if one retains only the $`\lambda ^2`$ corrections, then, there are no violations of the parton (perturbative) picture in the $`(V\pm A)`$ channels for the correlator measured in  while the violations are present in the $`(S\pm P)`$ channels.
Of course, this limiting case is not necessarily describing the reality and we proceed to quantitative fits to the data .

## 4 Analysis of the data

In Fig. 1 we confront the OPE predictions with the lattice data on the $`(S+P)`$ channel obtained in . The choice of the $`(S+P)`$ channel is motivated by the fact the single instanton contribution cancels from this channel 333The instanton contribution is not large also in the $`(V\pm A)`$ channels. However the $`\lambda ^2`$ terms are canceled from these channels, see the discussion in section 3, and we cannot add anything to the analysis of Ref. . and it was predicted in Ref.  that the $`\lambda ^2`$ correction will be manifested in this channel. The theoretical curves in Fig. 1 correspond to two sets of values of the condensates given in Table 2. The first set (SET 1) corresponds to the standard SVZ values of the gluon and four-quark condensate, the latter being obtained using the vacuum saturation assumption. The second set (SET 2), corresponds to the values of the condensates obtained in , where the value of the gluon condensate is two times the SVZ value and the four-quark condensate exhibits a violation of the vacuum saturation, first obtained from $`e^+e^{}`$ data in . In SET 3, one also accounts the presence of the new 1/$`Q^2`$-term first advocated in  and fitted from $`e^+e^{}`$ data in  444A common difficulty encountered in determining the quadratic corrections is that they usually compete with the standard perturbative radiative corrections. In , a suitable choice of the sum rules (e.g. ratio of moments) has been used such that the perturbative radiative corrections are eliminated to leading order and the contribution of the quadratic term becomes optimal. Moreover, the quadratic corrections corresponding to $`\lambda ^20.5GeV^2`$ do not affect in a significant way the determination of $`\alpha _s`$ from $`\tau `$-decay as has been explicitly shown in . On the contrary, the quadratic term appears to decrease very slightly $`\alpha _s`$ from the $`\tau `$-decay and bring it closer to the world average value at $`M_Z`$. In Ref.  bounds were obtained on the value of $`\lambda ^2`$ from the sum rules which have large perturbative terms. This turned possible due to a particular fixation of the of the perturbative terms in the complex $`q^2`$ plane. In particular, the common use of the running coupling would affect the procedure strongly .. Note also that in numerical fits we put $`\mathrm{ln}x^2=1`$, the same, as, say, in Ref..

The analysis indicates that a much better fit of the lattice data for the $`S+P`$ channel at moderate values of $`x`$ is achieved after the inclusion of the $`1/Q^2`$, or $`x^2`$ quadratic correction. A caveat is that we account only for the power corrections, not pure perturbative contributions. The reason is that the lattice data, in their present status, do not give any clear indication of the perturbative contributions. Note also that the data cannot discriminate between the values of the dimension four and six condensates entering in SET 1 and SET 2 as the effects of these two condensates tend to compensate each other for the choice $`\mathrm{ln}x^21`$. The agreement of the OPE with the lattice data at larger values of $`x`$ can be obtained by the inclusion of the $`D=8`$ condensate with a size $`+(x/0.58)^8`$ where we have used $`\mathrm{ln}^2x^21\mathrm{ln}x^2`$. This value can be compared with the one $`+(3395/30855168)\alpha _sG^2^2x^8(x/1.2)^8`$, which one would obtain from the evaluation of these contributions in  and where a modified factorization of the gluon condensates proposed in  has been used.

For completion, we show in Fig. 2, a fit of the lattice data in the V+A channel using SET 3 values of the gluon and quark condensates and quadratic term plus a $`D=8`$ contribution with the strength $`(x/0.7)^8`$ to be compared with the one $`\alpha _sG^2^2x^8/3428352(x/2.5)^8`$ which one would obtain using the results in . Both fits in Figs 1 and 2 might indicate that the vacuum saturation can be strongly violated for higher dimension condensates, a feature already encountered from different analysis of the $`\tau `$ and $`e^+e^{}`$ data . Therefore, we would also expect analogous large deviations in the V-A channel.

## 5 Discussions. Two-step QCD

While evaluating the emerging picture of the power corrections, one should face the possibility that the standard OPE (see, e.g., ) is valid only at very short distances. What is even more important, the mass scale where higher terms in the OPE become numerically comparable to the lowest ones is not necessarily the scale associated with the resonances but could be considerably higher. There is accumulating evidence to support such a view:
(1) A direct comparison of the OPE with the lattice data in the $`(VA)`$ channel demonstrates that the convergence radius of the OPE is no larger than 0.3 fm .
(2) Within the instanton liquid model , the distance between instantons is a few times larger than the size of the instantons. On the physical grounds, the OPE applies at distances smaller than the instanton size while the resonance properties are rather related to distance between the instantons (if encoded in the model at all). Respectively, neither the lattice data nor the predictions of the instanton liquid model exhibit any irregularity at the convergence radius of the OPE .
(3) Within the monopole-dominated-vacuum model the two scales are even more pronounced numerically: the monopole radius is about 0.06 fm while the distance between the monopoles is about 0.5 fm .
(4) If one replaces the local condensates of the standard OPE by their non-local counterparts (for a review see ), then the effect of the non-locality is strong already at $`(0.1÷0.3)`$ fm .

If, indeed, the validity of the standard OPE derived within the fundamental QCD is shrunk to very short distances, then improving within the standard OPE fits to the data obtained at presently available lattices might not be a proper criterion for selecting the right model.

Instead, there emerges a picture according to which the non-standard quadratic power corrections dominate the presently available “intermediate distances ” of about $`(0.1÷0.5)`$ fm. There are two pieces of evidence of the quadratic corrections dominating over the whole range of intermediate distances indicated above:
(1) non-perturbative contribution to the heavy quark potential is linear at all the distances $`r>0.1`$ fm (for discussion see ).
(2) Instanton density, as function of the instanton size $`\rho `$ is reproduced at all the distances $`\rho >0.1`$ fm .
Moreover,
(3) Introduction of the tachyonic gluon mass explains in a simple and unified way existence of the strong channel dependence of the scales of the violation of the asymptotic freedom .
To this list we add now a new observation:
(4) Instanton-liquid model plus the $`\lambda ^2`$ correction gives a reasonable fit in the $`(S+P)`$ channel at distances $`(0.1÷0.5)`$ fm.

As for the status of the $`\lambda ^2`$ correction it is most solid within the effective Higgs-like theories which are common within the monopole mechanism of confinement. Indeed, in the presence of a magnetically charged (effective) scalar field the symmetry of the theory is $`SU(3)_{color}\times U(1)_{magnetic}`$ . Upon the spontaneous breaking of the magnetic $`U(1)`$ the gauge boson acquires a non-vanishing mass and its mass squared is the only parameter of dimension $`d=2`$ consistent with the symmetry. Moreover, in exchanges between (color) charged particles the gauge-boson mass appears to be tachyonic mass as was demonstrated on the $`U(1)`$ example in Refs. . Detailed analysis of various power corrections within the Higgs-like models can be found in Refs. . Moreover, if the monopole size is indeed as small as indicated above, then the effective Higgs-like theories can apply at all the distances $`(0.1÷0.5)`$ fm.

To summarize, the $`\lambda ^2`$ corrections introduced in  rather drastically improve agreement of the theoretical predictions for the current correlator in the $`(S+P)`$ channel with the lattice data (without affecting the other channels measured), in the moderate $`x`$-region, and which can be extended to higher $`x`$-values by including higher dimension condensates. The main uncertainty of the analysis is due to neglect of the pure perturbative corrections, not detected so far on the lattice. Further checks could be provided by measuring the current correlators in the other channels discussed in . The success of the fits with $`\lambda ^2`$ corrections included can be understood within the effective theories.

## Acknowledgements

The authors are grateful to Th. De Grand, F.V. Gubarev, E.V. Shuryak for useful discussions and communications. S. Narison wishes to thank W-Y. Pauchy Hwang for the hospitality at CosPA-NTU (Taipei), where some part of this work has been done. V. Zakharov would like to acknowledge gratefully the hospitality of the Centre National de la Recherche Scientifique (CNRS) during his stay at the University of Montpellier when this was begun.