R. A. Quinn, C. Cui, J. Goree1,2 and J. B. Pieper
Department of Physics and Astronomy, University of Iowa, Iowa City, Iowa, 52242
H. Thomas and G. E. Morfill
Max-Planck-Institut für Extraterrestrische Physik, 85740 Garching, Germany
Pair and bond-orientational correlation functions and structure factors, as used in colloidal science, are applied as quantitative indicators of the phase of the newly discovered "plasma crystals." These Coulomb-lattice structures are formed by charged microspheres levitated in glow discharge plasmas and are imaged photographically. Static structural analysis is demonstrated on the experimental data of Thomas et al. [Phys. Rev. Lett., 73, 652, (1994)] and interpreted in the context of a 2-D dislocation-unbinding melting theory and a 2-D density-wave melting theory.
PACS: 82.70.Dd, 52.25.-b, 61.16.-d
1Author to whom correspondence should be addressed
2electronic address: email@example.com
A new cross-disciplinary area of study has emerged with the recent discovery of Coulomb lattice structures formed in dusty plasmas [1-6]. These macroscopic lattices, made of charged microspheres levitated in a glow discharge plasma, have been termed "plasma crystals." They serve as a physical model system to simulate atomic gases, liquids and solids. Model systems are valuable because experimental identification of actual atomic structure involves several difficulties . Among other physical model systems, ion crystals [8-10] have simple Coulomb interaction potentials, but rotate too rapidly for still images to be obtained, while colloidal crystals [7,11-15] are easily imaged, but have an equilibration time of months. The plasma crystal model system is attractive because it is easy to image directly, it reaches equilibrium in a few seconds, and it allows direct observation of lattice structure and dynamics.
The experimental identification of the phase of a plasma crystal has been based up to now mostly on a qualitative judgement of how much order an image appears to show, and on a comparison of histograms of the number of a particle's nearest-neighbor bonds to that of an isotropic gas [1,5,6]. Experimenters have also reported histograms of nearest-neighbor distance, the area and number of sides of Voronoi cells [1,5,6] and pair correlation functions [3,5,6]. Here we present a method of identifying liquid and solid phases quantitatively, using criteria that are founded on established melting theories and that have been used successfully already to identify the phase of 2-D colloids. We do this by analyzing the pair correlation function and two more indicators, the bond-orientational correlation function and the structure factor.
A plasma crystal is formed by charged micron-sized spheres levitated in a glow discharge plasma. The particles collect equal fluxes of electrons and ions, charging to a steady-state Q = 4[[pi]][[epsilon]]0a[[Phi]]s, where a and [[Phi]]s are the particle radius and surface potential respectively [16-18]. The individual particles attract Debye clouds of opposite charge and therefore interact via a screened Coulomb potential. When the particles' interaction potential energy exceeds their thermal kinetic energy, lattice structures can form . In most plasma crystal experiments, the particles are located in one to twenty layers above a horizontal electrode, where gravity is balanced by the sheath electric field. These layers act mainly as 2-D systems, with limited out-of-plane particle motion.
The 2-D nature of laboratory plasma crystals resembles 2-D colloids in several ways. Both are formed by charged micron-sized spheres, although colloids are in an aqueous suspension. For 2-D experiments, a colloid is confined to a single plane by an electrostatic force acting between two narrowly separated glass plates. This force, like the gravitational and sheath forces in the plasma crystal, restricts particles to motion that is mainly planar. In both systems, digitized images of particle locations show that the lattice microstructure is primarily hexagonal [1-7,11-15]. These similarities suggest applying 2-D colloidal science analysis tools to the plasma crystal.
One such tool is the pair correlation function, g(r), also known as the radial (pair) density distribution. This function represents the probability of finding two particles separated by a distance r, and it measures the translational order in the structure. For the case of a crystal at zero temperature, g(r) is a series of d functions, whose positions and heights can be determined from the particle separations in a perfect hexagonal lattice.
For this paper we calculated g(r) from experimentally measured particle locations by choosing one as a center point and counting the particles in concentric annular rings of radius r. This number is normalized by the annular ring's area and averaged using each particle as a center point. This gives g(r) in units of areal density. We normalize again by the average particle density so that the asymptotic value of g(r) is unity.
Another tool is the bond-orientational correlation function, g6(r), which is defined in terms of the nearest-neighbor bond angles of a lattice. It measures orientational order in the structure, based on the principle that all bonds in a perfect hexagonal lattice should have the same angle, modulo [[pi]]/3, with respect to an arbitrary axis. The bond angles and locations are determined from the particle locations by Delaunay triangulation (see Fig. 1). For a perfect hexagonal crystal at zero temperature, g6(r) is a constant equal to unity, while for other phases it decays with increasing r. A description of g6 in terms of a bond-orientational order parameter is given in Ref. .
Here we calculate g6(r) directly from the bonds, first by selecting a center point bond, having its midpoint at r0 and a bond angle [[theta]]0 with respect to an arbitrary axis. Next we perform a numerical average of cos [6 ([[theta]]i - [[theta]]0)] for all bonds i whose midpoints lie in an annular ring of radius r = |ri - r0| about the center point bond. This result is averaged again using each bond as a center point, yielding g6(r).
The g and g6 correlation functions can be interpreted within the context of a 2-D melting theory developed by Kosterlitz, Thouless, Halperin, Nelson, and Young, known as the KTHNY theory [20-24]. This theory predicts that 2-D crystals undergo two continuous melting transitions, in contrast to the single first-order transition of 3-D crystals. First, free dislocations (tightly bound pentagon-heptagon disclinations pairs) develop within a 2-D crystal, causing the crystal to melts into an intermediate phase, termed the hexatic, which has quasi-long-range orientational order, but short-range translational order. Second, the free dislocations unbind forming free disclinations so that the hexatic melts into an isotropic liquid phase, whose orientational and translational order are both short range. The predictions of KTHNY for translational order apply specifically to a translational correlation function gG defined in terms of reciprocal lattice vectors [20-24]. However, g is used more commonly than gG by colloidal experimenters for comparison with theory [7,11-14]. To use the KTHNY predictions, experimenters fit their g(r) and g6(r) data to exponential and power law functions.
A third tool is the static structure factor, S(k). This function is convenient when inverse methods such as Bragg scattering are used to study the structure of a solid or liquid. Plasma crystal experiments, however, yield a direct image of the structure of the particle positions, and this must be transformed into S(k). It is defined, using the experimental data for the particle positions ri, as S(k) [[equivalence]] N-1 < [[Sigma]]i,j exp [i k . (ri - rj)] >. Here the angle brackets denote an average over the direction of k, the sum is over all particle pairs i,j including i=j, and N is the number of particles. This can be re-written using the Bessel function as S(k) = N-1 [[Sigma]]i,j J0 (k |ri - rj | ), which can be easily evaluated for any value of k.
Here we list commonly-accepted criteria for identifying the phase. Those involving g and g6 are as follows . In the crystalline phase, g(r) [[proportional]] r-[[eta]]([[Tau]]) and g6(r) = const, where [[eta]]([[Tau]]) <= 1/3 is weakly temperature dependent. In the hexatic phase, g(r) [[proportional]] exp(-r/[[xi]]) and g6(r) [[proportional]] r- [[eta]] with 0 < [[eta]] <= 1/4. In the liquid phase, g(r) [[proportional]] exp(-r/[[xi]]), g6(r) [[proportional]] exp(-r/[[xi]]6) and [[xi]] = [[xi]]6. Here [[xi]] and [[xi]]6 are scale lengths for translational and orientational order respectively. In the liquid phase [[xi]] is smaller than in the hexatic. These criteria are based on the KTHNY theory. A criterion based on a density-wave theory of first-order freezing in two dimensions predicts that the first peak of the structure factor has a height of 5.25 to 5.75 at the phase transition . In 3-D systems , melting takes place at a height of about 2.85. In addition to these criteria, other empirical criteria have been used by colloidal researchers  .
We have calculated g, g6 and S(k) for the plasma-crystal image of Thomas et al. . This image was obtained in an experiment by shining a thin sheet of laser light on a horizontal lattice layer and recording the scattered light with a video camera. A single video frame was then digitized to identify the particle coordinates.
We analyzed the Thomas et al. data to find the nearest-neighbor bonds, as shown in Fig. 1. From the particle locations we calculated g(r), shown in Fig. 2. A fit to g(r), based on the d function peaks for an ideal hexagonal lattice, also appears in Fig 2. In fitting g(r), we assumed that the d functions become Gaussians of width [[sigma]], attenuated by an overall envelope that decays as exp (- r / [[xi]]). Using Eq. 4 of Ref. 15, we carried out a three-parameter non-linear fit, which for our data had a reduced chi-square of 3. These results were found: the mean interparticle spacing is [[Delta]] = 290 um, the translational scale length is [[xi]] = 2 [[Delta]], and [[sigma]] = 0.08 [[Delta]].
A plot of g6(r), computed from the Delaunay bonds, is presented in Fig. 3. To exclude incorrectly identified bonds near the image edge (see Fig. 1) from the calculation of g6(r), we only count bonds more than three mean bond spacings from the edge. Fitting g6(r) to an exponential decay yielded the scale length [[xi]]6 = (1.4 +/- 0.1) [[Delta]], while fitting to a power law decay yielded a coefficient [[eta]] = 1.0 +/- 0.05, as displayed in Fig. 3.
From these results, we find two indications of the phase that are in agreement. Both [[xi]] ~ [[xi]]6 and [[eta]] > 1/4 indicate that the structure is a fluid. The possibility of a random gas was already ruled out by Thomas et al. by comparing a histogram of their coordination data to that of a random lattice .
A third indication comes from the structure factor S(k), shown in Fig. 4. The first peak has a height of 3.2. Based on the 2-D density-wave theory , this value lies well on the liquid side of the single first-order phase transition.
An interesting issue is how the structure can be so disordered as to be a liquid when the charge on a grain is so large. Thomas et al. estimated the charge was between -9800 and -27 300, the Coulomb coupling parameter was [[Gamma]] > 20 700, and the ratio of the particle spacing to the linearized Debye length was in the range 0.6 <= [[kappa]] <= 4.8. Thus, the shielded coupling parameter was [[Gamma]] exp (-[[kappa]]) > 170, which is so large as to predict a crystalline state. The limitations of using OCP theory for plasma crystal experiments will be discussed in another paper .
Unlike the structure factor, the correlation functions are used by fitting them to a curve, and this introduces several kinds of errors into the results for [[xi]], [[xi]]6, and [[eta]]. These errors are small enough for the results and conclusions presented above, but in future applications such as analyzing a structure near its phase transition they may deserve close scrutiny, so we list them here. Particle and bond locations in the experimental images are assigned to radial bins that accumulate a finite number of counts, so that the g(r) and g6(r) data have error bars from counting statistics. These error bars are too small for the fitted curve to pass through, so for the purpose of calculating the uncertainties in the fit parameters [[xi]]6 and [[eta]] as reported above, we scaled the bins' error bars up by a constant factor so that the fit had a reduced chi-square of unity. This method of assigning an uncertainty to our results, however, does not reflect other uncertainties from the fitting. The results for [[xi]], [[xi]]6, and [[eta]] depend on which range of radii is included, and this range is chosen somewhat arbitrarily. They can also depend on the weighting that is used in the fitting; if we had not weighted data inversely by the square of the individual error bars, we would have found [[xi]]6 = 2.0 instead of 1.4 and [[eta]] = 1.0 instead of 1.05. In fitting data to curves, the curves model the data imperfectly; for example, we have found in various plasma crystal experiments that the first peak of g(r) is usually much higher in the experimental data than in the smooth curve. Finally, as for colloidal experiments and simulations , a single image contains a finite number of defects that fluctuates statistically, so that a definitive phase identification near a phase transition might require averaging over as many as 103 uncorrelated images. The single image used here seems to be sufficient to identify the phase, as all three indicators are in agreement and not very near the values where a phase transition is predicted.
As in the case of 2-D colloidal suspensions, plasma crystal experiments do not fully satisfy the assumptions of the melting theories that we have used to identify the phase. The KTHNY and 2-D density-wave theories assume identical particles, but in the experiment of Thomas et al. , there was a 3% dispersion in particle size, corresponding to a 3% variation in particle charge. Colloid experiments have comparable size dispersions. While the theories are 2-D, the experiments have some 3-D effects. An empirical measure of this is the average time [[tau]]plane that particles remain in the illuminated plane. For the Thomas et al. data, [[tau]]plane ~ 4 min, while in more recent experiments  out-of-plane motions have been reduced so that [[tau]]plane >> 100 min. (Phases have been identified in these other experiments that are more ordered and more disordered than the fluid phase reported here.) The cause of 3-D effects are believed to be slight curvature in the electric sheath that levitates the particles and local perturbations induced by particles in other planes, including aggregates of two spheres in the lowest lattice planes. Another limitation in comparing the experiments and theories is that a static analysis of an individual frame does not reveal whether a defect is stationary (as would be typical of a solid) or moving. As a future research direction, we believe that dynamical analyses can be of great usefulness in correcting these limitations. This is practical since the time scales (seconds) of plasma crystals are well-suited for temporal analysis of video images.
To conclude, we have demonstrated the use of static correlation functions and the structure factor in determining the phase of a plasma crystal. Applying these techniques represents a step in developing the plasma crystal as a viable model system for the study of atomic structure, dynamics and phase transitions.
The authors thank C.A. Murray for useful suggestions and literature. This work was supported by NASA (NAG8-292, NAGW-3126, NGT-51351) and NSF (ECS-92-15882).
 H. Thomas et al., Phys. Rev. Lett. 73, 652 (1994).
 J. H. Chu and Lin I, Phys. Rev. Lett. 72, 4009 (1994).
 J. H. Chu and Lin I, Physica A 205, 183 (1994).
 Y. Hayashi and K. Tachibana, Jpn. J. Appl. Phys. 33, 804 (1994).
 A. Melzer, T. Trottenberg, and A. Piel, Phys. Lett. A 191, 301 (1994).
 T. Trottenberg, A. Melzer, and A. Piel, Plasma Sources Sci. Technol. (submitted).
 C. A. Murray, in Bond Orientational Order in Condensed Matter Systems, edited by K. J. Strandberg (Springer-Verlag, New York, 1992), Chap. 4.
 F. Diedrich et al., Phys. Rev. Lett. 59, 2931 (1987).
 S. L. Gilbert, J. J. Bollinger, and D. J. Wineland, Phys. Rev. Lett. 60, 2022 (1988).
 I. Waki et al., Phys. Rev. Lett. 68, 2007 (1992).
 C. A. Murray and R. A. Wenk, Phys. Rev. Lett. 62, 1643 (1989).
 C. A. Murray, W. O. Sprenger and R. A. Wenk, J. Phys. Condens. Matter 2, SA385 (1990).
 C. A. Murray, W. O. Sprenger and R. A. Wenk, Phys. Rev. B 42, 688 (1990).
 Y. Tang et al., Phys. Rev. Lett. 62, 2401 (1989).
 D. G. Grier and C. A. Murray, J. Chem. Phys. 100, 9088 (1994).
 B. T. Draine and E. E. Salpeter, Astrophys. J. 231, 77 (1979).
 J. Goree, Plasma Sources Sci. Technol. 3, 400 (1994).
 O. Havnes, G. E. Morfill and C. K. Goertz, J. Geophys. Res. 89, 10999 (1984).
 H. Ikezi, Phys. Fluids 29, 1764 (1986).
 J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 (1973).
 B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121 (1978).
 D. R. Nelson and B. I. Halperin, Phys. Rev. Lett. 19, 2457 (1979).
 A. P. Young, Phys. Rev. B 19, 1855 (1979).
 D. R. Nelson, in Phase Transitions and Critical Phenomena, Vol. 7, edited by C. Domb and J. L. Leibowitz (Academic Press, London, 1983), p.1.
 T.V. Ramakrishnan, Phys. Rev. Lett. 48, 541-545 (1982).
 David G. Grier and Cherry Murray, in Ordering and Phase Transitions in Charged Colloids, edited by A. K. Arora (VCH Publishers, New York, 1995).
 C. Cui, J. Goree and R. A. Quinn, to be submitted.
FIG. 1 Particle locations and nearest-neighbor bonds, from the experimental data reported in Fig. 2 of Thomas et al. .
FIG. 2 Pair correlation function versus normalized distance. Here g(r) is normalized by the average areal density of particles, which is 1200 cm-2. Experimental and best-fit pair correlation functions are shown. The d functions for a perfect crystal are shown as peaks with correct positions and relative heights. The fit yielded the mean inter-particle spacing [[Delta]] = 290 um and a correlation length [[xi]] = 2 [[Delta]].
FIG. 3 Bond angle correlation function g6 versus normalized distance. The vertical axis is log-scale. Exponential and power law fits are shown, yielding [[xi]]6 and [[eta]] respectively.
FIG. 4 Structure factor S(k), computed from the measured particle locations.