text:
# Extensional rupture of model non-Newtonian fluid filaments

## I Introduction

The extensional dynamics of polymeric, non-Newtonian liquids controls processes ranging from commonplace activities such as chewing gum to large-scale commercial applications in extrusion processing of materials . Understanding the dynamics of the liquid in these situations  requires information beyond the usual shear viscosity characterization of flowing fluids, and has led to extensive experimental and theoretical work on “liquid bridges” . Here, one pulls on the ends of a cylinder of the liquid and observes shapes, forces, velocities and (to some degree) stresses, and attempts to relate the behavior to an appropriate rheological model. The ultimate fate of an extending liquid bridge is a rupture at some point, and one component of the efforts in this subject has attempted to understand the fracture mechanisms involved. Aside from practical relevance to materials failure, the thinning dynamics of a cylindrical neck of liquid has become an important subject in fluid interfacial dynamics . In almost all cases, these efforts have concerned macroscopic bodies of liquid, in the sense that the experiments do not generally have atomic resolution, and the calculations are based on continuum descriptions. The purpose of this paper is to consider the interfacial rupture of a liquid bridge configuration at an atomic scale. In addition to the intellectual question of exploring the late stages of rupture that are below the level of resolution considered previously, a further motivation comes from recent interest in micro- and nano-fluidics , where precisely these length and time scales can be reached in practice.

The focus of this paper is somewhat outside the issues usually considered in rheology. Previous studies of the rupture of liquid filaments  have concentrated on determining instability criteria based on continuum modeling. The focus has been on constant strain-rate deformations, and on identifying characteristic regions of the strain rate – maximum strain phase space, corresponding to viscous, elastic, rubbery and glassy behavior. It has proven useful to distinguish a recoverable component of the strain, representing the residual deformation of the material when the stretching forces are released and it is allowed to retract (i.e., omitting viscous flow deformations). It is inferred that true rupture occurs only when this recoverable strain exceeds 1/2. In the systems studied here, the chain molecules can distort in shape but not divide or combine, but one possible mechanism for a nonrecoverable strain could involve permanent disentanglement of the chains. A second focus of previous work has been the relation between time to rupture and applied stress; in our simulations the ends of the filament are moved at prescribed velocities, so that applied stress is not fixed and instead fluctuates during the course of the stretch. To discuss other configurations, empirical correlations or the results of non-Newtonian continuum modeling are used. The simulations performed here quickly lead to regions of filament which are only a few molecules in diameter, and not obviously described in these ways.

We have previously considered the rupture of fluid interfaces in two situations, using molecular dynamics (MD) simulations. In a study of the dynamics of a monatomic drop placed in an second immiscible fluid undergoing shear , we found fair quantitative agreement with experiments and continuum calculations in general, and in particular that the drop ruptures at a sufficiently large capillary number. The microscopic mechanism involved is a gradual flow of the drop’s molecules to one side or another induced by the external flow, coupled to surface tension (originating in the molecular interactions chosen to provide immiscibility) acting to reduce the area of the interface separating the drop and solvent liquids. In a previous paper on non-Newtonian fluids  we used MD simulations to study liquid bridge deformation, focusing on the shape evolution and internal velocity and stress fields. In that paper relatively small samples of several different model molecular fluids were considered – a Newtonian liquid based on short chain molecules, a solution of a few longer (20-40 monomer) molecules in a short-chain solvent, and melts of the longer chains. Although the simulations were continued to the point of rupture, these systems were too small to display any interesting behavior in that regard which could be quantitatively analyzed. The Newtonian liquid and melt systems ruptured due to fracture after only modest length increases, by less than a factor of two. The model solution system attained rather longer lengths, by virtue of the longer molecules becoming extended and providing a kind of backbone, and ruptured when the overall length exceeded the amount of molecular backbone available, essentially following the monatomic mechanism. Our goal here is to use a larger system based on a melt of longer molecules, and study the deformations of highly elongated liquid filaments.

## II Simulation Procedure

The simulational method and the molecular system closely follows our earlier papers on nanodrop dynamics , liquid bridge extensional flows , and single chains in a flow . The MD calculations themselves are of standard form , and consider a fluid made of monomers interacting with two-body forces of two types, a Lennard-Jones potential to provide an impenetrable core and chemical attraction at larger distances, along with a finitely-extensible nonlinear elastic force to link the monomers into model polymer chains. The explicit forms of the interaction potentials are

$$V_{\mathrm{LJ}}(r)=4ϵ\left[\left(\frac{r}{\sigma }\right)^{12}\left(\frac{r}{\sigma }\right)^6\right]V_{\mathrm{FENE}}(r)=\frac{1}{2}kr_0^2\mathrm{ln}\left(1\frac{r^2}{r_0^2}\right)$$
$`(1)`$
The Lennard-Jones interaction acts between any two monomers separated by a distance less than $`r_c=2.5\sigma `$, and a linear term is added to $`V_{\mathrm{LJ}}`$ so that the force vanishes at this cutoff. The FENE potential  acts only between adjoining monomers in any one molecule, and has the effect of limiting the bond length to $`r_0`$. The parameters in the FENE potential are taken to be $`r_0=1.5\sigma `$ and $`k=30ϵ/\sigma ^2`$, following Kremer and Grest , and are chosen to yield an entangled melt of non-crossing chains for the 100-monomer chain length used here. The characteristic time for monomer dynamics is defined by the LJ parameters, $`\tau =\sqrt{ϵ/m\sigma ^2}`$, where $`m`$ is the monomer mass. In the remainder of this paper, and in particular in all figures, when dimensions are not given we are implicitly using “MD units” where $`\sigma `$ is the unit of length, $`\tau `$ is the unit of time, and $`ϵ`$ the unit of energy.

A molecular liquid based on the FENE and LJ potentials (1) has been widely studied as an MD a model of polymer melts . In most previous work, the 6-12 potential is cut off at a smaller value than used here, $`r_c=2^{1/6}\sigma `$, which has the computational advantage of reducing the number of particles which are in interaction range, and speeding the calculation. The properties of this system are given in detail in . This value of $`r_c`$ corresponds to a purely repulsive potential which, while adequate for a confined material, has the disadvantage here that a fluid will expand to fill space uniformly, with no interface. Since an interface is crucial here, attraction is needed and we retain the larger value of $`r_c`$. However, because the interaction has been altered, the liquid studied here need not have the same transport coefficients. In separate simulations we find a static shear viscosity $`\mu (0)10^3m/\sigma \tau `$, defined as the limit of the dynamic shear viscosity $`\mu (\dot{\gamma })`$, the ratio of shear stress to strain rate in Couette flow, when the strain rate tends to zero. Typically polymer melts based on LJ and FENE potentials show shear thinning – a nearly constant viscosity up to a certain transition shear rate $`\dot{\gamma }^{}`$, and a systematic (roughly power-law) decrease at larger values. We use the inverse of the transition shear rate as a characteristic fluid relaxation time, $`\lambda 1/\dot{\gamma }^{}10^5\tau `$, in this case. This number and the quoted viscosity above are very crude estimates: at the very low velocities needed to reach the plateau for this chain length, the signal-to-noise ratio in a stress measurement is poor, and even in very long runs substantial fluctuations are found. The rheology of chain liquids and the definition of characteristic time in this context is discussed further in ; the latter paper considers MD simulations of the liquid bridge extensional dynamics of model melts (and other fluids) and also shows that they exhibit strain-hardening in extensional flow, in the sense that an appropriate extensional viscosity increases with strain rate. Lastly, the surface tension is measured to be $`\alpha =0.37\pm 0.01`$ using an independent equilibrium simulation of a slab of liquid with planar interfaces .

The MD calculations are performed in parallel on either T-3E, Origin-3000 or IBM SP platforms, using a spatial domain decomposition algorithm, a layered linked-cell list on each processor, and MPI routines. The numerical procedure is to integrate Newton’s equations of motion using the velocity Verlet algorithm  with time step 0.005$`\tau `$. Occasionally the computation was checked by re-running over a certain time interval with a smaller step length. A Nosé-Hoover thermostat  is used to maintain the liquid at a constant temperature $`T=1.5ϵ/k_B`$. The initial configuration, a roughly-circular cylinder of length $`L_0=80.5\sigma `$ and radius $`R_0=26.3\sigma `$, containing 1400 chain molecules 100 monomers in length, is in a disordered liquid state resulting from a somewhat convoluted preparation procedure. In our earlier simulations of length-100 chain molecule melts , we began with the monomers on regular lattice sites and randomized the configuration by cookingit at a high temperature and low density in a box with repulsive walls. This system was equilibrated to the point where its time-averaged monomer spatial density was uniform, the probability distributions of molecular sizes (end-to-end length and radius of gyration) reached a steady state, and the probability distributions of molecular orientations were rotationally invariant. A spherical drop was formed by cooling the liquid gradually to temperature 1.5, while simultaneously applying a weak central force until a stable sphere formed, and then gradually switching off the force. Next, half of the sphere was placed on a partially-wetting atomistic solid substrate and re-equilibrated, and finally the sphere was allowed to coalesce with a slightly randomized copy of itself. Because both hemispheres were placed on attracting substrates, the result was a liquid cylinder held between two solid surfaces, rather than a single larger sphere. Removing the solid substrate after coalescence gives the initial cylinder considered here. One might be concerned that when the cylinder is pulled it would tend to break at the mid-plane into its two original drops due to incomplete healing, but in fact we observe rupture at various other locations.

To stretch the filament, we simply pull on its ends. More precisely, each monomer in a disk of thickness 5$`\sigma `$ at each end is attached to a “tether site” by a stiff spring, and the tethers are moved outwards in a prescribed manner by displacing their position by the appropriate axial amount $`\pm v_0\mathrm{\Delta }t\widehat{z}`$ at each time step. In simulations of constant-velocity extension, $`v_0`$ is just the velocities at each end of the filament, while for constant-strain-rate extension where the length of the filament grows as $`L(t)=L_0e^{\dot{ϵ}_0t}`$, the fixed control parameter is the strain rate $`\dot{ϵ}_0`$, and the velocity is $`v_0(t)=dL(t)/dt`$. This procedure is similar to a common experimental one  where the liquid filament rests between solid end plates which are moved apart, but has the conveniences that the end plates need not be modeled and that flow along the end plates does not arise. The focus of this paper is the behavior of the liquid in the “middle” of the filament, so nothing is lost in this way. We also attempted an alternative procedure suggested by experiments on single DNA molecules  in which a constant axial force was applied to each monomer in the end regions, but the result was that the end regions contracted radially and ruptured occurred nearby, leaving the bulk of the filament unextended.

The operating conditions of the various simulations performed here are listed in Table I, in terms of the pulling velocity or strain rate. Note that in either case the initial velocity and strain rate are related by $`v_0=L_0\dot{ϵ}_0`$, where $`L_0`$ is the filament’s initial length. In terms of the usual dimensionless groups characterizing a flow, at the start of a run the Reynolds number $`Re=\rho v_0R_0/\mu (0)`$ is always much less than unity, while the Deborah number $`De=\lambda \dot{ϵ}_0`$ and the capillary number $`Ca=\mu (0)v_0/\alpha `$ are O($`10^2`$) or larger. Although the Reynolds, Deborah and capillary numbers describing this flow are all moderate by laboratory standards, it should be kept in mind that the physical velocities are relatively large, meters to hundreds of meters per second (depending on one’s interpretation of the Lennard-Jones monomers), and not at all small for the microscopically thin bodies of liquid considered here.

## III Constant Velocity Stretching

We first consider filament stretching at a relatively low constant velocity, displacing the two ends at $`\pm 0.05\sigma /\tau `$. Successive snapshots of the filament appear in Fig. 1, in a low-resolution side view intended to indicate the overall shape. At early times there is a mostly-uniform deformation, aside from some necking near the two ends which results from the tethering method used here – the monomers there are essentially pinned to their initial relative positions and when these are moved apart the melt attempts to fill in. A “dynamical” necking down appears near the center of the filament around $`1500\tau `$, followed by relatively uniform stretching of the central neck region through $`3100\tau `$, while the two end segments of the filament translate with the tethers with only minor changes of shape. Eventually the filament breaks when the neck region is so stretched as to run out of molecules. As seen in the close-up views in Fig. 2, the neck eventually thins down to a few molecular chains lined up along the stretch axis, and the breakage occurs when these chains are pulled in opposite direction and slide past each other. Note that the molecular model used here does not incorporate chain breaking, so the molecular “depletion” seen here is the only available rupture mechanism.

Calculations and laboratory experiments  on the thinning of a continuum Newtonian liquid cylinder typically find a power-law variation of minimum radius with time. Here, however, while we do observe a systematic decrease as seen in Fig. 3, the detailed time variation is not well fit by a power of time. Aside from the fact that the later stages of rupture are associated with a neck made of just a small number of highly extended molecules, it is not unreasonable to see quantitatively different behavior. Furthermore, at this stage the liquid is presumably highly viscoelastic so that Newtonian fluid predictions may not be applicable. Similar behavior is observed in all of the other cases presented below: although there is a systematic decrease in the filaments’ minimum radius with time, it does not fit a power-law variation. Figures 8 and 9 show a sequence of snapshots at higher (fixed) tether speeds, $`\pm 0.15`$ and $`\pm 0.5\sigma /\tau `$ respectively. In these various simulations, the initial liquid cylinder is always the same, and within each sequence the same length scale is used. First note that the rupture occurs at different points along the filament, indicating that the process is random and an indirect function of the particular sequence of fluctuating positions and velocities in each case, and furthermore that there is no special defect associated with the “sample preparation”. At these higher speeds, the stagnant end regions at the ends do not develop, and more molecules are available to the central portion of the filament, allowing it to attain longer strains before rupture. Here the ends are steadily drawn out, leaving a central region which does expand but at a distinctly slower speed, and begins to evolve into isolated drops. In the central region, relatively rapid thinning occurs at several intermediate points, and the configuration prior to rupture resembles the results of the quasi-static Rayleigh instability , where a cylinder of liquid without imposed flow breaks up into spherical drops. The final stages of rupture, again, amount to a very narrow neck only a few molecules in width attenuating down to nothing as its elongated chains are pulled past each other. The shape evolution is in a sense opposite to that of the low-velocity case in Fig. 1, where the end regions began to translate passively before becoming highly stretched. Finally, as indicated in Table I, note that the maximum length attained before breakup increases with pulling speed while the time to rupture decreases, presumably because the higher stress allows more molecules to be stretched.

## IV Constant Strain Rate Stretching

The more common experimental protocol in liquid bridge systems is stretching at constant strain rate, $`\dot{ϵ}_0=\dot{L}(t)/L(t)`$, leading to an exponential increase in length and velocity. Ostensibly the motivation is to provide a constant “strain environment” for the liquid, to facilitate study of its extensional rheology, although in practice  the strain can often be very unevenly distributed. In Fig. 10 we show some of the time-sequence of shapes observed at a strain rate $`\dot{ϵ}_0=0.001\tau ^1`$. The evolution combines features of the various constant-velocity deformations presented above: the filament narrows in the center and eventually ruptures there due to attenuation of molecules, but at the same time the end regions are not translating stagnantly but rather stretch at a slower rate than the bulk, and at the time of rupture are themselves in the process of thinning and presumably would eventually break up due to the Rayleigh instability.

The distinction between constant-velocity and constant-strain deformations is that in the former case the strain rate (and stress) decrease with time, so that after an initial stage where the system is pulled out of its equilibrium configuration, the molecules can relax back; see Fig. 4. In the present case, the stress is ostensibly constant, and the molecules are increasingly elongated and oriented with no configurational relaxation until nearly the time of rupture, as shown in Fig. 11. The behavior of the velocity profiles is also consistent with these remarks; as seen in Fig. 12, the velocity is closer to a linear variation along the axis for a considerably longer time interval.

A somewhat related MD study  of the elongation and relaxation of a model polymer melt considered chains of 30 monomers using the same interactions as in this paper (but with a short $`r_c=2^{1/6}\sigma `$ LJ cutoff). Instead of simulating the explicit stretching of filaments, a nonequilibrium MD method was used to apply a constant-strain uniaxial elongation to a homogeneous system in a box. Among other things, the behavior of the chain size (radius of gyration) and orientation was measured, and the qualitative behavior is similar to that seen here. The elongation and degree of orientation increased monotonically during the flow, at a rate which increased monotonically with the strain rate, and then relaxed monotonically when the flow was turned off.

We have considered the question of unrecoverable strain  in this simulation, by releasing the tethered monomers and allowing the liquid to take up a new equilibrium shape without external stress. If the filament is released well before rupture, it seems to withdraw into a sphere, whereas if released just before rupture the thinning neck continues to thin, rupture occurs anyway, and afterwards the liquid appears to be withdrawing into two spheres. Here the phrase “appears to be” means that as a function of time the liquid shape is tending towards one or two spheres at the original liquid density, but the process is so slow that we have not followed it all the way to equilibrium. The asymptotic strain does, however, go well below 1/2, and we do not see the values of unrecoverable strain suggested in  Similar behavior is found in the other cases as well. Given that the initial simulated filament was prepared by an equilbration procedure involving relaxation of a dilute high-temperature melt, it is difficult to see how the final state here would differ in any statistically significant way, other than perhaps the time needed to reach it.

At higher strain rates a different behavior appears, a cohesive failure of the melt in the region where the molecules are tethered to pulling sites. In Fig. 13, at strain rate $`\dot{ϵ}_0=0.004`$, there is a rapid thinning at the ends of the filament leaving a central stagnant region showing the onset of a Rayleigh instability, while the thin necks attaching this region to the ends attenuate until breakup. At the same time, some of the tethered molecules are pulled out of the filament, and the fluid bodies at the ends have not had time to relax to an equilibrium shape. At the still higher strain rate $`\dot{ϵ}_0=0.01`$ shown in Fig. 14, the tethered molecules are simply pulled free of rest of the the filament, and rapid rupture occurs. The rupture times and lengths for the various cases are recorded in Table I.

## V Force Measurements

As an alternative to the stress fields, we have studied the forces from a more direct viewpoint on the monomer scale. The quantity of perhaps most direct experimental relevance is the pulling force exerted on the tethers at the ends. While this quantity is easy to evaluate in the simulations, it is unfortunately quite noisy. Aside from the problem noted above of large atomic-level force fluctuations in general, here the procedure used to generate the motion exacerbates the problem. The atoms at the ends of the filament are attached to moving tether sites by stiff springs, whose force is very sensitive to small atomic displacements.

In the constant-velocity simulations the result is simple – the pulling force is roughly constant in time, fluctuating about a value which is approximately linearly proportional to the pulling velocity. At the three velocities studied, 0.05, 0.15 and 0.5$`\sigma /\tau `$, the force per pulled monomer is measured to be 2.72$`\pm `$0.16, 7.31$`\pm `$0.24 and 24.9$`\pm `$0.25$`ϵ/\sigma `$, respectively. (Typical values of the instantaneous force on an monomer at this density and temperature are O(1-10) in this unit, while its time averaged values over a 100$`\tau `$ interval are an order of magnitude smaller.) One can rationalize the time-constancy of the force by arguing that in this case the molecules are moving past each other with a roughly fixed value of relative velocity, and since the material is a viscous liquid in a low Reynolds number flow, there is a constant drag force on each molecule.

In contrast, in the constant strain rate cases there is a non-trivial time variation. For example, as shown in Fig. 15 for the lowest rate $`\dot{ϵ}_0=0.001\tau ^1`$, the tether force jumps to a finite value once the force is applied, then increases with time until about 1750$`\tau `$. After a sharp drop, the force again increases until rupture. Referring to Fig. 10, after this transition time the right and left sides of the filament are approximately translating without a significant shape change, while the central region thins down to nothing. The early increase can then be associated with strain hardening  and the shape deformation of the entire liquid filament; the sharp drop with relaxation of stress after some internal readjustment and formation of the central neck; the second increase with the strain hardening of the central thinner region; and the final drop-off with rupture. At the next highest strain value, 0.004$`\tau ^1`$, the numerical values of the tether force are comparable, and the tether force shows a similar variation, with a corresponding drop at about 600$`\tau `$, the time at which (see Fig. 13) the central region of the filament halts and ceases to change its shape while the regions to either side continue to elongate. At still highest strain rates there is no obvious pattern to the tether force evolution; the fragmentation of the regions near the tethers makes it difficult to identify simple mechanisms.

We can also compute the forces in the interior of the filament, and in Fig. 16 we show an example of the force in the axial direction in the $`ϵ_0=0.001`$ constant strain rate simulation at time 2150$`\tau `$, averaged over a 100$`\tau `$ time interval and a 6$`\sigma `$ spatial interval along the axis. The force on the end regions containing the tethered monomer shows large spatial and temporal fluctuations for the reasons given at the beginning of this Section, and these two regions have been suppressed in this plot. Aside from the fluctuations, the trend is a roughly-linear ramp with positive (negative) values corresponding to the halves of the filament moving to the right (left). As a function of time, the fluctuations shift around apparently randomly, and the height of the ramp is roughly constant in time. At higher strain rates, the analogous plot is again a fluctuating quasi-linear ramp in its central region, but shows rapid spatial fluctuations near the ends, corresponding to molecules being pulled out there. In contrast, in the constant-velocity simulations there is always a linear ramp as a function of $`z`$, and furthermore, at different pulling speeds the ramp height is roughly linear in velocity.

## VI Conclusions

We have discussed the results of molecular simulations of a liquid filament of model polymer melt stretched to the point of rupture. These model systems are composed of non-breakable, freely-jointed, entangled chains, and the rupture process typically involves ductile failure of some sort. Here we observed these filaments stretch into a very non-uniform shape, with one or more thin necked regions, and then fail when one of the necks simply attenuates down to nothing. More precisely, during the stretching process the thinning neck regions are populated with well-extended chain molecules lying roughly parallel to the extension axis, which are drawn past each other to one side or the other by their respective interacting neighbors. Eventually the filament is held together by only a few such adjacent elongated chains, and splits once these slide past each other. In continuum terms, the neck thins under a combination of mass transfer due to outflow and (Rayleigh) capillary instability. In analogous earlier simulations in systems of melts or solutions based on shorter length molecules (up to 20 monomers), presumably unentangled, a more familiar form of ductile fracture was observed in which a neck several molecules in thickness separated into two pieces . In that case, a model solution was able to attain the longest lengths prior to breakup, and there one could see that the fully-extended longer molecules acted as a backbone which stiffened the filament, and breakage occurred when the filament was extended to the point that these backbone molecules were no longer in contact with each other.

In term of the classification of failure zones by Malkus and Petrie , the constant velocity and lowest-rate constant strain runs are in the “flow zone”, where viscous deformation can procede indefinitely provided there is enough liquid. The two runs at higher constant strain appear to have “glassy states”, at least near the tethers where rupture occurs.

The filament shapes observed here seem unusually irregular compared to those studied elsewhere, and in no case do we observe a long neck of nearly constant diameter. Note however, that these simulations begin with a liquid sample whose diameter is already in the nanometer range, and so these results should at best be compared to the last, usually poorly-quantified stages of typical laboratory studies. The surface irregularities and fluctuations are, however, quite characteristic of molecular simulations of moving microscopic liquid bodies with free surfaces . In consequence of these fluctuations, the measured continuum fields show sufficiently large statistical fluctuations as to be non-trivial to understand even qualitatively. Even microscopic measurements which do yield a smooth result, such as the variation of minimum filament radius with time, are not in agreement with continuum calculations. The molecular level information is likewise fluctuating, but does seems to be readily interpretable at this scale.

A further issue in comparing these results to laboratory measurements is the pulling speed involved. Although the standard dimensionless groups characterizing these flows, the Reynolds, Deborah and capillary numbers, are similar to those considered elsewhere, there is quite a difference in velocities. In typical liquid flows the drift or peculiar velocity is a small fraction of the thermal velocity of the monomers, but here, particularly in the late stages of the constant strain simulations, the tether velocity approaches the thermal velocity, O(1$`\sigma /\tau `$). The severe fragmentation of the end regions of the filament is presumably one aspect of this high speed, but a second is that surface tension effects which would ordinarily smooth a free surface may not have time to act.

## Acknowledgments

We thank the NASA Office of Physical and Biological Sciences for financial support, and the NASA Advanced Supercomputing Division at the Ames Research Center, and NPACI at the San Diego Supercomputer Center for providing computational resources.

|  | Constant velocity | | | Constant strain rate | | |
| --- | --- | --- | --- | --- | --- | --- |
| $`v_0`$ or $`\dot{ϵ}_0`$ | 0.05 | 0.15 | 0.5 | 0.001 | 0.004 | 0.01 |
| length | 432 | 694 | 1080 | 701 | 1594 | 959 |
| time | 3520 | 2045 | 1000 | 2250 | 763 | 263 |

Table I. Final length and time to rupture (MD units).


format: text

text:
# Perturbative Inaccessibility of Conformal Fixed Points in Nonsupersymmetric Quiver Theories

## Abstract

The possibility that non-supersymmetric quiver theories may have a renormalization-group fixed point at which there is conformal invariance requires non-perturbative information.

The AdS/CFT correspondence of Maladacena provides a powerful tool to study non-gravitational gauge field theories in flat e.g. four-dimensional spacetime. In particular, it has been suggested  that nonsupersymmetric quiver gauge theories obtained from compactification of the Type IIB superstring on $`AdS_5\times S^5/\mathrm{\Gamma }`$ may, for finite number $`N`$ of $`D3`$branes, possess a conformal fixed point at the TeV scale and that one such theory may be the correct path to go beyond the standard model of particle phenomenology.

Three years ago, in a paper an argument based on perturbation theory was used to criticize this whole approach to string phenomenology. Our purpose here is not just to comment on that paper but on what we perceive to be a widespread belief that the way to approach the issue is through perturbative analysis.

Here we shall show that the S-duality of the underlying IIB superstring constrains the relevant RG beta function $`\beta (\alpha )`$ to suggest a conformal fixed point but that to any order of perturbation theory such a fixed point must remain inaccessible. We shall show this by an illustration of the possible form of $`\beta (\alpha )`$ which may be analytic around $`\alpha =0`$ and $`\alpha =\mathrm{}`$ as a function of $`\alpha `$ with arbitrary coefficients in the Taylor expansion for small $`\alpha `$.

Let us assume that the quiver theory has one independent coupling constant which is asymptotically free for small $`\alpha `$ so that the perturbative expansion is (here $`t=\mathrm{ln}\mu `$)

$$\beta =\frac{d\alpha }{dt}=b_0\alpha ^2+b_1\alpha ^3+b_2\alpha ^4+\mathrm{}$$
(1)
with $`b_0<0`$. As already mentioned in  the underlying S-duality under $`\alpha 1/\alpha `$ implies by

$$\beta (1/\alpha )=\frac{1}{\alpha ^2}\beta (\alpha )$$
(2)
that $`\beta (1)=0`$.

Let us illustrate the solution of Eq.(2) in the form

$$\beta (\alpha )=\alpha (\alpha \frac{1}{\alpha })F(\alpha )$$
(3)
where $`F(\alpha )=+F(1/\alpha )`$. In particular consider

$$F(\alpha )=\mathrm{\Sigma }_{n=0}^{n=\mathrm{}}C_n\alpha ^{n+2}(1+\alpha ^{2n+4})^1$$
(4)
Then

$$\beta (\alpha )=(1\alpha ^2)[C_0\alpha ^2(1+\alpha ^4)^1+C_1\alpha ^3(1+\alpha ^6)^1+\mathrm{}]$$
(5)
so that

$$\beta (\alpha )=C_0\alpha ^2+C_1\alpha ^3+(C_2C_0)\alpha ^4+\mathrm{}.$$
(6)
can reproduce any expansion of the perturbative type Eq.(1) by appropriate choice of the coefficients $`C_n`$. This form for $`\beta (\alpha )`$ can be analytic for both $`|\alpha |<1`$ and $`|\alpha |>1`$ but must be singular at least somewhere on the unit circle $`|\alpha |=1`$. In particular, $`\beta (\alpha )`$ can be analytic and vanishing at the points $`\alpha =0\mathrm{and}\mathrm{}`$. The specific illustration chosen has a simple zero on the real axis at $`\beta (1)=0`$. It is clear that analysis based on perturbative calculation of the coefficients $`b_n`$ in Eq.(1) could neither prove nor disprove the existence of this zero in $`\beta `$.

The point is that while calculation of the perturbative expansion for $`\beta (\alpha )`$ in Eq.(1) may well lead to supporting evidence (but not rigorous proof) of a fixed point for weak coupling, as in e.g. the well-known example of QCD with 16 flavors, such a calculation can never by itself confirm the presence or absence of a conformal fixed point.

The work of PHF was supported in part by the Office of High Energy, US Department of Energy under Grant No. DE-FG02-97ER41036. We thank Luis Alvarez-Gaume for a discussion.


text:
# Co-evolutionary games on networks

## I Introduction

Much research has been devoted to the statistical physics of complex systems with game theoretic interactions recently. One motivation are economic systems composed of a large number of agents with simple local interactions giving rise to complex global structures and dynamics Arthur (1999). In particular, the problem of stability and uniqueness of equilibria in economic systems has been readdressed in the context of the aggregate behavior of individual agents Kirman (1989). Game theory Neumann and Morgenstern (1953) and the theory of evolutionary games Maynard Smith (1982) provide a sufficient framework for modeling individual interactions whereas spatial structures has to be taken into account to tackle co-evolutionary dynamics of real-world systems Blume (1995); Nowak et al. (1994a); Thompson and Cunningham (2002).

Here, we will consider random networks of agents which face a social dilemma or, in physical terms, a frustrated interaction. Imagine a situation where each player can take two actions, say cooperating or defecting. The optimal global outcome would be achieved by all players cooperating. But an individual player can gain much more when exploiting the cooperators by defecting. Such a situation is called a social dilemma or a frustrated interaction. The central question is how social order is possible and how cooperative behavior can emerge. Examples for such spatially extended dilemmas are biological networks, where connected plants may or may not decide to share resources Oborny et al. (2000), the analysis of internet congestion Huberman and Lukose (1997), models for economic communication Miller et al. (2002), and, of course, many sociological problems from conflict research to public transportation Axelrod (1984); Schulz et al. (1994).

A simple model system is given by the iterated Prisoner’s dilemma (IPD) Axelrod and Hamilton (1981); Axelrod (1984) with co-evolutionary dynamics. The Prisoner’s dilemma game is probably the most prominent example of a basic model for the emergence of cooperative behavior in social, economic, and biological systems. It provides a frustrated two-particle interaction and has been extensively studied by physicists, economists, biologists, and mathematicians.

A spatially extended Prisoner’s dilemma was first proposed by Axelrod who concluded that territoriality strongly influences the evolution of cooperation Axelrod (1984). Extensive work on the spatial Prisoner’s dilemma started in 1992 when Nowak and May explored a cellular automaton based on this game on regular lattices. They and others found complex spatiotemporal dynamics and emergence of cooperation for strategy spaces confined to the strategies defecting and cooperating Nowak and May (1992, 1993); Huberman and Glance (1993); Nowak et al. (1994b, c); Herz (1994). For the Prisoner’s dilemma on lattices and strategy spaces confined to only cooperating and defecting (and Tit-For-Tat in Szabó et al. (2000)), methods from theoretical physics, as Monte-Carlo simulations, percolation theory, the theory of (nonequilibrium) phase transitions, and the concept of self-organized criticality, were used to understand why cooperators or defectors dominate or coexist in the system Szabó and Tőke (1998); Chiappin and de Oliveira (1999); Szabó et al. (2000); Vainstein and Arenzon (2001); Tomochi and Kono (2002); Lim et al. (2002). Lindgren and Nordahl introduced players which act erroneously sometimes, allowing a complex evolution of strategies in an unbounded strategy space Lindgren and Nordahl (1994). Others found that the limitation of a player’s memory to the last encounter, which translates to a bounded strategy space, does not provide a significant drawback for the players Hauert and Schuster (1997, 1998); Hauert (1999). Evolutionary games on networks, again with only two strategies, have been studied to ask questions how spatial organization influences the transition from defecting to cooperating Abramson and Kuperman (2001) and how the players themselves may influence the network topology Zimmermann et al. (2001). In the following, we will study the Prisoner’s dilemma on a network with a larger but bounded strategy space and co-evolutionary dynamics that lead to Nash equilibria as stationary states. It will be shown both numerically and theoretically that payoff matrix, strategy space, and topology are crucial to answer the question which equilibria will occur and how stable they are. In particular, critical avalanches of mutations are observed for such games and will be explained in detail.

This paper is organized as follows. In Section II, the spatially extended iterated Prisoner’s dilemma is described as well as its co-evolutionary dynamics. This is followed in Section III by numerical investigations of avalanche dynamics showing three distinct regimes due to changes in payoff matrix and topology. The observed Nash equilibria are described and explained in Section IV which enables us to understand the critical value of the control parameter of the payoff matrix. A confined branching process is introduced in Section V clarifying the relaxation mechanism and the emergence of scale-free behavior. Conclusions are drawn in Section VI.

## II A co-evolutionary spatially extended IPD

We start with a network with $`N`$ players as nodes where each player plays an iterated Prisoner’s dilemma game against each of its neighbors. The Prisoner’s dilemma is a two-person game with two possible actions in each encounter. The payoff matrix of the first player for the two strategies cooperating and defecting (denoted by $`\widehat{s}_1`$ and $`\widehat{s}_2`$) is given by

$$A=\left(\begin{array}{cc}3& 0\\ 5& 1\end{array}\right)=(a_{ij})_{i,j\{1,2\}},$$
(1)
with the entries $`a_{ij}=\widehat{\pi }_1(\widehat{s}_i,\widehat{s}_j)`$. $`\widehat{\pi }_1(\widehat{s}_i,\widehat{s}_j)`$ is the payoff of player 1 if player 1 plays strategy $`\widehat{s}_i`$ and player 2 plays $`\widehat{s}_j`$. The game is symmetric, i.e. $`\widehat{\pi }_1(\widehat{s}_i,\widehat{s}_j)=\widehat{\pi }_2(\widehat{s}_j,\widehat{s}_i)`$. Therefore, the corresponding payoff matrix of the second player is the transpose of the first player’s matrix. The Prisoner’s dilemma game in general is defined by the relations

$$a_{12}<a_{22}<a_{11}<a_{21}\text{ and }a_{12}+a_{21}<2a_{11}.$$
(2)
Hence, in one encounter of the Prisoner’s dilemma, defecting is the strategy that yields the best payoff regardless of the opponent’s strategy. This is no longer the case in the iterated game where mutual cooperation is more favorable than both players defecting or switching between defecting and cooperating. That is the reason why this system is a frustrated system. In each encounter, defecting would maximize a player’s payoff. But in the long run, when players will anticipate each other’s action, cooperation will in general do much better.

### II.1 IPD with memory on a network

Let us further specify the strategy space and payoff function of the spatial game. A strategy is viewed as a mapping of an agent’s “knowledge” to an “action”. “Knowledge” of an agent is given by the previous moves the agent can take into account to decide which action it will take next. We define the memory length $`m`$ of a player as the number of these previous moves and confine the memory of the agents to $`m=1`$, i.e. an agent remembers only its opponent player’s last action. If one player encounters another player it has to decide on its first move without any information about the opponent. Accordingly, the opening move is part of the strategy, too, which can be represented as a lookup table or a binary string (Tab. 1). The finite number of moves of one encounter is not known by any agent. In the course of the game, one player has to play against each of its neighbors on the network. Thereafter, its payoff is given by the average payoff per move and neighbor.

The strategy space of a player $`i`$ consists of up to 8 pure strategies $`S_i\{0,1,2,3,4,5,6,7\}`$ (cf. Tab. 2 for definition of the strategies). The pure-strategy space of the game is $`S=\underset{iI}{\times }S_i`$ with the set of the players $`I=\{0,1,\mathrm{},n\}`$. The (pure strategy) payoff function $`\pi _i:S`$ does not depend on the whole pure strategy profile $`s=(s_1,\mathrm{},s_n)`$ but only on the strategies of the neighboring nodes $`\pi _i=\pi _i(s_i,\text{neigh}(i))`$ and, of course, on the payoff matrix of the Prisoner’s dilemma game. Here, the set of the neighbors of a node $`i`$ is denoted $`\text{neigh}(i)`$. With $`\pi (s)=(\pi _1(s),\mathrm{},\pi _n(s))`$ the above defined game $`(S,I,\pi )`$ is a finite normal-form game. Such games in general have at least one Nash equilibrium Fudenberg and Tirole (1998). Here, only pure strategies will be considered neglecting possible mixed-strategy equilibria. In this setting, $`s_\text{D}=(0,\mathrm{},0)`$ and $`s_{\text{TFT}}=(6,\mathrm{},6)`$ are Nash equilibria for any payoff matrix $`A`$ obeying (2) and for a sufficiently high number of moves, which can be easily verified. The former equilibrium consists of players always defecting, whereas in the latter state each player repeats its opponent’s last move (Tit-For-Tat) with a cooperative opening move (cf. Tab. 2 ).

### II.2 Co-evolutionary dynamics

Let us now introduce mutations of a player’s strategy. The lookup table determining the strategy is viewed as a bit string of length $`2^m+1`$, where $`m`$ is the memory length as defined above. This bit-string will then be mutated during the iteration of the game.

At the beginning, a random network with a given mean degree $`k`$ is generated 111The resulting number of links $`L=kN/21`$ (rounded to an integer value) is distributed randomly with equal probability among all $`N(N1)/2`$ possible pairs of nodes leading to a constant probability $`k/(N1)`$ that an arbitrary pair of nodes is connected.. The strategies are assigned randomly, too. All agents play against each of their neighbors initially to update their payoffs. Thereafter, the following steps are iterated: (i) One agent $`i`$ is chosen randomly and its strategy is mutated from $`s_i`$ to a strategy $`s_i^{}S_i`$ picked out at random. (ii) The mutated agent plays again against its neighbors and its payoff is compared to the former result. The mutation is accepted in the case of a higher payoff, $`\pi _i(s_1,s_2,\mathrm{},s_i^{},\mathrm{},s_N)>\pi _i(s_1,s_2,\mathrm{},s_i,\mathrm{},s_N)`$ , and the payoffs of all neighbors are also updated. This corresponds to the assumptions that accepting any mutation is combined with some costs to the player and that mutations occur on a time scale slower than the time scale of the game. Iteration of this process leads to a stationary state with a fixed strategy distribution. In the stationary state, no agent can improve its payoff by changing its strategy whereas the other players’ strategies remain unchanged. This state corresponds to the game theoretic Nash equilibrium Nash (1950); Fudenberg and Tirole (1998). Note that, for its decisions, no more information than a player’s own payoff is required.

## III Perturbations and Avalanches

One essential property of evolutionary games is given by the equilibria or the evolutionarily stable states. All stationary states of the game are Nash equilibria. An interesting question is the stability of these equilibria against perturbations. In the following, we will study the dynamics of avalanches of mutation events following a perturbation of the Nash equilibrium. After the system has reached a stationary state, a new strategy is assigned to a random player. The insertion of a suboptimal strategy offers new opportunities for mutations to the perturbed player itself and to its neighbors. Since players are updated randomly, a perturbation leads to an avalanche of mutations until a stationary state is reached again. One quantity of interest is the avalanche size $`M`$ given by the number of mutations necessary to reestablish the equilibrium and its dependence on the payoff matrix $`A`$. We will first discuss the numerical results for the case of players on a random network. The second part of this section deals with a Prisoner’s dilemma on a ring, which will be the starting point for the theoretical treatment in the next two sections.

In the case of sparsely connected random networks, one observes three distinct regimes of the avalanche dynamics, with the temptation to defect $`a_{21}`$ as control parameter. For small temptations, $`3<a_{21}4`$, a subcritical regime occurs where large avalanches are suppressed exponentially (Fig. 1). For $`a_{21}>4`$, critical behavior occurs with avalanche sizes distributed according to a power law $`P(M)M^\gamma `$ with the scaling exponent $`\gamma =1.39\pm 0.10`$ (Fig. 2) and a cutoff scaling linearly with system size $`N`$. This critical regime is followed by a supercritical regime for $`4.70a_{21}<6`$ with an enhanced probability of very large events (Fig. 3).

Thus, above a critical value of the temptation to defect $`a_{21}^c=4`$, small perturbations of the system lead to long lasting avalanches that affect all players of the whole system with a mean avalanche size that diverges in the thermodynamic limit. The transition from a regime with small avalanches to a critical one with system-wide avalanches is robust in case of moderate changes of the strategy space $`S`$ and the mean degree $`k`$. It also occurs for smaller strategy spaces $`S_i`$ with $`\text{card}(S_i)5`$ and $`\{0,6,7\}S_i`$. The qualitative behavior remains even for indefinitely iterated games or with a very different payoff matrix 222 Since $`\widehat{a}_{12}=\widehat{a}_{22}`$ the game defined by $`\widehat{A}`$ is not a Prisoner’s dilemma in a strict sense (cf. Eqn. 2). Furthermore, getting the same payoff in case of a defecting opponent leads to coexistence of cooperative and defective domains in the subcritical regime. , which is sometimes used in the context of the Prisoner’s dilemma

$$\widehat{A}=\left(\begin{array}{cc}1& 0\\ \widehat{a}_{21}& 0\end{array}\right).$$
(3)
With $`S_i=\{0,7\}`$ and $`\widehat{A}`$ but quite different evolutionary dynamics, Lim, Chem, and Jayaprakash found critical avalanches on a two-dimensional square lattice, too Lim et al. (2002).

The different regimes of relaxation dynamics can be explained by a closer look on the structure of the Nash equilibria involved (Sec. IV) as well as on the relaxation mechanism which is given by a confined branching process (Sec. V). Before we start these considerations in the next two sections, we briefly discuss the case of a co-evolutionary Prisoner’s dilemma on a ring, i.e. on the regular network where each player’s degree is exactly $`k_i=2`$. Although this is not a very reasonable model for real spatially extended systems, it will give some useful insights and will allow us to calculate some properties of the spatially extended game analytically.

Like for random networks, subcritical, critical, and supercritical regimes occur, with subcritical avalanche distributions in the range of $`3<a_{21}4`$ and supercritical behavior for $`4.12a_{21}<6`$. However, this time the critical avalanche distribution has a scaling exponent of $`\gamma =1.04\pm 0.05`$ which significantly differs from the exponent obtained for random networks (Fig. 4).

## IV Nash equilibria and their dependence on the payoff matrix

The set of possible stationary states of the co-evolutionary Prisoner’s dilemma is the set of Nash equilibria which, as has been shown above, contains for all $`a_{21}(3,6)`$ the defective equilibrium $`s_\text{D}=(0,\mathrm{},0)`$ and the Tit-For-Tat equilibrium $`s_{\text{TFT}}=(6,\mathrm{},6)`$. One can also consider the macrostates of the system corresponding to the aggregated behavior of the agents. Identifying cooperative moves with “spin up” and defecting moves with “spin down” the macroscopic behavior at one instant of time is the magnetization of the system. Thus, the strategy profile $`s_\text{D}`$ of all agents playing strategy $`0`$ corresponds to the minimal magnetization $`1`$ whereas the Tit-For-Tat equilibrium $`s_{\text{TFT}}`$ leads to the maximal magnetization $`+1`$.

### IV.1 Equilibria on rings and random networks

Starting with the experimental findings for a ring topology, one observes three regimes in terms of adopted equilibria which are exactly matched by the three different regimes of avalanche dynamics. In the subcritical regime, the stationary states are a mixture of the strategies $`6`$ and $`7`$, i.e. generous-Tit-For-Tat players and unconditioned cooperators, respectively (Tab. 2). Only generous Tit-For-Tat prevails in the critical regime. With the onset of supercritical behavior, the defective equilibrium $`s_\text{D}`$ turns up. Its fraction of the equilibria reached by the game grows very fast with further increasing temptation to defect. The first transition can be explained by a simple calculation of the Nash equilibria. With $`a_{21}4`$, the cooperative macrostate is degenerated in many Nash equilibria since unconditioned cooperators are stabilized by neighbors with the strategy gTFT. Consider a player $`i`$ with its neighbors playing $`s_{i1}=7`$ (i.e. generous cooperate or gC) and $`s_{i+1}=6`$ (gTFT). Then player $`i`$ has to find a strategy being a compromise between exploiting the cooperator at $`i1`$ and maintaining cooperation with its other neighbor, the smarter Tit-For-Tat player at $`i+1`$. However, for $`a_{21}4`$, there is no such strategy yielding a better payoff than gTFT or even gC. This stabilization of the credulously cooperating agents gives rise to a degeneracy of the cooperative macrostate in many different strategy profiles that are Nash equilibria, diverging faster than $`2^{2N/3}`$ with the size of the ring. On the other hand, if $`a_{21}>4`$ there always exists such a compromise strategy and the degeneracy vanishes. That means that below the critical value $`a_{21}^c=4`$ the macrostate with magnetization $`+1`$ is strongly degenerated in many Nash equilibria whereas above $`a_{21}^c`$ there is only one Nash equilibrium with maximal magnetization left regardless of system size ($`s_{\text{TFT}}`$). The other macrostate with minimal magnetization is never degenerated since $`s_\text{D}`$ is the only Nash equilibrium that leads to such a defective macrostate. In case of regular lattices with different numbers of next neighbors $`k`$, the critical value $`a_{21}^c`$ is given by

$$a_{21}^c(k)=4\frac{2k+1}{k+3}.$$
(4)
For example, in the case of a two-dimensional lattice with periodic boundary conditions and a von-Neumann neighborhood, we find subcritical, critical, and supercritical behavior, with $`a_{21}^c5.14`$ and a critical exponent of $`\gamma =1.3\pm 0.2`$. Of course, for every payoff matrix $`A`$ satisfying (2) exists a finite value for the number of next neighbors $`k`$ above which the cooperative macrostate is always degenerated. However, $`a_{11}`$ and $`a_{21}`$ can be adjusted to increase this border to arbitrarily large values. Nonetheless, there is a reason why, for every payoff matrix, true critical behavior can only take place in sparsely connected networks which will be discussed in the next section. Why are only cooperative equilibria observed in both the subcritical and the critical range of $`a_{21}`$? Looking closer at the way to the equilibrium, there are $`(\sigma ^2(\sigma +1)(\sigma 1))/2`$ transition probabilities for a player $`i`$ changing its strategy, with $`\sigma :=\text{card}(S_i)`$. When increasing the temptation to defect $`a_{21}`$, some of these rules change from $`0`$ to a finite value and the respective inverse rule vice versa. At the transition to the supercritical regime, where the defective equilibrium is reached for the first time, exactly those rules change which govern the stability of the border between cooperative and defective domains. Below that threshold, the cooperative domains grow and above it defecting strategies can spread.

The situation is slightly different on a random network. Since there are always some nodes with a degree higher than the mean degree $`k`$, a small degeneracy of the cooperative macrostate can exist even for $`a_{21}>a_{21}^c(k)`$. Moreover, disconnected compounds may be in different equilibria at the same time. The highly connected nodes stabilize the cooperative equilibrium so that even for $`a_{21}6`$ cooperating strategies predominate. Yet these degeneracy does not compensate for the loss of cooperative equilibrium profiles for temptations larger than $`a_{21}^c`$ which is the reason for the transition from subcritical to critical behavior. The supercritical phase is again caused by the change of transition rules leading to increasing defective domains with their growth hindered by highly connected cooperative nodes.

### IV.2 Nash equilibria and ESS

As we are dealing with an evolutionary game, the question arises if any of the Nash equilibria is also evolutionarily stable. A strategy profile is called an evolutionarily stable state (ESS) if it is stable against the insertion of a small but finite fraction of mutants playing a different strategy Weibull (1996); Hofbauer and Sigmund (1998). Therefore, an ESS is a strict Nash equilibrium or a non-strict Nash equilibrium with the additional condition that other best replies play worse against themselves than the ESS strategy against them. Note that this concept is formulated for two-person games where two players encounter each other by chance. In this sense, both Nash equilibria $`s_\text{D}=(0,\mathrm{},0)`$ and $`s_{\text{TFT}}=(6,\mathrm{},6)`$ are no ESS for $`S_i=\{0,1,2,3,4,5,6,7\}`$ since other best replies score equally well as the equilibrium strategy (strategies $`2`$ and $`7`$, respectively). With respect to strategy spaces reduced by $`2`$ or $`7`$, the respective (now strict) Nash equilibrium becomes an ESS. But does this concept of stability apply to spatially extended games? Many approaches to evolutionary stability lead to the equivalence of ESS and strict Nash equilibria. So one may conjecture that in a game with $`S_i=\{0,1,5,6,7\}`$ the profile $`s_\text{D}`$ should be an ESS since it is a strict Nash equilibrium. Yet, as the experiments show, a small perturbation can cause the system to change from the strict Nash equilibrium $`s_\text{D}`$ to the non-strict cooperative equilibrium $`s_{\text{TFT}}`$. Thus, when applying the notion of evolutionary stability to spatially extended systems, one has to keep in mind that things may be different here since it is the local surrounding that decides whether an invader will overthrow the incumbent strategy or will fail instead.

## V Branching processes as a model of the relaxation process

Having understood the structure of the Nash equilibria and their connection to the transition between different regimes of avalanche dynamics, the question remains how to explain the distinct form of the probability distributions $`P(M)`$ and in particular the scaling exponents of the critical regimes. In fact, the relaxation process can be described by a type of branching process which very well predicts the scaling exponents for the different topologies as well as the subcritical and supercritical avalanche distributions.

### V.1 The Galton-Watson process

The starting point is a simple branching process, also known as Galton-Watson process, which will be reformulated in terms of mutated agents giving rise to future mutations of other players. Let $`Z_n`$ be the number of mutated players in the $`n`$-th generation. Each mutated player can cause other mutations in the next generation, with the probability $`p_m`$ that its mutation is succeeded by $`m`$ mutations in the next generation. The stochastic process $`(Z_n)_{n_0}`$ is called a branching process of the Galton-Watson type. Note that the number of generations constitutes a time scale completely different from the time scale of the game where at each instant of time one player is chosen to mutate its strategy. With the initial condition $`Z_0=1`$ and the total progeny $`Z:=_{i=0}^{\mathrm{}}Z_i`$, the quantity of interest is the distribution $`P(Z=r)`$ of total progeny or, in other words, the avalanche size. So far there is no spatial constraint to this process, i.e. $`Z_i`$ is not bounded by the system size, and mutations independently give birth to new mutations. The probability $`p_m`$ that a mutation of a player with $`k`$ neighbors will be followed by $`m`$ mutations in the next generation is given by

$$p_m=\left(\genfrac{}{}{0pt}{}{k+1}{m}\right)\alpha ^m(1\alpha )^{k+1m},$$
(5)
being the simplest choice if a player’s mutation can only affect its neighborhood including itself. Using generating functions Harris (1963) we calculate $`P(Z=r)`$ for this special Galton-Watson process with the same $`k_i=k`$ for all players to

$$P(Z=r)=\frac{\alpha 1}{\alpha k}\sqrt{\frac{k+1}{2\pi k}}\left(\frac{k^k}{\alpha (1\alpha )^k(k+1)^{k+1}}\right)^rr^{3/2}.$$
(6)
It is useful to introduce the mean number of a mutation’s “children” $`\overline{m}=\alpha (k+1)=EZ_0`$ and approximate (6) for $`\overline{m}1`$

$$P(Z=r)=Cr^{3/2}e^{\frac{r}{r_0}},$$
(7)
with

$$r_0=\frac{k+1}{2k}\frac{1}{(1\overline{m})^2},$$
(8)
and a constant

$$C=\frac{\overline{m}k1}{k\overline{m}}\sqrt{\frac{k+1}{2\pi k}}.$$
(9)
If the expectation of the numbers of descendants approaches one, i.e $`\overline{m}1`$, the exponential cutoff diverges with $`(1\overline{m})^2`$. The process becomes critical with a scale-free avalanche distribution $`P(Z=r)r^{3/2}`$. If $`\overline{m}>1`$, the probability is finite that $`Z`$ does not converge at all Jagers (1975). The branching process described above, which has no spatial constraints, is characterized by a subcritical, critical, and supercritical regime of its avalanche dynamics. Although this is very similar to the IPD on a random network, in the case of a ring it yields a wrong scaling exponent of $`\gamma =3/2`$. Such behavior could be gained equally well from a random walk of the number of mutated sites with drift to a reflecting boundary. In the following, we will show that it is the restriction of the branching process to the network topology that completely explains the dynamics and leads to the correct scaling exponents.

### V.2 Confined branching processes

The confinement of the branching process leads to two effects. First, $`Z_n`$ will be bounded by the system size $`N`$; second, the mutation events caused by mutated players are no longer stochastically independent. We will denote a branching process as confined or restricted to a network (i) if there exists a one-to-one mapping of players and nodes and (ii) if a mutated player can only give birth to mutations in its neighborhood including itself (Fig. 5 A).

This corresponds to the fact that if a player changes its strategy only the payoffs of its neighbors and of the player itself will be affected. We assume that each neighbor and the mutated site itself has the same probability $`\alpha `$ of mutation in the next generation. With the random variable $`X_\nu ^{(n)}`$ being $`1`$ if the player at node $`\nu `$ is mutated in generation $`n`$ and $`0`$ otherwise, the confined process $`(Z_n^{})_{n_0}`$ is defined by

$$Z_n^{}=\underset{\nu =1}{\overset{N}{}}X_\nu ^{(n)}.$$
(10)
The probability of a mutation at site $`\nu `$ in generation $`n`$ is

$$P(X_\nu ^{(n)}=1)=1(1\alpha )^\lambda ,$$
(11)
with

$$\lambda =\underset{\mu \text{neigh}(\nu )\{\nu \}}{}X_\mu ^{(n1)}.$$
(12)
The confined branching process $`(Z_n^{})`$ can now be used to calculate the avalanche distributions of the spatially extended Prisoner’s dilemma numerically. Applying it to random networks, both the subcritical and supercritical distributions are matched well (dashed curves in Fig. 1 and Fig. 3). The distribution of the confined branching process agrees even better with the experimental data in the critical regime (dashed curves in Fig. 2). The same is true for the Prisoner’s dilemma on a ring (Fig. 4). In both critical cases, the branching process shows the correct finite-size scaling of the cutoff which is proportional to the system size. Note that the critical regimes of the game have different scaling exponents due to network topology which are both correctly obtained by the confined branching process. The critical exponents depend only on the topology rather than on the parameter $`\alpha `$ of the process. Therefore, the relaxation mechanism of the spatially extended co-evolutionary Prisoner’s dilemma is a confined branching process.

Mean-field approaches can be applied successfully to explain the parameter $`\alpha `$ of the confined branching process in the subcritical and critical regime (Tab. 3). To calculate a mean-field approximation $`\alpha _{\text{mf1}}`$ of the branching parameter, the transition probabilities of a mutated agent’s neighbors are determined using a random neighborhood for both the player and its neighbors. The structure of the game is taken into account only by assuming that the stable strategies are absorbing states. A second approach is to average the transition probabilities over game realizations numerically, yielding $`\alpha _{\text{mf2}}`$. Both values, $`\alpha _{\text{mf1}}`$ and $`\alpha _{\text{mf2}}`$, agree well with the parameter $`\alpha `$ obtained from the avalanche distributions of the subcritical and critical regime. This corresponds to the explanation that this transition occurs solely because of the change in the degeneracy of the cooperative macrostate at the critical value $`a_{21}^c`$. The supercritical case is not matched by the mean-field approaches which may be due to the fact that here the dynamics are governed by local effects, i.e. the competitive growth at the boundaries between cooperative and defective domains. The dynamics on a ring topology can be explained by a similar mean-field approach, too, if one assumes that the effective maximal number of a player’s descendants is approximately two and not three. This reduction of potential progeny is caused by the strong overlap of the neighborhoods in this regular lattice (Fig. 5 B).

Although the definition of the confined Galton-Watson process is quite intuitive and simple, its analytical treatment is not. The reason is that mutation events has become dependent on each other. Two mutations can affect the same site in the next generation (Fig. 5 B) leading to dependent recursive equations (11, 12) for the mutation probability. With the simplification that the $`Z_n^{}`$ mutated sites of generation $`n`$ are randomly distributed over the network, one can shed some light on the critical behavior of the confined branching process. The conditional expectations of the number of mutated players are with this assumption

$$E(Z_{n+1}^{}|Z_n^{})=\overline{m}Z_n^{}(1+\xi )$$
(13)
with

$$\xi =\left(\overline{m}\frac{Z_n^{}}{N}\right)^1\left[1\left(1\frac{1}{k+1}\overline{m}\frac{Z_n^{}}{N}\right)^{k+1}\overline{m}\frac{Z_n^{}}{N}\right].$$
(14)
If $`\xi 1`$ and $`\overline{m}1`$ the confined process approximately is a martingale for all values of $`Z_n^{}`$ and should show critical behavior. For $`\overline{m}1`$ the avalanche dynamics are subcritical as the process becomes a supermartingale. With $`\overline{m}1`$ obviously resulting in supercritical dynamics, the remaining case of interest is $`\overline{m}1`$. In the event of highly connected networks with $`k1`$ the correction $`\xi `$ is of the order $`1`$ suppressing large avalanches. Thus critical avalanche dynamics are expected only for sparsely connected networks, for too strong dependencies of mutation events lead to either subcritical or supercritical distributions of avalanche sizes.

## VI Conclusions

In this paper, we have introduced a spatially extended Prisoner’s dilemma game with co-evolutionary dynamics that lead to Nash equilibria as stationary states. We have shown that critical avalanche dynamics are characteristic for a broad range of these games. The observed intermittent evolution with sudden avalanches of activity is reminiscent of self-organized criticality Bak (1996); Jensen (1998). Depending on the payoff matrix, subcritical, critical, and supercritical regimes can be observed. Calculating the Nash equilibria and introducing a confined branching process, we were able to quantitatively explain the critical value of the control parameter, i.e. the temptation to defect, and the avalanche distributions. Therefore, investigations on the spatially extended Prisoner’s dilemma, which has become a widely used toy model for the emergence of cooperation, have to take into account the stability of possible equilibria depending on chosen payoff matrix, strategy space, and topology. Complex behavior should only be found for subcritical or critical dynamics whereas in the supercritical regime small perturbations will totally mix up the whole system preventing the evolution of local structures. The results on the stability of the Nash equilibria and their connection to evolutionarily stable states indicate that the concept of equilibrium, originating from classical mechanics and brought into the fields of game theory and evolution Ekeland (2002), has to be further specified to take into account co-evolution on networks and other spatial structures.

###### Acknowledgements.

We thank Christoph Hauert, Thomas Lux, and Uwe Rösler for useful discussions and comments. H. E. gratefully acknowledges support by the Studienstiftung des deutschen Volkes (German National Merit Foundation).


text:
# Search for neutrinoless 𝜏 decays involving the 𝐾_𝑆⁰ meson

(August 6, 2002)

## Abstract

We have searched for lepton flavor violating decays of the $`\tau `$ lepton with one or two $`K_S^0`$ mesons in the final state. The data used in the search were collected with the CLEO II and II.V detectors at the Cornell Electron Storage Ring (CESR) and correspond to an integrated luminosity of 13.9 $`fb^1`$ at the $`\mathrm{{\rm Y}}(4S)`$ resonance. No evidence for signals were found, therefore we have set 90% confidence level (C.L.) upper limits on the branching fractions $`(\tau ^{}e^{}K_S^0)<9.1\times 10^7`$, $`(\tau ^{}\mu ^{}K_S^0)<9.5\times 10^7`$, $`(\tau ^{}e^{}K_S^0K_S^0)<2.2\times 10^6`$, and $`(\tau ^{}\mu ^{}K_S^0K_S^0)<3.4\times 10^6`$. These represent significantly improved upper limits on the two-body decays and first upper limits on the three-body decays.

preprint:   CLNS 02/1795 CLEO 02-12

S. Chen,1 J. W. Hinson,1 J. Lee,1 D. H. Miller,1 V. Pavlunin,1 E. I. Shibata,1 I. P. J. Shipsey,1 D. Cronin-Hennessy,2 A.L. Lyon,2 C. S. Park,2 W. Park,2 E. H. Thorndike,2 T. E. Coan,3 Y. S. Gao,3 F. Liu,3 Y. Maravin,3 R. Stroynowski,3 M. Artuso,4 C. Boulahouache,4 K. Bukin,4 E. Dambasuren,4 K. Khroustalev,4 R. Mountain,4 R. Nandakumar,4 T. Skwarnicki,4 S. Stone,4 J.C. Wang,4 A. H. Mahmood,5 S. E. Csorna,6 I. Danko,6 G. Bonvicini,7 D. Cinabro,7 M. Dubrovin,7 S. McGee,7 A. Bornheim,8 E. Lipeles,8 S. P. Pappas,8 A. Shapiro,8 W. M. Sun,8 A. J. Weinstein,8 R. Mahapatra,9 R. A. Briere,10 G. P. Chen,10 T. Ferguson,10 G. Tatishvili,10 H. Vogel,10 N. E. Adam,11 J. P. Alexander,11 K. Berkelman,11 V. Boisvert,11 D. G. Cassel,11 P. S. Drell,11 J. E. Duboscq,11 K. M. Ecklund,11 R. Ehrlich,11 R. S. Galik,11 L. Gibbons,11 B. Gittelman,11 S. W. Gray,11 D. L. Hartill,11 B. K. Heltsley,11 L. Hsu,11 C. D. Jones,11 J. Kandaswamy,11 D. L. Kreinick,11 A. Magerkurth,11 H. Mahlke-Krüger,11 T. O. Meyer,11 N. B. Mistry,11 E. Nordberg,11 J. R. Patterson,11 D. Peterson,11 J. Pivarski,11 D. Riley,11 A. J. Sadoff,11 H. Schwarthoff,11 M. R. Shepherd,11 J. G. Thayer,11 D. Urner,11 G. Viehhauser,11 A. Warburton,11 M. Weinberger,11 S. B. Athar,12 P. Avery,12 L. Breva-Newell,12 V. Potlia,12 H. Stoeck,12 J. Yelton,12 G. Brandenburg,13 D. Y.-J. Kim,13 R. Wilson,13 K. Benslama,14 B. I. Eisenstein,14 J. Ernst,14 G. D. Gollin,14 R. M. Hans,14 I. Karliner,14 N. Lowrey,14 C. Plager,14 C. Sedlack,14 M. Selen,14 J. J. Thaler,14 J. Williams,14 K. W. Edwards,15 R. Ammar,16 D. Besson,16 X. Zhao,16 S. Anderson,17 V. V. Frolov,17 Y. Kubota,17 S. J. Lee,17 S. Z. Li,17 R. Poling,17 A. Smith,17 C. J. Stepaniak,17 J. Urheim,17 Z. Metreveli,18 K.K. Seth,18 A. Tomaradze,18 P. Zweber,18 S. Ahmed,19 M. S. Alam,19 L. Jian,19 M. Saleem,19 F. Wappler,19 K. Arms,20 E. Eckhart,20 K. K. Gan,20 C. Gwon,20 T. Hart,20 K. Honscheid,20 D. Hufnagel,20 H. Kagan,20 R. Kass,20 C. Morris,20 T. K. Pedlar,20 J. B. Thayer,20 E. von Toerne,20 T. Wilksen,20 M. M. Zoeller,20 H. Muramatsu,21 S. J. Richichi,21 H. Severini,21 P. Skubic,21 S.A. Dytman,22 J.A. Mueller,22 S. Nam,22 and V. Savinov22

1Purdue University, West Lafayette, Indiana 47907

2University of Rochester, Rochester, New York 14627

3Southern Methodist University, Dallas, Texas 75275

4Syracuse University, Syracuse, New York 13244

5University of Texas - Pan American, Edinburg, Texas 78539

6Vanderbilt University, Nashville, Tennessee 37235

7Wayne State University, Detroit, Michigan 48202

8California Institute of Technology, Pasadena, California 91125

9University of California, Santa Barbara, California 93106

10Carnegie Mellon University, Pittsburgh, Pennsylvania 15213

11Cornell University, Ithaca, New York 14853

12University of Florida, Gainesville, Florida 32611

13Harvard University, Cambridge, Massachusetts 02138

14University of Illinois, Urbana-Champaign, Illinois 61801

15Carleton University, Ottawa, Ontario, Canada K1S 5B6

and the Institute of Particle Physics, Canada M5S 1A7

16University of Kansas, Lawrence, Kansas 66045

17University of Minnesota, Minneapolis, Minnesota 55455

18Northwestern University, Evanston, Illinois 60208

19State University of New York at Albany, Albany, New York 12222

20Ohio State University, Columbus, Ohio 43210

21University of Oklahoma, Norman, Oklahoma 73019

22University of Pittsburgh, Pittsburgh, Pennsylvania 15260

In physics, all fundamental conservation laws are expected to have associated symmetries. Lepton flavor conservation, however, is an experimentally observed phenomena with no associated symmetry in the Standard Model. Lepton flavor violation (LFV) is expected in many extensions of the Standard Model such as lepto-quark, supersymmetry, superstring, left-right symmetric models, and models that include heavy neutral leptons Cvetic02 . Experimentally, both Super Kamiokande Fukuda98  and SNO SNO01  observe neutrino oscillation, which may imply LFV in the neutrino sector; therefore LFV is expected to occur in charged lepton decay at some branching fraction, albeit very small. The $`\tau `$ lepton provides a clean laboratory for such searches. Ilakovac Ilakovac95  has calculated upper limits on branching fractions for many neutrinoless LFV modes within a model involving heavy Dirac neutrinos. The branching fractions depend on the heavy neutrino masses and mixings. For the decays $`\tau ^{}\mathrm{}^{}K_S^0`$ note1 , the branching fractions are of $`𝒪(10^{16})`$, where $`\mathrm{}`$ can be $`e`$ or $`\mu `$. For the decays $`\tau ^{}\mathrm{}^{}K_S^0K_S^0`$ or $`\tau ^{}\mathrm{}^{}K^+K^{}`$, the branching fractions are of $`𝒪(10^7)`$. The decays with two kaons in the final state are therefore of particular experimental interest. Previous published upper limits on the branching fractions for the decays $`\tau ^{}\mathrm{}^{}K_S^0`$ are of $`𝒪(10^4)`$ Hayes82 . There are no previous results for the decays $`\tau ^{}\mathrm{}^{}K_S^0K_S^0`$. In this paper, we present the results of a search for the decays into one lepton and one or two $`K_S^0`$ mesons, with the $`K_S^0`$ decaying into two charged pions.

The data used in this analysis were collected using the CLEO detector CLEOII  from $`e^+e^{}`$ collisions at the Cornell Electron Storage Ring (CESR) at a center-of-mass energy $`\sqrt{s}10.6`$ GeV. The total integrated luminosity of the data sample is 13.9 $`fb^1`$ corresponding to the production of $`N_{\tau \tau }=1.27\times 10^7`$ $`\tau ^+\tau ^{}`$ events. The CLEO detector is a general purpose spectrometer with excellent charged particle and shower energy detection. The momenta of charged particles are measured with three drift chambers between 5 and 90 cm from the $`e^+e^{}`$ interaction point (IP), with a total of 67 layers. For $``$63% of the data collected, the innermost tracking chamber was replaced by a three-layer silicon vertex detector CLEOIIV . The specific ionization ($`dE/dx`$) of charged particles is also measured in the main drift chamber. The tracking system is surrounded by a scintillation time-of-flight system and a CsI(T1) calorimeter with 7800 crystals. These detector systems are installed inside a superconducting solenoidal magnet (1.5 T), surrounded by an iron return yoke instrumented with proportional tube chambers for muon identification.

The $`\tau ^+\tau ^{}`$ candidate events must contain four or six charged tracks with zero net charge. The polar angle $`\theta `$ of each track with respect to the beam must satisfy $`|\mathrm{cos}\theta |<0.90`$. To reject beam-gas events, the distance of closest approach of each non-$`K_S^0`$ track to the IP must be within 0.5 cm transverse to the beam and 5 cm along the beam direction. Photons are defined as energy clusters in the calorimeter with at least 60 MeV in the barrel ($`|\mathrm{cos}\theta |<0.80`$) or 100 MeV in the endcap ($`0.80<|\mathrm{cos}\theta |<0.95`$). We further require every photon to be separated from the projection of any charged track by at least 30 cm unless its energy is greater than 300 MeV. In order to diminish QED background such as radiative Bhabha and $`\mu `$-pair events with photon conversion, we require each event to have total energy less than 95% of the center-of-mass energy. This requirement rejects most of the QED background, while incurring a small loss in detection efficiency.

We divide each event into two hemispheres (signal and tag), one containing one charged track and the other containing three or five charged tracks, using the plane perpendicular to the thrust axis Farhi77 . The thrust axis is calculated from both charged tracks and photons. The invariant mass of the tag hemisphere must be less than the $`\tau `$ mass, $`M_\tau =1.777`$ GeV/$`c^2`$ PDG00 . The signal hemisphere must contain an electron or a muon and one or two $`K_S^0`$ mesons. The electron candidate must have shower energy to momentum ratio in the range, $`0.85<E/p<1.10`$, and when available, the specific ionization lost must be consistent with that expected for an electron. The muon candidate must penetrate at least three absorbtion lengths of iron. The $`K_S^0`$ candidate is reconstructed in the $`\pi ^+\pi ^{}`$ final state with a detached vertex, and the invariant mass must be within approximately three standard deviations of the nominal mass, $`485<m_{\pi \pi }<510`$ MeV/$`c^2`$, as determined from a signal Monte Carlo simulation (see below). In order to diminish radiative Bhabha and $`\mu `$-pair events further, neither pion should be consistent with identification as an electron.

Since there is no neutrino in the signal hemisphere while there is at least one neutrino undetected in the tag hemisphere, the missing momentum of the event must point toward the tag hemisphere, $`0<\mathrm{cos}\theta _{tagmissing}<1.0`$. In order to suppress the background from radiative Bhabha and $`\mu `$-pair events, the direction of the missing momentum of the event is required to satisfy $`|\mathrm{cos}\theta _{missing}|<0.90`$. For the decay $`\tau ^{}e^{}K_S^0`$, $`\mathrm{cos}\theta _{tagmissing}`$ is further restricted to be less than 0.99 to reduce the radiative backgrounds. This corresponds to the minimum ratio of background to detection efficiency. The background is estimated from the sidebands in the invariant mass vs. total energy distribution of the decay candidates.

To search for decay candidates, we select $`\tau `$ candidates with invariant mass and total energy consistent with the expectations. The following kinematic variables are used to select the candidate events:

$`\mathrm{\Delta }E`$ $`=`$ $`EE_{beam}`$
$`\mathrm{\Delta }M`$ $`=`$ $`MM_\tau ,`$
where $`E_{beam}`$ is the beam energy, and $`E`$ and $`M`$ are the reconstructed $`\tau `$ energy and mass. The $`\mathrm{\Delta }E`$ vs. $`\mathrm{\Delta }M`$ distributions of the decay candidates in the data and signal Monte Carlo samples are shown in Figs. 1 and 2. The center of the signal region is slightly shifted from zero in the $`\mathrm{\Delta }E`$ vs. $`\mathrm{\Delta }M`$ plane to account for initial state radiation and shower leakage. The signal region is defined as the area within three standard deviations ($`\sigma `$) of the expectation for both kinematic variables, as determined from a signal Monte Carlo simulation.

In the Monte Carlo simulation, one $`\tau `$ lepton decays according to two- or three-body phase space for the mode of interest, and the other $`\tau `$ lepton decays generically according to the $`KORALB/TAUOLA`$ $`\tau `$ event generator taugen85-91 . The phase space model is appropriate for an unpolarized tau. If the Lorentz structure of the neutrinoless decay is V-A, as in the Standard Model, correlations between the spins of the two $`\tau `$’s in the event will lead to slightly higher detection efficiency than phase space, while V+A decays will lead to lower detection efficiency. The detector response is simulated using the $`GEANT`$ program Brun-GEANT . The estimated detection efficiencies ($`ϵ`$eff  are summarized in Table 1.

The upper limit on the branching fraction is related to the upper limit $`\lambda `$ on the number of signal events by

$`={\displaystyle \frac{\lambda }{2ϵN_{\tau \tau }_1(_{K_S^0\pi ^+\pi ^{}})^n}},`$
where $`_1=(84.71\pm 0.13)\%`$ is the inclusive 1-prong branching fraction PDG00 , $`_{K_S^0\pi ^+\pi ^{}}=(68.61\pm 0.28)\%`$ is the branching fraction for $`K_S^0`$ to decay to two charged pions, and $`n`$ is the number of $`K_S^0`$ mesons in the final state. No candidate decays are observed, so we take $`\lambda `$ as 2.44 events in each mode at 90% confidence level, according to the frequentist method Feldman98 . The upper limits on the branching fractions with and without systematic uncertainties are shown in Table 1. The systematic uncertainties include the $`\tau ^+\tau ^{}`$ cross section (1%), luminosity (1%), track reconstruction efficiency (1% per charged track), $`K_S^0`$ detection efficiency (2% per $`K_S^0`$), lepton identification (1.5% for electron and 4% for muon), and the statistical uncertainties in the detection efficiencies due to limited Monte Carlo samples (1-2%).

Black et al. Black02  have analyzed the constraints on the new physics scale for dimension-six fermionic effective operators involving $`\tau \mu `$ mixing, motivated by the observed $`\nu _\mu \nu _\tau `$ oscillation. The most stringent lower limits from exotic heavy quarks and $`\tau `$ decays on the physics scale of the operators involving quarks are $``$10 TeV. The new upper limit on $`(\tau ^{}\mu ^{}K_S^0)`$ presented in this paper yields a lower limit of 17.3 and 18.2 TeV for the axial vector and pseudoscalar operators, respectively.

In conclusion, we have searched for $`\tau `$ decays involving $`K_S^0`$ mesons that violate lepton flavor, but find no evidence for a signal. This results in improved upper limits for the decays $`\tau ^{}\mathrm{}^{}K_S^0`$ and first upper-limits for the decays $`\tau ^{}\mathrm{}^{}K_S^0K_S^0`$. The upper limits for the $`\tau ^{}\mathrm{}^{}K_S^0K_S^0`$ final states are more stringent than those found previously for $`\mathrm{}^{}K^+K^{}`$ Bliss98 .

We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. M. Selen thanks the PFF program of the NSF and the Research Corporation, and A.H. Mahmood thanks the Texas Advanced Research Program. This work was supported by the National Science Foundation, and the U.S. Department of Energy.