# Velocity statistics in excited granular media

## I Introduction

Granular systems are often treated statistically, since the large number of degrees of freedom and the complexity of interparticle forces limit analyses on the particle level. What are the statistical properties of particles in an excited granular medium? In many respects the dynamics of granular media strongly resemble the dynamics of ordinary fluids and solids. Yet there are fundamental differences between granular solids and fluids and their molecular counterparts that make the answer to this question intriguing . Ordinary temperature does not lead to measurable velocity fluctuations due to the large mass of granular particles. When granular particles are excited by an external energy source, their thermal energies are much smaller than their kinetic energies. During collisions between particles, frictional forces and deformation near the points of contact lead to dissipation of kinetic energy. Therefore it is not surprising to find that kinetic energy is transformed into thermal energy when particles interact in an excited granular medium .

For strongly excited granular media, assumptions similar to those in kinetic theory are often made . Among these assumptions are that velocity distributions are Gaussian and that the mean energy is shared equally among the various degrees of freedom. However, recent numerical and theoretical research indicates that excited inelastic hard spheres can exhibit non-Gaussian velocity distributions , though the predicted velocity distributions differ considerably from each other. Some thermodynamic descriptions of granular media make use of the concept of entropy , or separate the dissipative degrees of freedom from conservative ones .

Vibrated granular media have played a special role in efforts to understand the dynamics of granular materials, in part because vibration is a convenient method of replacing the energy lost to friction and inelasticity. A variety of novel phenomena have been discovered in the past fifteen years : heaping and convection rolls , standing and traveling waves , oscillons , and fluidization .

In this paper, we study the simple case of a single layer of particles with vertical excitation, varying the fractional coverage of this layer and the amplitude of the driving waveform. We determine the velocity distributions precisely and compare them to various functional forms, an issue that previous experiments have left unresolved for granular particles that are free to move in three dimensions. Quantitatively studying particles in three dimensions has only recently been achieved in rotated cylinders using NMR techniques; it yielded important information about the segregation process . By observing a three dimensional excited granular medium from above, we are able to focus on the shapes of the distributions of the horizontal velocity components. A combination of otherwise identical white and black particles allows us to track individual particles even for high fractional coverages at high excitations where particles frequently move over each other. We measure the dependence of the variances of the distributions on the excitation amplitude, frequency and the fractional coverage. Sharing of energy in mixtures of granular materials is also investigated.

In Sec. II, we discuss the background for our experiments. In Sec. III we present the experimental setup. In Sec. IV, we show that, in most cases we treated experimentally, the velocity distribution $`P(v)`$ deviates measurably from a Gaussian and can be described by $`P(v)\mathrm{exp}(|v/v_c|^{1.5})`$. We show precisely how the variance increases with driving amplitude and changes with fractional coverage $`c`$, where $`c=1`$ is the coverage for a compact crystal of beads. The energy sharing between particles of different types is also described here.

## II Background

A number of studies on excited granular media have addressed the issues of clustering and energy sharing. Numerical simulations and theoretical derivations have shown that the presence of inelasticity in granular flows can lead to the formation of clusters; as a consequence, equipartition of energy fails . Experiments performed in a horizontal two-dimensional layer were consistent with this predicted clustering effect . Vertical one-dimensional experiments and simulations were performed earlier ; a crossover from a condensed (clustered) to a fluidized state was found as a function of the driving acceleration, the number of beads, and the coefficient of restitution.

Deviations from equipartition due to clustering are straightforward to understand on physical grounds. Inelastic collisions imply a loss of energy each time a collision occurs. When particles begin to gather in a certain region of space, the rate of their collisions increases. The rate of energy loss for this group of clustered particles is thus greater, and the distribution of their velocities becomes narrower than that of the particles in less dense regions.

Clustering can affect the velocity statistics. Assuming a Gaussian velocity distribution for a nearly homogeneous granular medium, Puglisi et al.  predict a non-Gaussian velocity distribution due to clustering: a superposition of Gaussian velocity distributions with different widths. For inelastic particles this leads to high velocity tails in the velocity distribution which decrease more slowly than a Gaussian function but faster than an exponential. The velocity distributions obtained in a simulation by Taguchi and Takayasu  are power laws, resulting from clustering. These distributions have diverging variance; this calls into question the notion of a granular temperature, which is generally derived from a variance. In a two-dimensional simulation, Peng and Ohta  found that the velocity distributions deviate from Gaussian behavior under the influence of gravity unless $`g\delta h<<T`$, where $`\delta h`$ is the height of the region of observation. $`T`$ is the granular temperature (see e.g. Ref. ), defined as the variance of the velocity (minus the mean velocity)

$$T=<(v^2<v>^2)>.$$
(1)   
The high energy tails of velocity distributions in a homogeneous granular fluid were investigated theoretically by Esipov and Pöschel  for the unforced case, and by Noije and Ernst  for both the unforced and the heated case, based on the Enskog-Boltzmann equation. In the free case, i.e. without energy input into the system, the velocity distribution approaches an exponential at high velocities. When energy is added to the granular medium randomly and uniformly throughout the system, the high velocity approximation becomes $`P(v)\mathrm{exp}(|v|^{1.5})`$. In order to compare our experimental results to this prediction we have to assume that energy input into the measured horizontal motion occurs randomly.

While measurements of velocity statistics in an excited granular medium have been carried out , few measurements precise enough to distinguish between different functional forms of the velocity distribution exist to our knowledge. One exception is a recent study of clustering and ordering near a peak acceleration of $`a=1`$ g by Olafsen and Urbach , who found significant deviations from Gaussian distributions, at low and especially at high velocities, where the distributions become exponential.

Deviations from equipartition can occur for reasons other than clustering. For example, Knight and Woodcock  studied a vibrationally excited granular system theoretically (without gravity) and concluded that equipartition need not be observed at high amplitude of excitation due to the anisotropy of the energy source. A two dimensional system (a vertical Hele Shaw cell) was studied experimentally by Warr, Huntley, and Jacques . Velocity distributions, though roughly Gaussian, were checked and found to exhibit anisotropy between the vertical and horizontal motion: the horizontal velocity distribution was narrower than the vertical one. Grossman, Zhou, and Ben-Naim  considered a two-dimensional granular gas (without gravity) from a theoretical quasi-continuum point of view and found the density to be nonuniform and the velocity distributions to be asymmetric for “thermal” energy input from one side. McNamara and Luding  considered the sharing of energy between rotational and translational motion of the particles, and found a significant violation of equipartition.

The scaling of the granular kinetic energy with vibration amplitude is also of interest in connection with the experiments to be discussed in the present paper. It has been considered experimentally by Warr et al., numerically by Luding, Hermann, and Blumen , and theoretically by Kumaran in the nearly elastic limit of weak dissipation  and also by Huntley  in a simple model. The results of these different studies do not seem to be mutually consistent with each other, perhaps because different regimes were explored; the situation is unclear.

## III Experimental setup and methods

The experiments are conducted in a circular container of diameter 32 cm made of delrin. It is driven vertically with sinusoidal acceleration at a single frequency  using a VTS500 vibrator from Vibration Test Systems Inc. A computer controlled feedback loop keeps the vibration amplitude constant and reproducible. The frequency $`f`$ used in most experiments is $`100`$ Hz and the peak acceleration $`a`$ of the plate is in the range $`38g`$. The peak plate velocity $`v_p=a/2\pi f`$ lies between $`4.7`$ cm/s and $`14.1`$ cm/s. The particles are glass beads $`4`$ mm in diameter (from Jaygo Inc.) with fractional coverage $`c`$. A glass cover at $`2`$ cm height prevents the beads from escaping from the container. Collisions with the cover are rare for most experimental parameters but lead to measureable changes in the velocity distributions at the largest accelerations, if the coverage is low. Some charging of the glass beads is noticeable when the beads are at rest, but electrostatic forces are negligible at the range of accelerations investigated here.

Our objective is to measure the horizontal velocity distributions in an excited three dimensional granular medium for a large range of particle densities. However, particle tracking becomes increasingly difficult as particle tracks approach each other and cross frequently at large $`c`$, thus hindering reliable identification of the horizontal positions of individual particles. However, if only a modest number of particles are reflective, the frequency of collisions between them is small and nearly independent of coverage or excitation. By tracking these test particles we can also directly compare the velocity distributions for different coverage.

Tracking a subset of the particles is accomplished by using some white glass beads among black glass beads. Except for the color, the black and white beads have identical physical properties. Stainless steel beads of three different diameters replaced the white glass beads in some experiments. The material properties of all particles we used are listed in table I.

Images of an area $`16.74\times 15.70`$ cm at the center of the container are taken at a resolution of $`512\times 480`$ pixels using a fast camera (SR-500, Kodak Inc.) operated at $`250`$ or $`500`$ frames/s. At a vibration frequency of 100 Hz this ensures that images are taken at 5 different phases relative to the phase of the plate vibration yielding the average energy throughout the cycle. The images are analyzed using IDL (Research Systems Inc.) software. First each image is enhanced using a bandpass filter and thresholding to eliminate noise. The positions of all bright particles are found from the enhanced image by calculation of their centroid; this defines the particle positions reproducibly to within less than $`0.1`$ pixels. Small effects due to the finite pixel size are noticeable sometimes as slightly increased probabilities of particle displacements that are multiples of the pixel width ($`0.0327`$ cm in the physical system). This displacement corresponds to $`v=8.18`$ cm/s in most experiments. The high frame rate ensures that even the fastest beads move less than one particle diameter between images. This allows accurate tracking of all bright beads for all $`546`$ sequential frames (the maximum available with our camera) with a typical precision of $`\pm 2\%`$ for the velocity measurements. We project each step onto two perpendicular directions, and study the statistics of each velocity component. One concern was that as a black bead moves over a bright bead, the centroid position could move away from the center of the bead and possibly alter the measured velocity distribution. For a test we thus eliminated particles whose integrated greyscale intensity changed rapidly. We found that eliminating these points does not measurably alter the velocity distribution. However, if a collision occurs between frames, our measurements indicate the average of the velocity prior to and after a collision. The measured distribution of velocities will therefore probably be slightly closer to a Gaussian than the real distribution of velocities. To improve the data in this regard, it would be necessary to sample faster while retaining the same relative accuracy in the velocity measurement. It would therefore be necessary to use a faster camera, to zoom in closer to the sample and to extend the measurement over a significantly larger number of frames to obtain the same statistics.

## IV Experimental Results

### A Granular Temperature

Extracted particle tracks are shown in Fig. 1 for $`c=0.28`$ (half of the beads are white) and $`a=5`$g. For clarity, only tracks longer than $`200`$ images are shown, which eliminates some tracks close to the edge. A total of $`546`$ frames is acquired at $`250`$ frames/s (i.e. for $`2.18`$ seconds) with approximately $`200`$ tracked particles in each frame.

The particle velocities are determined from the particle displacement between consecutive frames. This does not always represent the true velocity of the particle though. If a collision occurs between frames, the apparent velocity will be lower than the true velocity. We can define an apparent temperature based on displacements along either horizontal coordinate denoted here as $`x`$; it depends on the time between frames $`\mathrm{\Delta }t`$:

$$T(\mathrm{\Delta }t)=<(x_j(t_k+\mathrm{\Delta }t)x_j(t_k))^2>_{j,k}/\mathrm{\Delta }t^2,$$
(2)   
where the average is taken over all particles (j) and frames (k). The particle tracks obtained at $`250`$ frames/s during an interval of $`2.18`$ s allow us to determine $`T(\mathrm{\Delta }t)`$ for $`1/250\mathrm{s}\mathrm{\Delta }t<<2.18\mathrm{s}`$. It is often useful to express the velocities in units of the peak plate velocity $`v_p`$, which yields a dimensionless temperature $`\stackrel{~}{T}(\mathrm{\Delta }t)=T(\mathrm{\Delta }t)/v_p^2`$. Fig. 2 shows the dimensionless temperature $`\stackrel{~}{T}(\mathrm{\Delta }t)`$ vs. frame rate ($`1/\mathrm{\Delta }t`$) at $`c=0.42`$ and $`f=100`$Hz for different accelerations. For small frame rates, i.e. large $`\mathrm{\Delta }t`$, $`\stackrel{~}{T}(\mathrm{\Delta }t)`$ increases approximately linearly with the frame rate $`\stackrel{~}{T}(\mathrm{\Delta }t)1/\mathrm{\Delta }t`$. This indicates that the particle motion may be described by an ordinary diffusion law when many collisions occur in the sampling interval $`\mathrm{\Delta }t`$. Assuming such a diffusion process, one expects that

$$\stackrel{~}{T}(\mathrm{\Delta }t)\frac{\stackrel{~}{T}\tau _c}{\mathrm{\Delta }t}(\mathrm{for}\mathrm{\Delta }t>>\tau _c).$$
(3)   
The dashed and the solid lines in Fig. 2 are linear in $`1/\mathrm{\Delta }t`$ and give upper and lower limits to $`\stackrel{~}{T}\tau _c`$. We estimate that $`\stackrel{~}{T}1.0`$ based on the high frame rate limit. We can now estimate the collision time to be $`0.02\mathrm{s}<\tau _c<0.04\mathrm{s}`$ for $`f=100`$ Hz and $`c=0.42`$, which corresponds to one collision every $`24`$ oscillations of the vibrator. On the other hand, for large frame rates (small $`\mathrm{\Delta }t`$), $`\stackrel{~}{T}(\mathrm{\Delta }t)`$ approaches a constant. This occurs when $`\mathrm{\Delta }t`$ is much smaller than $`\tau _c`$; in this limit $`T(\mathrm{\Delta }t)T`$, so the granular temperature $`T`$ (for displacements along one axis) can be defined as

$$T=\underset{\mathrm{\Delta }t0}{lim}T(\mathrm{\Delta }t)=<v^2>.$$
(4)   
The approach of the measured granular temperature to a constant cannot be fitted by a simple exponential or power law. For an excited granular material, this shape could be influenced by correlations in velocity between neighboring particles and correlations between the local density and particle velocity. Starting from the fact that the sampling time is roughly $`10\%`$ of the mean collision time at an intermediate coverage, we estimate that the true granular temperature $`T`$ might be up to $`10\%`$ higher than the measured value.

All values of $`T`$ and $`\stackrel{~}{T}`$ presented in this paper are the measured temperatures obtained at the highest frame rate (usually $`250`$ frames/s), where the limit $`\mathrm{\Delta }t0`$ is justified. This limit is approached in a very similar way for different accelerations (see Fig. 2). This allows us to determine the acceleration dependence of $`\stackrel{~}{T}`$. Note that the collision time can be roughly independent of acceleration since the height of the average bounce increases with increasing peak plate acceleration. For the fixed number of particles in our system this leads to an increase in the mean free path with increasing acceleration. When the number of particles is changed, the collision time also changes. As shown in Fig. 3, the limit of Eqn. (4) is approached fastest at the lowest energy of $`c=0.14`$ and slower at $`c=0.42`$ and at $`c=0.98`$. However, we are close enough to the limit at all coverages (since the mean time between collisions is always significantly longer than the time between frames) to observe qualitatively how $`\stackrel{~}{T}`$ changes with $`c`$.

We can therefore determine the dependence of the temperature on plate acceleration $`a`$ and coverage $`c`$, shown in Fig. 4. The temperature increases with $`a`$ for all coverages. As a function of $`c`$, $`T`$ increases at low coverage, exhibits a maximum around $`c=0.30`$ and decreases with $`c`$ at high coverage. This trend reflects changes in the true granular temperature as shown in Fig. 3. No measurable change in the temperature dependence occurs around $`c=1`$. A fractional coverage above unity is meant to indicate that more than one close packed layer of beads is used. The highest granular temperature indicates that the average potential energy of the particles corresponds to a mean height above the plate comparable to one particle diameter $`d`$. All experimental results are therefore limited to the regime of particle energies smaller than or comparable to the only characteristic energy of a granular material, the potential energy of raising one particle by one diameter $`mgd`$. When scaled by the peak plate velocity $`v_p`$ as in Fig. 5, the granular temperature becomes independent of acceleration to within $`\pm 10\%`$ for most data points. Remarkably, the dependence on coverage follows approximately the same behavior at all accelerations.

The scaling of velocities by $`v_p`$ implies that $`T1/f^2`$. Fig. 6 shows, on a log-log plot, that $`T`$ does indeed decrease approximately $`1/f^2`$ for the two accelerations and the three coverages shown. The plot covers two orders in magnitude of the granular temperature, ranging from conditions where beads rarely hit the container lid to conditions where frequent collisions with the lid occur. We conclude that the scaling of the bead velocity by the plate velocity is very robust and is not significantly influenced by additional contacts with the container lid.

### B Velocity Distributions

The velocity distribution along one axis, obtained from particle tracks of white beads, is shown in Fig. 7. Fits to a Gaussian distribution $`F_g(v)=A_g\mathrm{exp}(|v/\sqrt{2}v_c|^2)`$ are shown as dashed lines, and fits to the prediction of Ref. 

$$F_2(v)=A_2\mathrm{exp}(|v/1.164v_c|^{1.5})$$
(5)   
are shown as solid lines. The characteristic velocitiy $`v_c`$ is defined as the square root of the variance, so that for both $`F_g`$ and $`F_2`$:

$$v_c=\sqrt{<v^2>}=\sqrt{T}.$$
(6)   
The data points for the fits are weighted equally on a linear scale in Fig. 7(a,b) and equally on a logarithmic scale in Fig. 7(c,d). The characteristic velocity $`v_c`$ in (a,b) is $`v_c=6.19`$ cm/s for $`F_g`$ and $`v_c=6.99`$ cm/s for $`F_2`$. The Gaussian fit underestimates the probability of both low and high ($`v>3v_c`$) velocities, while the fit to $`F_2`$ describes the probability distribution quite well over three orders of magnitude in probability. The increased weight of the high velocity experimental data in (c,d) leads to $`v_c=7.66`$ cm/s for $`F_g`$ and $`v_c=6.81`$ cm/s for $`F_2`$. The Gaussian fit again underestimates high and low velocity probabilities, while the fit to $`F_2`$ proves to be insensitive to the weighting of data points, indicating a robust fit. The fit to $`F_2`$ is also insensitive to the choice of the fitted range of velocities, while the characteristic velocity decreases with decreasing fitting range for $`F_g`$. We conclude that $`F_2`$ provides a better fit than $`F_g`$.

The velocity distribution for a large range of accelerations $`3\mathrm{g}a8\mathrm{g}`$ can be described accurately by $`F_2`$ as shown in Fig. 8(a). The data and the best fitting lines are shifted vertically as needed for clarity; this amounts to multiplication by a constant on a log-linear plot. The probabilities are plotted against $`|\stackrel{~}{v}|^{1.5}`$, where $`\stackrel{~}{v}=v/v_p`$. On this log-linear scale, $`F_2`$ is a straight line, in good agreement with the experimental data for a large range of accelerations. The dependence of the velocity distribution on coverage is more complex. The coverage dependence is shown in Fig. 8(b). At relatively high coverage, above $`c=0.28`$, the experimental data can be described well by $`F_2`$ (solid lines). However, at lower coverage (for approximately the same range of $`c`$ where $`T`$ increases with $`c`$) the probability of high velocities is underestimated by $`F_2`$ (dashed line).

While the granular temperature is proportional to $`1/f^2`$ to a good approximation, measureable differences are apparent in the distribution of non-dimensional velocities $`\stackrel{~}{v}`$, as a function of frequency. Fig. 9 compares the distributions for $`40`$ Hz and $`140`$ Hz. The distribution falls more slowly with velocity as $`f`$ is increased. The cause of this behavior is probably the smaller (unscaled) velocities at higher frequency, which decreases the collision rate (i.e. the rate of energy loss through inelastic collisions), and the higher frequency of vibration (i.e. roughly the rate of energy input). For the highest unscaled particles velocities and lowest excitation frequency ($`a=5`$ g and $`f=40`$ Hz) deviations from fits to $`F_2`$ become observable.

At low accelerations $`a=2`$ g we observe that the velocity distribution has exponential tails and an approximately Gaussian central component. The crossover from a Gaussian distribution to an exponential distribution is shown in Fig. 10. These results are similar to the velocity distributions around $`a=1`$ g obtained by Olafsen and Urbach , and possibly related to clustering effects observed at low accelerations. In addition, the system is nearly two dimensional since most beads bounce lower than one particle diameter..

It is possible to investigate the case of a freely cooling granular medium to some extent by looking at a time averaged velocity distribution. We shut off the vibrator abruptly, and simultaneously trigger the camera at $`250`$ frames/s. We then extract particle tracks from the image sequence in the same way as we did for image sequences of continuously excited granular media. In order to analyse the functional form of the velocity distribution, we need to accumulate velocities over $`150`$ frames, i.e. $`0.6`$ s. Fig. 11 shows such a velocity distribution for $`c=0.84`$, where the vibrator was operated at $`f=100`$ Hz and $`a=5`$ g prior to shut off. The velocity distribution is exponential, in agreement with calculations for the free cooling case . However, the measured velocity distribution represents an average over almost the entire free cooling process, since the instantaneous velocities indicate that the granular temperature decreases by more than one order of magnitude during the averaging time. We have also made measurements over $`0.1`$ s; although the statistics are not as good, the distributions still appear to be exponential.

### C Equipartition

In a final set of experiments, the white glass beads are replaced by steel beads of different size (grade 100 stainless steel 316 for the two smaller sizes, grade 100 stainless steel 302 for the largest beads). At a total coverage $`c=0.42`$, $`14\%`$ of the glass beads were replaced by steel beads (for the smallest bead size only $`5\%`$ of the glass beads were replaced). In all experiments most collisions of steel beads therefore occur with glass beads. The distributions of $`\stackrel{~}{v}`$ for the tracked particles are shown in Fig. 12. Larger beads have smaller non-dimensional characteristic velocities $`\stackrel{~}{v_c}`$ and thus a smaller granular temperature $`T`$ than smaller beads. The velocity distributions for the two smaller steel bead sizes are described well by $`F_2`$, while that for the largest beads is better described by a Gaussian (dashed line). This could indicate that the large beads effectively prevent clustering since they act as a source of momentum for the other beads. While the granular temperature of the largest steel beads is lower than the temperature of the glass beads, their energy $`mT`$ is larger, while the energy of the smaller steel bead sizes is smaller, as shown in Fig. 13.

## V Summary and conclusion

We have reported experimental studies of velocity statistics for a (fractional) layer of glass beads subjected to vertical vibration. The horizontal motion of a small subset of beads was measured using a high speed camera at a frame rate sufficiently high to measure instantaneous velocities accurately. The measurements were acquired over an interval substantially longer than the time between interparticle collisions. These capabilities allowed us to determine particle statistics of the horizontal motion in detail. We analyzed granular temperatures and velocity distributions for a large range of excitation frequencies, amplitudes, and coverages.

The variance of the particle velocity distribution (or the granular temperature $`T`$ of horizontal motion) varies approximately in proportion to the plate velocity (Fig. 5). It increases with increasing coverage $`c`$ at low $`c`$ and decreases at higher $`c>0.42`$. On the other hand, the mean energy associated with the vertical motion probably declines with increasing $`c`$ for all $`c`$ due to additional dissipation from collisions. In this interpretation the decrease in $`T`$ with $`c`$ at high $`c`$ mirrors the decrease in energy associated with vertical motion. In contrast, the smaller values of $`T`$ found at low $`c`$ likely indicate that the energy transfer from vertical to horizontal motion becomes less efficient at low $`c0.42`$.

We have shown that particles of different mass do not have the same kinetic energy or the same granular temperature when both are present simultaneously (Fig. 13), an apparent violation of equipartition. The most reasonable explanation is that all particles acquire similar vertical velocity fluctuations from the container. The more massive particles therefore obtain a larger vertical kinetic energy. This excess vertical energy is transferred to the horizontal motion, leading to a violation of equipartition.

An important result of this investigation is that in the steady state the velocity distributions deviate measurably from a Gaussian (Fig. 7), but can be described well by $`P(v)\mathrm{exp}(|v/v_c|^{1.5})`$ for broad ranges of frequencies $`f`$, accelerations $`a`$ and coverages $`c`$, in agreement with the theory of Ref . In most of our experiments the time between interparticle collisions $`\tau _c`$ is comparable to the time between contacts with the plate. However, if forcing collisions become significantly less frequent than interparticle collisions (as for free cooling, Fig. 11), then the distribution approaches an exponential. This experimental observation is consistent with the numerical results of Puglisi et al. . This quantitative correspondence with experiment indicates that the theoretical and numerical approaches of Refs.  to describing the statistical properties of granular particles are promising.

Puglisi et al.  suggest that non-Gaussian behavior and clustering are indications of essentially the same particle dynamics. If so, our experimental results indicate that clustering must occur for a very large range of excitation amplitudes and frequencies.

Acknowledgments: This work was supported in part by the National Science Foundation under Grant No. DMR-9704301. We thank Eric Weeks and John Crocker for providing particle tracking software. We thank Yuhai Tu for valuable discussions. Technical support was provided by Bruce Boyes.