text:
# Time evolution of damage under variable ranges of load transfer

## I Introduction

The phenomenon of fracture in heterogeneous materials is a complex physical problem which has been of great interest to scientists for quite a long time r1 ; r2 ; r3 . A disordered system is understood to be one with random time or/and space dependent breaking properties r1 . The random nature of this dependence arises because modeling the fracture of heterogeneous materials entails dealing with systems made up of many constituents, each one having mechanical properties that can be considered as being independent of the properties of the other constituents, but with many body interactions among the different parts of the system r1 ; r2 ; r3 . In heterogeneous materials the evolution of the rupture process is radically different from the single crack growth mechanism that occurs in homogeneous materials. Though there is neither a complete numerical solution nor analytic solution to the fracture problem, we have a better understanding of it due to some recent algorithms that have been developed to simulate the fracture process r4 ; r5 ; r6 ; r7 ; r8 .

Fiber Bundle Models (FBM) form a fundamental class of approaches to the fracture problem. They arose in close connection with Daniels’ and Coleman’s seminal works on the strength of bundles of textile fibers r9 ; r10 , and have harbored an intense research activity in recent years r1 ; r2 ; r3 ; r4 ; r11 ; r12 ; r13 ; r14 ; r15 ; r16 ; r17 ; r18 ; r19 ; r19b ; r19c . Fiber Bundle Models are important, despite their very simple nature, because they exhibit most of the essential aspects of material breakdown. In addition, the deep understanding of fracture processes they provide has served as a starting point for more complex models r13 ; r14 ; r15 ; r19a ; r20 ; r21 . In FBMs the material is made up of a set of parallel fibers, each having a statistically distributed threshold strength or life-time. Besides the classical static FBM one can also introduce a dynamic version r18 . In the static FBM, the failure process is simulated according to a quasi static loading of the material in which the output of the simulation is the ultimate strength of the material, i.e, the maximum load above which it breaks down. In the dynamic FBM, a constant load is maintained on the system and the fibers break by fatigue after some time.

After a fiber breaking, the load acting on it is redistributed among the intact fibers according to the stress redistribution rule associated with the particular load-transfer model. There are two standard load transfer models comprising the FBM and they correspond to the extreme limits of stress redistribution. The global load sharing (GLS) redistributes the load of a failed fiber equally among the active fibers remaining in the system. It is known as the global fiber bundle model (GFBM), and assumes long-range interaction among the fibers which makes it a mean-field approximation that can be solved analytically r11 ; r12 . On the other hand, the local load sharing (LLS) redistributes the load of a failed fiber among the intact fibers that are nearest neighbors to the failed ones, and is thus known as the local fiber bundle model (LFBM). This assumes short-range interaction among the fibers and it has not, in general, been solved analytically. However, in actual heterogeneous materials the stress redistribution is expected to fall in between the two regimes of load transfer. Very recently, some of us have proposed a load transfer scheme with variable range of interaction among the fibers, which interpolates between the two extreme cases. Since most of the physics of the fracture problem is hidden in the stress redistribution, considering within the same model the possibility of varying the range of interaction is a relatively greater improvement toward a better understanding of the fracture problem.

Motivated by the results obtained for the static setting, we have studied a stochastic dynamic fiber bundle model, in which fibers break by fatigue over a period of time, such that the range of interaction among fibers is variable through an adjustable stress transfer function. We have observed a crossover from mean-field to short-range behavior for the macroscopic quantities describing the fracture dynamics as in the static case. In addition, a more detailed inspection of the failure process revealed that the microscopic behavior of the system is also different as the range of interaction varies. Finally, we have also studied the effect of heterogeneity on the failure process from macroscopic as well as microscopic perspectives. The rest of the paper continues with the description of the stochastic model in the next section. Section III is devoted to study the lifetime of the bundle and the rate at which fibers break when the range of interaction changes. The damage spreading for several load transfer modalities is addressed in section IV while the last section is devoted to conclusions rounding off the paper.

## II The stochastic model.

We assign each fiber to the sites of a square lattice, with periodic boundary conditions. Assuming elastic interaction among the fibers, we may state that the stress redistribution obeys a power law r21 ,

$$\sigma _{inc}\frac{1}{r_{ij}^\gamma },$$
(1)
where $`\sigma _{inc}`$ is the load increment on an intact fiber $`i`$ at a distance $`r_{ij}`$ from the failed fiber $`j`$, and $`\gamma `$ is a variable parameter that controls the effective range of interaction among the fibers. $`\gamma `$=0 corresponds to GLS since the additional load on each intact fiber as a result of a fiber breaking is the same irrespective of its distance from the broken fiber. On the other hand, $`\gamma =\mathrm{}`$ corresponds to LLS since in this limit only the nearest neighbors get equal portions of the load of a failed fiber, and $`\sigma _{inc}`$=0 for $`r>1`$. Assuming that there is no dissipation of the load of a failed fiber, Eq. (1) leads us to the normalized stress transfer function,

$$S(\gamma ,r_{ij})=\frac{1}{r_{ij}^\gamma }\frac{1}{_{iA}r_{ij}^\gamma }$$
(2)
$`r_{ij}`$ being the distance between an active fiber $`i`$, with coordinates $`(x_i,y_i)`$, and a failed fiber $`j`$, with coordinates $`(x_j,y_j)`$, and $`A`$ denotes the set of intact fibers.

One can consider two equivalent approaches to the dynamic FBM r22 ; r23 ; r24 . The time elapsed until the final collapse of the system is the lifetime or time to failure of the bundle. In the first approach r22  the lifetime of each element is an independent identically distributed random variable, i.e, each fiber has a different random lifetime, taken from the same statistical distribution (the Weibull distribution is a good empirical distribution in materials science) and each one breaks if the time elapsed exceeds its individual lifetime. This is a quenched model of fracture where the disorder is introduced once for all at the beginning of the process and thus the growth mechanism is deterministic. In this version, the effect of the increase in stress for a particular fiber $`i`$ due to the redistribution of load from failed fibers is the reduction of its initially assigned lifetime $`t_i`$ to a new lifetime $`\tau _i`$ given by

$$t_i=_0^{\tau _i}\left(\frac{\sigma _i(t)}{\sigma _0}\right)^\rho 𝑑t,$$
(3)
where $`\sigma _0`$ is the initial stress on fiber $`i`$ at $`t=0`$ and $`\rho `$ is the Weibull index, $`2\rho 50`$. $`\rho `$ gives the degree of heterogeneity of the system; as $`\rho `$ increases the system becomes more homogeneous. In this model the next element to break is exactly the one with the lowest lifetime at time $`t`$.

In the probabilistic approach r23 ; r24 , it is considered that in a time step all the intact elements have the same mean time to failure. The element that breaks in the time interval between two successive failures is selected by chance and thus the fiber whose probability of breaking (a function of the load it bears) is the largest is more likely to fail. The probabilistic approach is an example of the so called annealed disorder since the algorithm is stochastic and thus randomness is uncorrelated in time. Thus, we start at time $`t=0`$ with $`N`$ fibers loaded with an initial common stress of $`\sigma _i(t)=\sigma _i(0)=1`$. The mean time interval $`\delta `$ for an individual element to break by fatigue is,

$$\delta =\frac{1}{_{iA}\sigma _i^\rho (t)},$$
(4)
where $`A`$ is the set of intact fibers and $`t`$ is the time elapsed up to $`k1`$ breakings, i.e., $`t=_{p=1}^{k1}\delta _p`$. In the first time interval where $`\sigma _i=1`$ $`i`$, with all the $`N`$ elements active, $`\delta _1=\frac{1}{N\sigma _i^\rho }`$ so that the first fiber breaking is completely random. This changes with time due to stress transfers and fiber failures. The probability of a particular fiber breaking $`j`$ in one sweep of the lattice in the time interval $`\delta _k`$, is given by

$$P_j(t)=\sigma _j^\rho (t)\delta _k.$$
(5)
When a fiber breaks, its load is transfered and as a consequence, the probability of failure of the receptors is increased. In this sense there are no weak or strong fibers, all the fibers are equivalent but carrying different loads (except for mean field approaches) and thus with different breaking probabilities. It is worth stressing that both approaches are equivalent. This can be intuitively understood by noting that in both models the load history plays a key role. For the quenched setting the fiber that breaks is that with the lowest lifetime which in its turn depends on the load history. Since the individual times to failure are reduced each time a fiber receives load from a failed element, the more stressed a fiber is, the more likely its lifetime is the lowest. On the other hand, for the annealed version the magnitude that is modified by the load redistribution is the probability of breaking and the load history affects the failure process just as in the deterministic version: the more stressed a fiber is, the more likely it breaks. Additionally, we note that this algorithm is the same as the one used in polymer failure pol  and in describing dielectric breakdown db  with the main distinction that here the broken fibers need not to be connected to the single growing cluster as in dielectric breakdown.

Now, the load born by the fiber that has just failed is redistributed according to Eq. (1) such that in a time interval $`\mathrm{\Delta }t=t_kt_{k1}=_{p=1}^k\delta _p_{p=1}^{k1}\delta _p=\delta _k`$, the load increase on the active fiber $`i`$ reads as

$$\sigma _i(t_k)=\sigma _j(t_{k1})S(\gamma ,r_{ij})+\sigma _i(t_{k1}),$$
(6)
where $`\sigma _j`$ is the load of the element that has failed in the time interval $`\mathrm{\Delta }t`$ after $`k`$ fiber failures. After the redistribution of load the rupture process continues by applying again Eq. (5) that will point to the next fiber to break. The lifetime of the material is then given by the sum of all the $`N`$ $`\delta `$s, $`T_f=_{i=1}^N\delta _i`$.

The complete analytic solution to the dynamic FBM as defined above in the GLS ($`\gamma =0`$) case is feasible. The lifetime of the material can be computed as

$$T_f=_1^N\frac{(N+1x)^{\rho 1}}{N^\rho }𝑑x=\frac{1}{\rho }\left[1\frac{1}{N^\rho }\right]$$
(7)
that includes also the dependence of the lifetime on the system size $`N`$ giving the mean field result $`1/\rho `$ in the thermodynamic limit.

We should note that there is no avalanche here, unlike to the static case r21 , since there is no external driving on the system. Once the fibers start breaking by fatigue, they continue breaking with time until the final collapse of the system, but within a time interval $`\delta _k`$, only one fiber may break in a single sweep of the lattice.

We note that both Eqs. (3) and (5) assume a power-law breaking rule $`\kappa _p(\sigma )=\nu _0(\sigma /\sigma _0)^\rho `$. This stems from the former assumption that the lifetimes satisfy the Weibull distribution. An alternative breaking rule could be and exponential hazard rate of the form $`\kappa _e(\sigma )=\varphi \mathrm{exp}(\eta (\sigma /\sigma _0))`$ mainly used for thermally activated processes. On the other hand, $`\kappa _p`$ has the same functional form as the Charles power law that describes the subcritical crack growth induced by stress corrosion in geological materials at constant temperature: $`v=AK_I^n`$, where $`v`$ is the crack velocity and $`K_I`$ is the stress intensity factor for mode $`I`$ fracture. Sometimes, $`n`$ is known as the stress corrosion index. Moreover, we emphasize that Eq. (4) is not a real measure of time so that a quantitative comparison between the results here shown and those from experiments is not feasible. Additionally, Eq. (4) is a simple form one can consider for the probability of breaking. In principle, one could also include more realistic rules. Nevertheless, as we shall see, this simple model is very rich and might help understand physical effects present in real materials.

## III Lifetime of the bundle

We simulate the failure process by large scale numerical simulations. We first explore the behavior of the lifetime as a function of the stress transfer range for a fixed level of heterogeneity ($`\rho =2`$). Then we vary the Weibull index such that we get a more homogeneous system with smaller times to failure.

The results obtained for the lifetime are depicted in Fig. 1 for different values of the range of interaction among the fibers and several system sizes up to $`N=6400`$ fibers. A crossover from mean field behavior to a regime dominated by the short range interactions among the fibers is clearly appreciated. Furthermore, the critical value of the parameter $`\gamma _c`$ defines a region where the results for global load sharing models hold beyond $`\gamma =0`$, the true mean-field regime for which the load of a failed element is shared equally among the surviving elements. For $`\gamma \gamma _c`$ the material behaves macroscopically as for $`\gamma =0`$, that is, the lifetime is independent on the system size and does not depend on the actual value of the exponent $`\gamma `$. It is not a simple numerical task to determine accurately the exact value of the transition point due to the stochastic nature of the model and the fluctuations of the lifetime of the system. In fact, the time to failure of the bundle follows a Gaussian distribution in what concerns its frequency distribution. The width of the distribution depends on the level of heterogeneity (controlled by the Weibull index) that in turn also influences the lifetimes that become shorter as we move to high levels of homogeneity. Within this numerical uncertainty we have found that $`\gamma _c2`$. Interestingly, the same value was found to characterize the transition of the bundle’s ultimate strength from long to short range behavior in the static case of the model r21 .

The influence of the disorder on the crossover behavior can be studied by changing the value of $`\rho `$. Fig. 2 shows the time to failure of the material as a function of the effective range of interaction for several system sizes. We observe that the transition is still present but the range where the mean field regime applies is reduced. In particular, as the system gets more homogeneous $`\gamma _c`$ shifts leftwards to smaller values. Additionally, the true local load sharing regime appears to be also slightly shifted to the left. When the range of interaction $`\gamma `$ falls below a second transition point $`\gamma _{c2}`$, the lifetime of the bundle becomes again independent on the effective interaction among the fibers but it is still size dependent. This later behavior can be easily understood by noting that for local load sharing schemes the time to failure of the system in the thermodynamic limit is zero. We have checked that in our model this is actually verified, although the drop of the lifetime as the system size is increased is slow. On the other hand, wherever the GLS regime arises, the time to failure of the system is not size dependent for large $`N`$ (see Eq. (7).)

Another way to characterize the evolution of the fracture process is to inspect the rate at which fibers fail. We expect two different asymptotic regimes. For long range interaction, i.e., below $`\gamma _c`$ the system should behave in a mean field manner. This means that damage is gradually accumulated in the material up to a point in which the load is too high as to be carried by the remaining fibers. Only at this stage of the damage process the rate of fiber failures will speed up owing to the small values of the very last $`\delta `$s. On the contrary, in the region where short range interaction prevails, the system does not accumulate damage uniformly. In this case there appear regions within the material in which stress enhancement takes place making the fibers along the crack tips to support much more load than other active fibers placed far from the clusters of broken fibers.

Accordingly, the $`\delta `$s would be modified and there would be more breakings for the same time interval. In other words, the breakdown of the material occurs suddenly for the very localized regime where about $`50\%`$ of the fibers breaks in a time interval of the order of $`0.1T_f`$. In Fig. 3 we have represented the evolution with time of the number of broken fibers $`N_f`$ for different values of $`\gamma `$. Two distinct groups of curves corresponding to the extreme cases can be clearly seen. For intermediate values of the effective range of interaction, the behavior is more like the case of long range interaction and may correspond to other load sharing schemes such as the hierarchical fiber bundle model r18 ; r25 .

Consider that the breaking rate of the bundle is defined as

$$r(t)=\frac{dN_f(t)}{dt},$$
(8)
with $`N_f(0)=0`$ and $`N_f(T_f)=N`$. Upon approaching the complete failure, the breaking rate scales with the lifetime of the material as r26

$$r(t)(T_ft)^\zeta ,$$
(9)
where the exponent $`\zeta `$ depends on both the range of interaction and the Weibull index. However, we can again identify two limiting groups of curves for the same value of $`\rho `$. Figure 4 shows the rate of fiber breaking for several load transfer ranges and a Weibull index $`\rho =2`$. These results confirm the behavior observed for the evolution of the number of broken fibers, namely that there is a sharp increase in the failure rate when approaching the lifetime for the case where short-range interaction dominates the damage spreading. The fit to the curves gives $`\zeta 1/2`$ for $`\gamma \gamma _c`$ and $`\zeta 2/3`$ when $`\gamma `$ is in the range where the effective interaction among the fibers can be assumed to be very localized. Note that the curves for intermediate values of $`\gamma `$ and for the GLS regime have been shifted for the sake of clarity. The numerical results are not very smooth because of fluctuations but the general trend of $`r(t)`$ confirms the validity of relation (9). As to the dependence of the above results on the heterogeneity level we have observed that the less heterogeneous the material is, the sharper the failure acceleration is in all cases.

The inset in Fig. 4 shows the breaking rate for the case of long range interaction and $`\rho =10`$. The higher value of $`\zeta 0.86`$ indicates that much of the fibers break in a very small time interval close to the lifetime of the material. As the range of interaction gets more localized, the exponent $`\zeta `$ approaches unity.

## IV Crack growth

A further characterization of what is going on in the fracture process can be carried out by focusing on the properties of the clusters of broken fibers. Specifically, we have monitored the growth of the cracks inside the bundle. Cracks are defined as connected clusters or regions of broken fibers. Here, we consider a coordination number $`q=8`$, that is, each fiber has 8 neighbors. Similar results are obtained if we take into account only nearest neighbors ($`q=4`$). At the very initial stages, regardless of the range of interaction among the fibers, the failure of fibers can be assumed to be random, that is, the initial cracks are randomly nucleated inside the material. This situation changes with time for different load transfer schemes. By studying the growth of the largest crack area $`C_m`$ at each time step, one could distinguish the different mechanisms leading to the rupture of the material as the range of interaction varies. It is worth noting that this is just a way that allows to discern between different ranges of interaction and levels of heterogeneity. For instance, one can consider instead the linear size of the largest crack and study the fractal dimension of the crack distribution for different $`\rho `$.

Figure 5 shows a microscopic aspect of the material breakdown with time for the two limiting cases of load redistribution and a system consisting of $`N=900`$ fibers ($`\rho =2`$). Initially, it can be seen that cracks are randomly nucleated in the material. At later times the individual microcracks coalesce thereby causing a jump in the largest crack area. For long range interaction, the nucleation of cracks continues to be random because all the fibers carry the same load and thus they break by chance. Therefore the largest broken cluster does not grow linearly with the number of broken fibers. In this case, there are isolated cracks inside the material that grow essentially by coalescence when they meet each other. This is the reason why sudden jumps in the area of the largest crack are observed in the intermediate stages of the damage spreading. At the end of the process, the material has accumulated many of these cracks giving rise to the linear crack growth shown in the inset.

For localized range of interaction, the mechanism of damage spreading is radically different. Again, at small times, the cracks are randomly nucleated inside the material. However, as time goes on and more fibers get broken, the load is redistributed to the fibers located along the crack tips provoking the accumulation of stress in these fibers and the appearance of regions where fibers bear a huge amount of load. It is thus expected that the newly broken fibers add to already existing cracks and that a dominant crack appears. From this perspective, the largest crack area should scale linearly with the number of broken fibers, i.e., $`C_mN_f`$. This is indeed the case for $`\gamma =20`$ as represented in Fig. 5. The straight line is a linear fit to a time window in which more than 150 fibers have been broken. The value $`0.90`$ of the slope confirms the above picture. Note that in this case the number of coalescence events is smaller than for the GLS regime and that after one of such events the linear growth of the largest crack is recovered. We also note that, up to the intermediate stages of the rupture process, the largest crack for a given number of failed fibers is much higher in the localized case than the global case indicating the formation of a (few) dominant crack(s) in the later case. Additionally, at the end of the process there are no differences between the two extreme load transfer schemes since more than a half of the material is already broken and is very unlikely to find a region where isolated cracks can be formed and grow. Thus, at the final stage each additional breaking event occurs at the crack tips of existing single dominant cracks. Nevertheless, as stated in the previous section the rate of fiber failures is quite different in both asymptotic regimes. We shall note here that the same behavior as for the mean field regime is observed for any value of $`\gamma `$ provided it is below the transition point $`\gamma _c=2`$.

Figure 6 further substantiates our previous arguments by showing how the largest crack area varies as a function of the number of broken fibers for several levels of heterogeneity. While for $`\gamma =0`$ the picture is always the same (inset), for $`\gamma =10`$, that is, in the localized regime, the time taken for cracks to become dominant decreases with increasing $`\rho `$. Furthermore, when the system is very homogeneous ($`\rho =10`$) and local interactions prevail, a dominant crack which grows until the material collapses is formed almost instantaneously confirming that the mechanism of rupture and crack growth for homogeneous materials is radically different from that of heterogeneous systems. Nevertheless, for the global load case, the change in system’s homogeneity does not alter the way dominant cracks appear and grow (see the inset, where no changes, apart form statistical fluctuations, can be observed). The reason why this happens is given by the way the system gets broken. Equation (5) tell us that the more stressed a fiber is, the more likely it breaks. This always holds except for the global load sharing case, where the fibers share the same amount of load and thus all of them have the same probability to break. As the system is more homogeneous (larger $`\rho `$), for local load sharing schemes, the probability of breaking for the same load is higher so that the appearance of a dominant crack is enhanced. Therefore, for long range interaction there is no correlated crack growth while for short range regimes this is precisely the dominant mechanism since the first stages of the damage spreading.

## V Conclusions

We have extended the fiber bundle model with variable range of interaction between fibers to the dynamic setting. As for the static version, two very different regimes are identified as the exponent of the stress transfer function varies. The lifetime of the material for $`\gamma <\gamma _c`$ does not depend on both the system size and $`\gamma `$ qualifying for a mean field behavior. On the contrary, for the short range regime, the time to failure of the system systematically decreases when increasing the size of the bundle. The analysis performed also showed that the exact value of the crossover point depends on the level of heterogeneity of the system. Besides, we investigated how the material approaches its point of total breakdown at both sides of the crossover point. The crossover from one regime to the other also influences the behavior of the rate at which fibers break, explicitly manifested in a power law divergence as $`T_f`$ is approached, but with an exponent that depends on the range of interaction. This result is relevant from a practical point of view since for the localized regime the acceleration of the failure process takes place at the very final stages of the rupture process. Although fibers break by fatigue, one by one, they do break in very different time intervals according to the range of interaction. In this sense, a global load sharing regime is safer, since we get more warnings before the material breakdown. On the other hand, the precursory activity when the range of interaction gets localized is almost absent leading to a sudden breakdown of the bundle in a very short time interval.

The numerical exploration of the damage spreading mechanisms under different load transfer regimes further supported the results obtained for the lifetime of the system. Regardless of the range of interaction, the breaking of fibers is a completely random process at the initial stages of the failure process. After some time, the mechanism of failure propagation radically changes when the exponent $`\gamma `$ varies. In the limiting case of global load sharing there is no correlated crack growth in the system, whereas for the short range regime the damage spreading is driven by a dominant crack, and thus, the crack growth is strongly correlated with high stress concentration at the fibers located along the perimeter of the dominant cluster of broken fibers.

These differences are clearly appreciated in Fig. 7, where we represent snapshots of the lattice when nearly $`40\%`$ of the system is broken for several values of $`\gamma `$ and $`\rho `$. For the long range interaction limit the material’s level of heterogeneity does not influence the random nucleation and growth of cracks and there are no clearly distinguishable dominant cracks. This continues to be so until coalescence drives the further breaking of the material. On the contrary, for localized regimes (left column in Fig. 7), when the system gets more homogeneous the dominant crack appears at early stages and damage spreading is strongly correlated resembling a single crack growth mechanism typical of homogeneous materials. When the bundle is heterogeneous, crack growth is still correlated, but in this case we can identify more than one large and dominant crack. Finally, our results suggest that actually there are only two limiting cases relevant to experiments. The one in which mean field assumptions apply could be of great importance since this will allow the extension of known analytic results to ranges of interaction beyond $`\gamma =0`$.

## VI Acknowledgments

One of the authors (Y. M. ) would like to thank A. F. Pacheco and J. B. Gómez for valuable comments on this work. O. E. Y. thanks the ICTP and UNESCO for their financial support and hospitality. Y. M. acknowledges financial support from the Ministerio de Educación Cultura y Deportes (Spain) and of the Spanish DGICYT Project BFM2002-01798. F. K. acknowledges financial support of the Bólyai János Fellowship of the Hungarian Academy of Sciences and of the Research Contracts FKFP 0118/2001 and T037212. This work was partially supported by the project SFB381, and by the NATO grant PST.CLG.977311.


format: text

text:
# Aspects of the Cosmic Microwave Background Dipole*footnote **footnote *Dedicated to the memory of David Wilkinson.

## I Introduction

The primary aim of NASA’s recently launched Microwave Anisotropy Probe (MAP)  and the European Space Agency’s Planck satellite  will be to measure small-scale fluctuations in the cosmic microwave background (CMB) in order to better determine cosmological parameters  and test inflationary cosmology , and thus improve upon the already remarkable results from recent balloon and ground experiments . However, the most salient feature in the CMB is not this small-scale structure; it is the much more prominent dipole  of amplitude $`\mathrm{\Delta }T=3.365\pm 0.027`$ mK.

The simplest explanation for the dipole is a Local Group velocity $`v=627\pm 22`$ km sec-1 toward $`(l,b)=(276^{}\pm 3^{},33^{}\pm 3^{})`$ with respect to the CMB rest frame. Some measurements of galaxy velocities on large scales suggest independently, although far less precisely, a similar value . It has thus been generally accepted that the CMB dipole is due to our peculiar velocity, but theorists have occasionally speculated that the dipole might be due, at least in part, to an intrinsic temperature fluctuation .

If we had sufficiently precise multi-frequency observations to measure the frequency spectrum at each point on the sky, it would alway be consistent with a blackbody spectrum, although with an angle-dependent temperature. In practice, the observation frequencies and sensitivities are limited. The temperature at each point on the sky is thus determined by measuring the intensity fluctuation, and then converting it to a temperature fluctuation, generally assuming that the temperature fluctuation is sufficiently small that the intensity fluctuation can be related to the temperature fluctuation by the first term in a Taylor expansion of the Planck function, the relation between the temperature and intensity. With prior experiments, the temperature fluctuation has always been sufficiently small that this approximation has been warranted.

In this paper we point out that with the observed dipole amplitude and Planck’s (and possibly MAP’s) sensitivity, this will no longer be a good approximation—the second-order term in the Taylor expansion will be non-neglible. To quadratic order in the dipole amplitude, there is an intensity quadrupole with a frequency dependence given by the second derivative of a Planck function. Since this frequency dependence differs from that of the first-order term, this dipole-induced intensity quadrupole can be disentangled from the temperature quadrupole, even after taking into account a number of foregrounds, as we show below. This higher-order effect will provide a robust pre-determined signal that Planck and possibly MAP should be able to detect. It can be used to calibrate large-angle-fluctuation measurements and.or benchmark multi-frequency foreground-subtraction techniques. It may thus be a useful addition to the astrophysical point sources and annual modulation of the dipole that have until now served as calibration sources for CMB experiments.

In the next Section we discuss the dipole-induced intensity quadrupole, which we refer to simply as the “dipole quadrupole” and abbreviate DQ. In Section III, we then show that the DQ should be detectable by Planck and distinguishable from the CMB quadrupole, even after subtracting several foregrounds. Our analysis in Sections II and III assumes that the dipole is due entirely to the velocity of the solar system with respect to the CMB rest frame. To indicate the magnitude and detectability of of the DQ, we evaluate how well the magnitude and orientation of the peculiar velocity could be determined—assuming the dipole were due to a peculiar velocity—from the DQ alone. In Section IV we show that the same DQ arises even if the CMB dipole is an intrinsic temperature fluctuation, rather than the result of a peculiar velocity. Section V provides some closing remarks.

## II The Dipole-Induced Intensity Quadrupole

The specific intensity of a blackbody of temperature $`T_{\mathrm{CMB}}`$ is

$$I_\nu ^{}^{}=C\frac{x^3}{e^x^{}1},$$
(1)
where $`C=2(k_BT_{\mathrm{CMB}})^3(hc)^2`$ is a constant, $`x^{}=h\nu ^{}(k_BT_{\mathrm{CMB}})^1`$, $`\nu ^{}`$ is the photon frequency, and $`h`$, $`k_B`$, and $`c`$ are, respectively, the Planck and Boltzmann constants and the speed of light. The photon frequency measured by an observer moving with velocity $`v`$ relative to the blackbody is $`\nu =\gamma \nu ^{}(1+\beta \mu )`$, where $`\beta =v/c`$, $`\gamma =(1\beta ^2)^{1/2}`$, and $`\mu `$ is the cosine of the angle between the velocity and the photon direction. The specific intensity transforms as $`I_\nu ^{}^{}/\nu ^3=I_\nu /\nu ^3`$, and thus the observer-frame specific intensity is

$$I_\nu =C\frac{x^3}{e^{x\gamma (1+\beta \mu )}1},$$
(2)
where $`x=h\nu (k_BT_{\mathrm{CMB}})^1`$. Expanding in $`\beta `$,

$$I_\nu =C\frac{x^3}{e^x1}\left\{\left[1+O(\beta ^2)\right]f(x)\beta \mu +f(x)g(x)\beta ^2(\mu ^21/3)+\mathrm{}\right\}.$$
(3)
where $`f(x)xe^x(e^x1)`$ and $`g(x)(x/2)\mathrm{coth}(x/2)`$.

The term linear in $`\beta `$ is the dipole, with the appropriate frequency spectrum $`f(x)`$ for the dipole , which is also more generally the frequency spectrum for a thermal fluctuation. However, to order $`\beta ^2`$, there is another term, the DQ, with a frequency spectrum that differs from that from small thermal fluctuations. With COBE-DMR’s sensitivity and frequencies (in the Rayleigh-Jeans regime), the frequency dependence was negligible. However, with the improved sensitivity and expanded frequency range of new satellite experiments, the frequency dependence will become detectable. This frequency dependence was noted in Refs. , but seems to have since been overlooked. Refs.  identified the effect (sometimes referred to as the “kinematic quadrupole”) in calculations of CMB polarization induced by reionization, but they did not consider the CMB dipole. Ref.  discusses peculiar-velocity effects on the CMB, but only on small-scale fluctuations.

## III Detectability of the DQ

To assess the detectability of the DQ, we determine how well the amplitude and direction of the solar-system velocity can be determined with MAP and Planck under the assumption that the velocity is to be reconstructed entirely from the DQ.

Since the frequency dependences of the CMB quadrupole and the DQ differ, they can be distinguished if the quadrupole moments of the intensity are determined in several frequencies. However, the CMB quadrupole and the DQ will also have to be distinguished from the quadrupole moments due to unsubtracted, or imperfectly substracted, foregrounds, each with its own frequency dependence.

In what follows we discuss the brightness temperature (which is proportional to the intensity) rather than the intensity so that we work with units that are familiar in CMB studies. If the brightness temperature $`T_\nu (\widehat{𝐧})`$ is measured in several frequencies $`\nu `$ as a function of position $`\widehat{𝐧}`$ over the entire sky, then the five (for $`m=2,1,0,1,2`$) quadrupole components $`T_{m,\nu }`$ can be constructed from $`T_{m,\nu }=Y_{2m}(\widehat{𝐧})T_\nu (\widehat{𝐧})𝑑\widehat{𝐧}`$, where $`Y_{lm}(\widehat{𝐧})`$ are spherical harmonics. The possible contributions to these $`T_{m,\nu }`$ from several unsubtracted foregrounds (dust, synchrotron radiation, bremsstrahlung, the far infrared background, and the Sunyaev-Zeldovich effect), which will have frequency dependences shown in Fig. 1, must also be considered. The quadrupole moment at each frequency is then the sum of the contributions from all of these sources: $`T_{m,\nu }=_iT_m^if_i(\nu )`$, where $`i=\{\mathrm{CMB},\mathrm{DQ},\mathrm{dust},\mathrm{synch},\mathrm{brem},\mathrm{firb},\mathrm{SZ}\}`$, and $`f_i(\nu )`$ is the frequency dependence of source $`i`$, referenced to the CMB frequency dependence. Thus, $`f_{\mathrm{CMB}}(\nu )=1`$; $`f_{\mathrm{DQ}}(\nu )=(x/2)\mathrm{coth}(x/2)`$; $`f_{\mathrm{dust}}(\nu )=x^2e^{x_dx}(e^x1)^2(e^{x_d}1)^2`$ (where $`x_d=h\nu (k_BT_{\mathrm{dust}})^1`$ and we have taken a thermal dust spectrum with $`T_{\mathrm{dust}}=20`$ K and emissivity index $`\alpha =2`$); $`f_{\mathrm{synch}}(\nu )=x^{4.8}(e^x1)^2e^x`$; $`f_{\mathrm{brem}}(\nu )=x^{4.16}(e^x1)^2e^x`$; $`f_{\mathrm{firb}}(\nu )=x^{1.64}e^{x_fx}(e^x1)^2(e^{x_f}1)^2`$ (where $`x_f=h\nu (k_BT_{\mathrm{firb}})^1`$ and we have taken a thermal FIRB spectrum with $`T_{\mathrm{firb}}=18.5`$ K and emissivity index $`\alpha =0.64`$ as a fit to data ); and $`f_{\mathrm{SZ}}(\nu )=2f_{\mathrm{CMB}}(\nu )f_{\mathrm{DQ}}(\nu )`$. For a review of foreground properties see Ref. .

For each of the five components $`m`$, we will have five (for the MAP frequencies of 22, 30, 40, 60, and 90 GHz) or nine (for the Planck frequencies of 30, 44, 65, 100, 143, 217, 353, 545, and 857 GHz) data points $`T_{m,\nu }`$ that we model with a vector of parameters $`s_i=T_m^i`$ again for $`i=\{\mathrm{CMB},\mathrm{DQ},\mathrm{dust},\mathrm{synch},\mathrm{brem},\mathrm{firb},\mathrm{SZ}\}`$. Following Ref. , the standard error ($`1\sigma `$) to the quadrupole component $`T_m^{\mathrm{DQ}}`$ will be $`\sigma _m^{\mathrm{DQ}}=[\alpha ^{1/2}]_{\mathrm{DQ},\mathrm{DQ}}`$, where

$$\alpha _{ij}=\underset{\nu }{}\left[\frac{T_{m,\nu }}{s_i}\frac{T_{m,\nu }}{s_j}\right]\frac{1}{\sigma _{m,\nu }^2}=\underset{\nu }{}\frac{f_i(\nu )f_j(\nu )}{\sigma _{m,\nu }^2},$$
(4)
is the covariance (Fisher) matrix. Here, $`\sigma _{m,\nu }^2`$ is the variance to the quadrupole from instrumental noise at frequency $`\nu `$; we use 4-year values for MAP from Ref.  and 2-year values for Planck from Ref. .

Once the DQ components $`T_m^{\mathrm{DQ}}`$ have been determined, the peculiar-velocity components can be found. Expanding the angular dependence, $`(\mu ^21/3)`$, where $`\mu `$ is the cosine of the angle between the velocity and the photon direction, in spherical harmonics gives us $`(T_{m=0}^{\mathrm{DQ}}/T_{\mathrm{CMB}})=\sqrt{4\pi /45}(\beta _x^2+\beta _y^22\beta _z^2)`$, $`\sqrt{2}(\mathrm{Re}T_{m=1}^{\mathrm{DQ}}/T_{\mathrm{CMB}})=\sqrt{16\pi /15}\beta _x\beta _z`$, $`\sqrt{2}(\mathrm{Im}T_{m=1}^{\mathrm{DQ}}/T_{\mathrm{CMB}})=\sqrt{16\pi /15}\beta _y\beta _z`$, $`\sqrt{2}(\mathrm{Re}T_{m=2}^{\mathrm{DQ}}/T_{\mathrm{CMB}})=\sqrt{4\pi /15}(\beta _x^2\beta _y^2)`$, and $`\sqrt{2}(\mathrm{Im}T_{m=2}^{\mathrm{DQ}}/T_{\mathrm{CMB}})=\sqrt{16\pi /15}\beta _x\beta _y`$. Let us suppose that these moments have been measured precisely. Then the equations for $`T_{m=0}^{\mathrm{DQ}}(a)`$, $`\mathrm{Re}T_{m=2}^{\mathrm{DQ}}(b)`$, and $`\mathrm{Im}T_{m=2}^{\mathrm{DQ}}(c)`$ can be inverted (after introducing the shorthand $`a`$, $`b`$, and $`c`$) to give us the components $`\beta _x^2=\sqrt{15(16\pi )^1}(\sqrt{b^2+c^2}+b)`$, $`\beta _y^2=\sqrt{15(16\pi )^1}(\sqrt{b^2+c^2}b)`$, and $`\beta _z^2=\sqrt{15(16\pi )^1}(\sqrt{b^2+c^2}\sqrt{3}a)`$. Since we have now determined the squares of three components of the velocity, we are still left with a residual eightfold degeneracy in the velocity. However, the signs of the components $`\mathrm{Re}T_{m=1}^{\mathrm{DQ}}`$ and $`\mathrm{Im}T_{m=1}^{\mathrm{DQ}}`$, which are proportional, respectively, to $`\beta _x\beta _z`$ and $`\beta _y\beta _z`$, can be used to break two of these degeneracies, leaving us with only a twofold degeneracy, the sign of the velocity, which is undetermined by the quadrupole.

When the measurements are done, we will then want to check whether the peculiar velocity induced from the DQ is consistent with that inferred from the dipole itself. If the data analysis and our understanding of the instrument are reliable, then the two velocities should be consistent. To quantify the degree to which MAP and Planck can test this consistency we calculate, assuming $`v=0`$, the expectation value of $`\chi ^2`$ for the hypothesis that $`v=627`$ km sec-1. If, e.g., $`\chi ^2=9`$ then the $`v=0`$ hypothesis can be distinguished from the $`v=627`$ km sec-1 hypothesis at $`3\sigma `$. We find, for arbitrary $`v`$, $`\chi ^2=9(4.2\mu \mathrm{K}/\sigma _m^{\mathrm{DQ}})^2\left[v/(627\mathrm{km}\mathrm{sec}^1)\right]^4`$. From this $`\chi ^2`$ we also derive the minimum $`v`$ distinguishable at $`3\sigma `$ from $`v=0`$ to be $`v_{\mathrm{min}}=305`$ km sec-1 $`(\sigma _m^{\mathrm{DQ}}/\mu \mathrm{K})^{1/2}`$.

As another indicator of the precision with which the DQ can be reconstructed, we now determine the precision with which the magnitude and orientation of the velocity can be inferred from the DQ alone. There will be a measurement error of $`\sigma _m^{\mathrm{DQ}}`$ to each of the $`T_m^{\mathrm{DQ}}`$, as described above, and the three velocity components, $`v_x`$, $`v_y`$, and $`v_z`$, will be fit to all five of these quadrupole components. To estimate the errors on these components we once again calculate a Fisher matrix, this time for these five data points and the parameters $`v_x`$, $`v_y`$, and $`v_z`$. We then invert the Fisher matrix and choose the $`z`$ axis along the inferred velocity so that $`v=v_z`$. We then find that the measurement error (at $`3\sigma `$) to the velocity is $`\mathrm{\Delta }v=3\sqrt{45(64\pi )^1}\sigma _m^{\mathrm{DQ}}c^2(v_zT_{\mathrm{CMB}})^1`$, where $`v`$ is the best-fit velocity. Under the null hypothesis that the peculiar velocity is that inferred from the dipole, we find $`\mathrm{\Delta }v=225(\sigma _m^{\mathrm{DQ}}/\mu \mathrm{K})`$ km sec-1. The Fisher-matrix analysis tells us that the measurement errors to the $`x`$ and $`y`$ components (those perpendicular to the best-fit velocity) are each $`\sqrt{4/3}`$ times that for the $`z`$ component. Thus, under the same null hypothesis, the ($`3\sigma `$) error to the velocity orientation will be $`24^{}(\sigma _m^{\mathrm{DQ}}/\mu \mathrm{K})`$.

Table I shows results of our calculations of $`\chi ^2`$ and $`\mathrm{\Delta }v`$ assuming a variety of combinations of components will be fit to the data. If the dust contribution to MAP must be determined from the MAP data itself, then the smallest detectable peculiar velocity is $`2100`$ km sec-1, too big to be interesting. This is because the highest-frequency channel becomes a dust monitor, and the remaining MAP channels at $`60`$ GHz do not provide enough leverage to disentangle the DQ and the CMB quadrupole. However, if the dust quadrupole can be determined precisely from other observations, then MAP might be able to isolate the DQ to better than $`7\sigma `$. The synchrotron and bremsstrahlung foregrounds are strongest at low frequencies, where the DQ and CMB frequency dependences are similar. Thus, if synchrotron and bremsstrahlung emission are included in the analysis, the lowest-frequency channels become foreground monitors and the ability to separate the DQ is not degraded significantly.

Because of the improved detector sensitivity, and especially the broadened frequency coverage, the outlook for Planck is much better. Even without the highest-frequency channels, which act as dust monitors, there is still a sufficiently broad spectrum of frequencies $`100`$ GHz where the spectral dependences of the DQ and the CMB quadrupole differ the most. Again, the synchrotron and bremsstrahlung foregrounds contribute mostly at low frequencies and thus do not degrade significantly the DQ signal. Quite remarkably, even if we marginalize over a number of uncertain foreground amplitudes, Planck should be able to detect the DQ as long as the velocity is greater than roughly 140 km sec-1. Moreover, Planck should be able to detect a deviation from the dipole-inferred velocity as small as 45 km sec-1 and determine its direction (modulo the sign) to better than $`5^{}`$.

It turns out that the Sunyaev-Zeldovich frequency dependence is a linear combination of the CMB-quadrupole and DQ frequency dependences. Thus, the DQ is formally degenerate with the cosmological and Sunyaev-Zeldovich quadrupoles. However, the quadrupole moment due to the Sunyaev-Zeldovich effect will probably be negligible. Calculations of the Sunyaev-Zeldovich power spectrum near $`l100`$ find $`l^2C_l/(2\pi )10^{12}`$ . Doing shot-noise extrapolation to the quadrupole, we find $`\sigma _m^{\mathrm{SZ}}0.15\mu \mathrm{K}`$, corresponding to a velocity uncertainty (at $`3\sigma `$) $`35`$ km sec-1, which is just a bit below the expected statistical uncertainty indicated in Table I for Planck. The actual SZ number will probably be much smaller, and certainly much more will be known by the time Planck flies. The SZ octupole will be useful in constraining the possible SZ quadrupole contribution.

Our foreground modeling makes several simplifying assumptions. For example, we do not actually know the frequency dependences perfectly and they can vary spatially (due to, e.g., spatial variation in chemical composition of the dust). There may also be components we have not yet considered, such as spinning dust. Further, about 25% of the sky will be lost to a Galactic-plane cut. It is thus possible that Planck may not be able to achieve the tiny velocity errors quoted in the table. However, the values in the final row may be achievable. Although uncertain, the frequency dependences of the foregrounds are all (with the exception of the SZ effect) considerably different than those of the the DQ and the CMB quadrupole. The foregrounds are all likely to have low amplitudes near 100 GHz: $`3\mu `$K for dust, $`1\mu `$K for bremsstrahlung and less for everything else . Thus the DQ is not far below foreground quadrupoles, so multi-frequency foreground subtraction need not be done to better than about 10%. Of course, only the measurements themselves will answer these questions definitively.

## IV An Intrinsic Dipole?

Although a peculiar velocity is the simplest explanation for the dipole, and certainly that most consistent with the prevailing inflationary paradigm, it is also possible that the dipole could be due, at least in part, to an intrinsic temperature dipole produced by a super-horizon entropy perturbation . In this Section, we show that even if the temperature dipole is intrinsic, the same DQ still arises. We also clarify the distinction between an intrinsic CMB temperature fluctuation, the thermal quadrupole induced by a peculiar velocity, and the DQ, which unlike the other two, is not a temperature quadrupole.

A velocity $`v`$ induces a temperature pattern $`T(\theta )=T_{\mathrm{CMB}}(1+\beta \mu )^1\sqrt{1\beta ^2}`$ . To second order in $`\beta ^2`$, this can be written as

$$T(\theta )=T_{\mathrm{CMB}}\{[1+O(\beta ^2)]\beta \mu +\beta ^2(\mu ^21/3)\}.$$
(5)
Thus, the peculiar velocity induces a temperature quadrupole of magnitude $`O(\beta ^2T_{\mathrm{CMB}})`$, in addition to the dipole of amplitude $`O(\beta T_{\mathrm{CMB}})`$. However, the quadrupole in the intensity \[Eq. (3)\] arises from a combination of both the temperature dipole and the temperature quadrupole. To see this, we rewrite Eq. (5) as

$$T(\theta )=T_{\mathrm{CMB}}\{[1+O(\beta ^2)]\beta _1\mu +\beta _2^2(\mu ^21/3)\},$$
(6)
so that we can see where the DQ comes from. Doing so, we find that Eq. (3) becomes

$$I_\nu =C\frac{x^3}{e^x1}\left\{\left[1+O(\beta ^2)\right]f(x)\beta _1\mu +f(x)\left[\beta _2^2+\beta _1^2\left(g(x)1\right)\right](\mu ^21/3)+\mathrm{}\right\}.$$
(7)
We thus see that what we have been calling the DQ consists of two parts: The first is the $`\beta _2^2`$ term, which is due to the temperature quadrupole induced by our peculiar velocity. Since this is an honest-to-goodness temperature quadrupole, it has the frequency dependence of the usual lowest-order thermal fluctuation, and it cannot be distinguished from an intrinsic CMB temperature fluctuation. The second term, proportional to $`\beta _1^2`$, arises from the term in the Taylor expansion of the intensity that is second order in the dipole amplitude. This second-order term has a frequency dependence that differs from the usual lowest-order thermal fluctuation.

It is also clear from Eqs. (6) and (7) that even if the temperature dipole were intrinsic (that is, due to an entropy perturbation), then there would still be a DQ. Although the frequency dependence is $`g(x)1`$ rather than $`g(x)`$, only the part proportional to $`g(x)`$ can be distinguished by multi-frequency maps from a temperature quadrupole. For this reason, the amplitude and orientation of the DQ would be exactly the same as if the dipole were due to a velocity, and the detectability would be exactly as determined above. Likewise, the DQ cannot be used to tell whether the dipole is due to a peculiar velocity or due to an entropy perturbation, as previously suggested .We thank E. Wright for illuminating discussions on this point and for pointing out an error in an earlier draft.

Finally, we mention one last point. Strictly speaking, an intrinsic CMB temperature quadrupole will produce a contribution to the intensity quadrupole with a frequency dependence proportional to $`g(x)`$, just like the DQ. However, if the temperature-quadrupole amplitude is $`\mathrm{\Delta }T`$ ($`10^5`$ for the intrinsic fluctuation and $`10^6`$ for the velocity-induced temperature quadrupole), then this contribution will be of order $`(\mathrm{\Delta }T)^210^{10}`$, much smaller than the $`10^6`$ expected for the DQ.

## V Conclusions

We have shown that the dipole amplitude is sufficiently large that the discrepancy between the exact frequency dependence of a thermal fluctuation and the lowest-order frequency dependence usually assumed will be detectable by Planck and possibly MAP. To second-order in the dipole amplitude, this discrepancy is manifest as an intensity quadrupole that can be distinguished with multi-frequency measurements from an intrinsic temperature quadrupole. This provides a robust pre-determined target for CMB experiments and it may prove to be a useful tool for calibration for forthcoming space experiments.

We thank the Santa Barbara KITP for hospitality. This work was supported at Caltech by NSF AST-0096023, NASA NAG5-9821, and DoE DE-FG03-92-ER40701, at Davis by NASA NAG5-11098, and at the KITP by NSF PHY99-07949.


text:
# Smallest small-world network

Phys. Rev. E 66, 046139 (2002)

## Abstract

Efficiency in passage times is an important issue in designing networks, such as transportation or computer networks. The small-world networks have structures that yield high efficiency, while keeping the network highly clustered. We show that among all networks with the small-world structure, the most efficient ones have a single “center”, from which all shortcuts are connected to uniformly distributed nodes over the network. The networks with several centers and a connected subnetwork of shortcuts are shown to be “almost” as efficient. Genetic-algorithm simulations further support our results.

The small-world network models have received much attention from researchers in various disciplines, since they were introduced by Watts and Strogatz Watts and Strogatz (1998) as models of real networks that lie somewhere between being random and being regular. Small-world networks are characterized by two numbers: the average path length $`L`$ and the clustering coefficient $`C`$. $`L`$, which measures efficiency footnote1  of communication or passage time between nodes, is defined as being the average number of links in the shortest path between a pair of nodes in the network. $`C`$ represents the degree of local order, and is defined as being the probability that two nodes connected to a common node are also connected to each other.

Many real networks are sparse in the sense that the number of links in the network is much less than $`N(N1)/2`$, the number of all possible (bidirectional) links. On one hand, random sparse networks have short average path length (i.e., $`L\mathrm{log}N`$), but they are poorly clustered (i.e., $`C1`$). On the other hand, regular sparse networks are typically highly clustered, but $`L`$ is comparable to $`N`$. (All-to-all networks have $`C=1`$ and $`L=1`$, so they are most efficient, but most expensive in the sense that they have all $`N(N1)/2`$ possible connections and so they are dense rather than sparse.) The small-world network models have advantages of both random and regular sparse networks: they have small $`L`$ for fast communication between nodes, and they have large $`C`$, ensuring sufficient redundancy for high fault tolerance. Many networks in the real world, such as the world-wide web (WWW) Albert et al. (1999), the neural network of *C. elegans* Watts and Strogatz (1998); Watts (1999), collaboration networks of actors Watts and Strogatz (1998); Watts (1999), networks of scientific collaboration Newman (2001), and the metabolic network of *E. coli* Wagner (2001), have been shown to have this property. The models of small-world networks are constructed from a regular lattice by adding a relatively small number of shortcuts at random, where a link between two nodes $`u`$ and $`v`$ is called a *shortcut* if the shortest path length between $`u`$ and $`v`$ in the absence of the link is more than two Watts (1999). The regularity of the underlying lattice ensures high clustering, while the shortcuts reduce the size of $`L`$.

Most work has focused on average properties of such models over different realizations of random shortcut configurations. However, a different point of view is necessary when a network is to be designed to optimize its performance with a restricted number of long-range connections. For example, a transportation network should be designed to have the smallest $`L`$ possible, so as to maximize the ability of the network to transport people efficiently, while keeping a reasonable cost of building the network. The same can be said about communication networks for efficient exchange of information between nodes. We fix the number of shortcuts here and as a result the clustering coefficient $`C`$ for any configuration of shortcuts is approximately as high as that of the underlying lattice. The problem we address in this paper is: *given a number of shortcuts in a small-world network, which configuration of these shortcuts minimizes $`L`$?* footnote2 .

Most random choices of shortcuts result in a suboptimal configuration, since they do not have any special structures or organizations. On the contrary, many real networks have highly structured configurations of shortcuts. For example, in long-range transportation networks, the airline connections between major cities which can be regarded as shortcuts, are far from being random, but they are organized around hubs. Efficient travel involves ground transportation to a nearest airport, then flights through a hub to an airport closest to the destination, and ground transportation again at the end.

In the following, we show that the average path length $`L`$ of a small-world network with a fixed number of shortcuts attains its minimum value when there exists a “center” node, from which all shortcuts are connected to uniformly distributed nodes in the network footnote3 . An example of such a configuration is illustrated in Fig. 1(a).

We also show that if a small-world network has several “centers” and its subnetwork of shortcuts is *connected*, then $`L`$ is almost as small as the minimum value. An example of such configuration is shown in Fig. 1(b). We then derive an explicit formula for the minimum average path length in the case of the small-world network models constructed from a one-dimensional lattice by adding a fixed number of shortcuts. Finally, we verify the results by performing genetic-algorithm simulations for minimizing $`L`$.

Our general argument proceeds as follows. A small-world network is composed of two parts: the underlying network (e.g., a regular lattice) and the subnetwork of shortcuts containing only the shortcuts and their nodes. Let $`m`$ denote the number of shortcuts. First, for $`L`$ to be as short as possible, the subnetwork of shortcuts must be connected. This connectivity is unlikely to happen if the shortcuts are chosen at random, since the network is sparse. Indeed, the probability is less than $`m!/N^{m1}`$, where $`N`$ is the number of nodes in the network. For example, for $`N=1000`$ and $`m=10`$, the probability is smaller than $`10^{22}`$. Having a disconnected component in the subnetwork of shortcuts increases the value of $`L`$. In particular, consider the configuration of shortcuts as shown in Fig. 2(a), where one of the shortcuts in Fig. 1(a) is disconnected from the rest of the subnetwork of shortcuts.

If the shortest path between a pair of nodes involves going from the disconnected shortcut to the rest of the subnetwork, then its length is increased by 2 compared to the path length between the corresponding pair in Fig. 1(a). This increases the average path length $`L`$.

Next, observe that the nodes in the subnetwork of shortcuts must be uniformly distributed over the network. This can be seen by noting that the average length of the shortest path from a node to its nearest shortcut is smallest when these nodes are uniformly distributed.

Finally, among all possible configurations of connected subnetworks of shortcuts with uniformly distributed nodes, ones with a single center involve the largest number of nodes (namely, $`m+1`$). Figure 2(b) shows some examples of connected subnetworks with $`m=6`$. Obviously, increasing the number of nodes involved in the shortcut subnetwork reduces $`L`$, since it reduces the average path length to the nearest shortcut node. Among all connected configurations of shortcuts having $`m+1`$ nodes, the ones having a single center give the shortest value for $`L`$, since the average path length of the shortcut subnetwork is the smallest in that case.

These arguments indicate that given a fixed number of shortcuts, the networks with a connected subnetwork of shortcuts having nodes uniformly distributed have smaller $`L`$ than a typical random configuration, and among those the ones with a single center minimize $`L`$. In other words, the “smallest” small-world networks are characterized by these structures.

Now we will compute explicitly the average path length for a configuration with a single center in the case of small-world networks constructed from a one-dimensional lattice. Consider $`N`$ nodes arranged uniformly on a circle of unit circumference, where each node is connected to its two nearest-neighbor nodes. In addition, consider shortcuts connecting $`m`$ arbitrary pairs of nodes. To make the calculation simple, we take the continuum limit $`N\mathrm{}`$ with $`m`$ fixed, in which the network becomes a continuous graph composed of a circle corresponding to the lattice and chords representing the shortcuts. Let us define the distance $`d(P,Q)`$ between points $`P`$ and $`Q`$ on the continuous graph as the length of the shortest continuous path along the graph, *regarding the length of a chord as zero*. In other words, a shortcut is regarded as identifying two points on the circle, rather than merely connecting them. Then, the number of links in the shortest path between nodes $`P`$ and $`Q`$ in the original network, normalized by $`N`$, tends to $`d(P,Q)`$ as $`N\mathrm{}`$. This one-dimensional model, despite being one of the simplest models of small-world networks, captures basic features of many real networks. In Ref. Newman et al. (2000), a mean-field-type argument was used to derive an analytical expression for an average of $`L`$ over random configurations of shortcuts, which was later improved in Ref. BarBour and Reinert (2001). In the following, we derive an analytical expression for the configuration with a single center.

Consider the configuration of shortcuts with a center node connected to $`m`$ other points on the circle, as shown in Fig. 3.

The $`m+1`$ points including the center point are equally spaced with $`\xi 1/(m+1)`$, and they divide the circle into $`m+1`$ arcs of the same length. We will compute the average $`\overline{d(P,Q)}`$ taken over all pairs $`(P,Q)`$. Without loss of generality, we may consider $`P`$ as fixed. Let $`A_P`$ be the arc in which $`P`$ lies. Suppose first that $`QA_P`$ as in Fig. 3(a). Because the end points of $`A_P`$ are connected to each other by two shortcuts via the center, the distance in $`A_P`$ is equivalent to the distance on a circle of circumference $`\xi `$. Therefore, the average of $`d(P,Q)`$ over all pairs $`(P,Q)`$, such that $`QA_P`$, is equal to the average distance between two points on a circle of circumference $`\xi `$, which is $`\xi /4`$. Suppose now that $`QA_P`$ as in Fig. 3(b). Let us denote the distance from $`P`$ to its closest shortcut connection by $`\alpha `$, and the distance from $`Q`$ to its closest shortcut by $`\beta `$. Since the shortest path between $`P`$ and $`Q`$ must pass through two shortcuts of length zero, we have $`d(P,Q)=\alpha +\beta `$. Averaging this over all possible choices of $`\alpha `$ and $`\beta `$, which can take any value between 0 and $`\xi /2`$ independently, we obtain $`\xi /2`$. Noting that the probabilities that $`QA_P`$ and that $`QA_P`$ are $`1/(m+1)`$ and $`m/(m+1)`$, respectively, the normalized average path length $`l`$ can be calculated as

$$l=\overline{d(P,Q)}=\frac{1}{m+1}\left(\frac{\xi }{4}\right)+\frac{m}{m+1}\left(\frac{\xi }{2}\right)=\frac{2m+1}{4(m+1)^2}.$$

Let us now consider more general situation where each node in the network has connections to its neighboring nodes, up to $`k`$th nearest neighbors. Because of the connections to $`k`$th nearest neighbors, following the shortest path between nodes $`P`$ and $`Q`$ takes $`1/k`$ times less steps compared to the case discussed above. Hence, we must also scale $`l`$, the normalized average path length of the network, by a factor $`1/k`$ yielding

$$l=\frac{1}{k}\overline{d(P,Q)}=\frac{2m+1}{4k(m+1)^2}.$$
(1)
An important observation about Eq. (1) is that it can be written as $`l=f(m)/k`$, where $`f(m)`$ is a function that depends only on the number of shortcuts. The formula derived in Ref. Newman et al. (2000) for the average $`l_r`$ of normalized path length over random configuration of shortcuts also has the same form with different function for $`f`$, namely,

$$l_r=\frac{1}{2k\sqrt{m^2+2m}}\mathrm{tanh}^1\left(\frac{m}{\sqrt{m^2+2m}}\right).$$
(2)
Note also that since the shortcuts are considered to have length zero, the derivation above remains correct as long as the subnetwork of shortcuts is connected and has uniformly distributed nodes, suggesting that in the continuum limit these two conditions are sufficient to achieve the minimum of $`L`$.

Figure 4 compares the calculation summarized in Eq. (1) (continuous curve) with numerical computation of $`l`$ for a single center (circles) and of $`l_r`$ over 10 random configurations of shortcuts (squares).

This shows an excellent agreement of Eq. (1) with the simulation. In fact, the error in the Eq. (1) due to the approximation $`N\mathrm{}`$ is of order $`1/N`$, mainly because the normalized length of a shortcut is considered to be zero rather than $`1/N`$. The inset in Fig. 4 shows the ratio $`l_r/l`$ as a function of the number ($`m`$) of shortcuts. Here the ratio is computed from numerical simulations (circles) and from the theoretical results (1) and (2) (continuous curve). Since Eq. (1) is valid for $`mN`$ and Eq. (2) is valid for $`1mN`$, the curve in the inset is exact in the limit $`N\mathrm{}`$ with $`m1`$ fixed. Using the asymptotic form $`l_r(\mathrm{log}2m)/4m`$ of Eq. (2) for $`m1`$, one sees that $`l_r/l\mathrm{log}m`$, explaining the fact that the curve in the inset is almost a straight line for large $`m`$. Numerical results in the inset indicate that the effect of finite size and large shortcut density actually increases the ratio, making the benefit of optimizing the shortcut configuration to a single-center model even larger than the theoretical prediction.

Finally, we simulate optimization of the shortcut configuration for a one-dimensional array of nodes using the genetic-algorithm (GA) methodology Mitchell (1996). An initial population is described as being a collection of various shortcut configurations specified by $`m`$ pairs of integers representing the locations of nodes connected by shortcuts. The fitness of each configuration is defined to be $`L^1`$, where $`L`$ is the average path length. A new population of shortcut configurations is created from the old one in analogy with reproduction in population genetics: a configuration is viewed as being the genome of an individual in the population, and in creating a new population, we allow there to be one-point crossovers (i.e., interchanging subsets of shortcuts) and mutations (i.e., changes in the location of end points by Gaussian random numbers). This creation process is continued until the fitness of the best individual in the population is constant over 100 generations. This gives a candidate for the optimal solution. The program for the simulation was developed using a C++ library called GAlib gal .

Ten best solutions (here best means having shortest average path length) resulting from 254 independent runs with $`m=10`$, $`k=1`$, and $`N=1000`$, and the population size of 100 are shown in Fig. 5.

First, observe that in each case the subnetwork of shortcuts is connected. This was the case in every solution found using the genetic algorithm. Second, in each case there are centers from which many shortcuts emanate. Moreover, the nodes in the subnetwork are approximately equally spaced around the circle. These observations are consistent with the argument used above to establish our results. All solutions in Fig. 5 have the average path length within 2% of the average path length achieved by the single-center configuration (which is 44.577). In contrast, the corresponding value for random shortcuts ($`88`$) is almost double the single-center solution. Although the single-center solution was not found by the genetic algorithm due to the limited number (254) of simulation runs, the results show that configurations with several centers are almost as efficient as the single-center configuration, as long as the subnetwork of shortcuts is connected and its nodes are uniformly distributed. The single-center solution was found for smaller networks with $`N=100`$ and $`m=5`$, for which the computation is less demanding.

Any other values of $`k`$ should lead to similar results. The case of $`k=2`$ is shown in Fig. 6. In fact, due to the generality of the argument given earlier, we expect that the results can be extended to the case where the shortcuts are added to a lattice of higher dimension, or to a regular network of another type.

The result of these simulations using the GA methodology shows that design elements for efficient networks are (1) connectedness of the shortcut subnetwork, (2) uniform distribution of nodes in the subnetwork, and (3) existence of centers.

We expect to see many examples of real networks with such structures. Our computations on the neural network of *C. elegans* (which has 285 nodes, 2347 links, and 112 shortcuts) show that the structures are indeed present: (i) the shortcut subnetwork has much fewer (= 15) connected components than the average ($``$ 47) for randomly chosen shortcuts, and the size of its giant component (= 75) is significantly larger than the average ($``$ 12) over random shortcuts; (ii) most ($``$ 88%) of the nodes are within one step of a shortcut; (iii) there are a few nodes having many shortcuts (11 shortcuts in the main center). In general, a network with such structures is robust against random failures, although it is sensitive to deliberate attacks to the centers. This property, which is shared by scale-free networks Barabási and Albert (1999), is shown to characterize many real networks such as the Internet and the WWW Albert et al. (2000). However, some biological networks may be robust even against attacks on the centers since loss of a center can result in shortcuts reconnecting to nearby nodes followed by the optimization process that quickly recovers the smallest configuration.

We have shown that among the small-world networks having a fixed number of shortcuts, the average path length is smallest when there exists a single center through which all of the shortcuts are connected and shortcut nodes are uniformly distributed in the network. We have also shown that the average path length is almost as small when the shortcuts are connected and have a few centers, which was supported by the result of the GA simulations. Our results have important consequences in situations where the efficiency of information flow over a large network is required. The fact that the architecture of connected shortcuts with centers arises through genetic algorithms suggests the possibility that such a structure could emerge in networks in natural organisms (e.g., the neural network of *C. elegans*), although the fitness used in GA here is not necessarily related to that of natural selection in biology. In particular, it provides a potential mechanism for the appearance of highly connected nodes while keeping high clustering in networks that are evolving but not necessarily growing, such as neural and metabolic networks.

T.N. was supported by DARPA/ONR through Grant N00014-01-1-0943. A.E.M. thanks Fapesp for financial support. Y.C.L. was supported by AFOSR under Grant No. F49620-01-1-0317. F.C.H. was supported in part by NSF through Grant DMS-0109001. We are grateful for helpful suggestions from A. P. S. de Moura.