R. A. Quinn, C. Cui, J. Goree^{1,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

^{1}Author to whom correspondence should be addressed

^{2}electronic address: john-goree@uiowa.edu

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 [7]. 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 [19]. 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. [15].

Here we calculate g6(r) directly from the bonds, first by selecting a center
point bond, having its midpoint at **r**0 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 = |**r**i - **r**0| 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
**r**i, as S(k) [[equivalence]] N^{-1} < [[Sigma]]i,j exp [i
**k** . (**r**i - **r**j)] >. 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 |**r**i -
**r**j | ), 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 [7]. 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
[25]. In 3-D systems [26], melting takes place at a height of about 2.85. In
addition to these criteria, other empirical criteria have been used by
colloidal researchers [13] [26].

We have calculated g, g6 and S(k) for the plasma-crystal image of Thomas *et
al.* [1]. 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 [1].

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 [25], 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
[27].

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
[7], 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 10^{3} 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.* [1], 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 [27] 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).

**References**

[1] H. Thomas *et al*., Phys. Rev. Lett. **73**, 652 (1994).

[2] J. H. Chu and Lin I, Phys. Rev. Lett. **72**, 4009 (1994).

[3] J. H. Chu and Lin I, Physica A **205**, 183 (1994).

[4] Y. Hayashi and K. Tachibana, Jpn. J. Appl. Phys. **33**, 804 (1994).

[5] A. Melzer, T. Trottenberg, and A. Piel, Phys. Lett. A **191**, 301
(1994).

[6] T. Trottenberg, A. Melzer, and A. Piel, Plasma Sources Sci. Technol. (submitted).

[7] C. A. Murray, in* Bond Orientational Order in Condensed Matter
Systems*, edited by K. J. Strandberg (Springer-Verlag, New York, 1992),
Chap. 4.

[8] F. Diedrich *et al*., Phys. Rev. Lett. **59**, 2931 (1987).

[9] S. L. Gilbert, J. J. Bollinger, and D. J. Wineland, Phys. Rev. Lett.
**60**, 2022 (1988).

[10] I. Waki *et al*., Phys. Rev. Lett. **68**, 2007 (1992).

[11] C. A. Murray and R. A. Wenk, Phys. Rev. Lett. **62**, 1643 (1989).

[12] C. A. Murray, W. O. Sprenger and R. A. Wenk, J. Phys. Condens. Matter
**2**, SA385 (1990).

[13] C. A. Murray, W. O. Sprenger and R. A. Wenk, Phys. Rev. B **42**, 688
(1990).

[14] Y. Tang *et al*., Phys. Rev. Lett. **62**, 2401 (1989).

[15] D. G. Grier and C. A. Murray, J. Chem. Phys. **100**, 9088 (1994).

[16] B. T. Draine and E. E. Salpeter, Astrophys. J. **231**, 77 (1979).

[17] J. Goree, Plasma Sources Sci. Technol. **3**, 400 (1994).

[18] O. Havnes, G. E. Morfill and C. K. Goertz, J. Geophys. Res. 89, 10999 (1984).

[19] H. Ikezi, Phys. Fluids **29**, 1764 (1986).

[20] J. M. Kosterlitz and D. J. Thouless, J. Phys. C **6**, 1181 (1973).

[21] B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. **41**, 121 (1978).

[22] D. R. Nelson and B. I. Halperin, Phys. Rev. Lett. **19**, 2457
(1979).

[23] A. P. Young, Phys. Rev. B **19**, 1855 (1979).

[24] 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.

[25] T.V. Ramakrishnan, Phys. Rev. Lett. 48, 541-545 (1982).

[26] David G. Grier and Cherry Murray, in *Ordering and Phase Transitions in
Charged Colloids*, edited by A. K. Arora (VCH Publishers, New York, 1995).

[27] C. Cui, J. Goree and R. A. Quinn, *to be submitted*.

**Figure Captions**

FIG. 1 Particle locations and nearest-neighbor bonds, from the
experimental data reported in Fig. 2 of Thomas *et al.* [1].

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.