Large-scale steady-state structure of a 2D plasma crystal

An analytical model, based on a simple physical analogy, and a non-linear analytic/numerical model of the large-scale structure (averaged over the size of a single cell) of a two-dimensional (2D) lattice have been developed. In the first model, a physical analogy between the lattice layer steady state and the stressed state of a rotating solid body was used to derive the model equations and the model is shown to be in good agreement with the results of a 2D simulation. The non-linear model is derived from the force balance between external and internal forces in the continuum limit and compares favourably with an experimental example.


74.3
on the other hand, have been obtained under the conditions where κ = a/λ D ∼ 0.5-2 and the above-mentioned analytical models need to be extended to this interaction range, a task we undertake a first step towards here.
In the next section, we propose a model based on a simple physical analogy, developing it for a 2D cylindrically symmetric plasma crystal as well as a 1D planar system. This model is compared favourably with simulation results in section 3. In section 4, a non-linear model is developed and compared with experimental results in section 5. Both of these models take into account the effectively long-range nature (even for these screened Coulomb systems) of the interaction potential.

Simple model
We develop our simple model for the macroscopic density distribution n in a plasma crystal through an analogy with the stretched state of a solid cylinder rotating around its symmetry axis with a fixed angular velocity. The equation for the cylinder is [9]: where u is the displacement vector defined at each lattice site, C s is the compressional sound speed defined in terms of the solid body elasticity modulus, and rot is the angular velocity of the rotation. The equation describing plasma crystal steady state can be written down immediately by analogy: where we have changed the sign on the right-hand side since the plasma crystal is compressed by the parabolic confinement potential. Now, C s is the dust-lattice sound speed of the plasma crystal and is the 'harmonic frequency' of the horizontal confinement (i.e. a dust grain having mass m has a potential energy m 2 r 2 /2 in the confining field). Note that the right-hand side 2 r is just the first non-zero term in the expansion of a general confinement force F ext , assuming azimuthal symmetry.
Equation (2) can be formally derived in the same way as equation (1), i.e. by introducing the stretch modulus. Since this procedure is already well known, we use the simple physical analogy of a stretched rotating solid to illustrate the principles of the full derivation. Note that equations (1) and (2) are the linear elastic deformation relations while we have a non-linear system. Fortunately, however, near the crystal centre the deformations of the particle positions from the 'ideal' crystal case are small and the linear theory may, in a first approximation, be applied.
For a plasma crystal, C s is given by the well-known relation [10] where Q and m are the charge and mass per grain and κ is defined above. The function F(κ) is computed as a sum over individual grain positions for an ideal 2D lattice, and can be usefully fitted by the following expression: Equation (4) fits with better than 0.5% accuracy for κ < 5, the range considered here (see figure 1).
The open circles give F(κ) from the numerical summation, based on the results in [10], and the curve is the fit. The next step is to obtain an expression for the density distribution in a lattice in terms of the displacements u. Specifically, u represents the steady-state displacements of the particles in the plasma crystal from their ideal locations in a uniform 2D crystal. When the plasma crystal is treated as a continuous medium (see [11,12] and section 4 below), the grain number density can be written: where n = n(r ) is the spatial density distribution and n 0 is the density at the crystal centre. Equation (5) becomes an exact equality in the continuum limit, while here the approximate equality is due to coarse graining. Note also that n = 2/(3 1/2 a 2 ) for the perfect 2D triangular lattice. In this paper, we consider density changes on scales much larger than a.

Solution for the cylindrically symmetric case
Next we want to solve equations (2) and (5) for the cylindrically symmetric case. To do so, we first need to construct proper boundary conditions. Here, the analogy between the rotating solid body and the plasma crystal breaks down due to the different types of interparticle potentials: attractive-repulsive van der Waals type of interaction for the first but only repulsive for the second 1 . Hence, the boundary conditions can be different. For example, in the case of rotationinduced deformation, the surface of the solid body should be free from stresses [9]. Usually, it is assumed that the normal component of the stress tensor should be equal to zero on the surface and that the displacement should be finite throughout. Undoubtedly, the second condition should be fulfilled in the present case as well. However, the first does not apply for the present case, so we take as the second boundary condition n = n 0 at the crystal centre.

74.5
Applying these boundary conditions (finite density at r = 0 and n 0 /n = 1 at r = 0) to equation (2), the solution is straightforward and takes the form According to equation (6), observation of the structure near the crystal centre immediately gives information on the confinement parameter, if the sound speed was found independently. For example, C s can be found by analysing the naturally occurring phonon spectrum [13], or by shock or soliton propagation [11].

Solution for the 1D case
Under the action of non-symmetric external forces, and/or perturbing internal forces, such as chains of defects, the crystal structure may be approximately 1D planar (see [11] and below). Here, the density in the crystal varies primarily along the x direction. The solution is straightforward and, as in equation (6), one obtains

Simulation results
To check our model, a molecular dynamics simulation was performed. The simulation consisted of a 2D crystal layer containing 721 particles, each having a charge of Q = 16 000e and interacting via a Yukawa-type (screened Coulomb) interparticle potential. The grain charge and screening length were kept constant during the simulation. The particles were confined by a cylindrically symmetric parabolic potential well and cooled by interaction with an Epsteintype drag force. The particles were placed randomly into the simulation box and the system was allowed to gradually approach equilibrium. The drag force was only sufficient to suppress numerical noise in the system and remove the initial kinetic energy. Finally, in equilibrium, the particles had a kinetic temperature of about 0.0001 eV. See [11] for further details about the simulation. The equilibrium distribution is shown in figure 2, along with the Voronoi map for the lattice. The Voronoi map defines the set of points around each lattice site which are closer to that lattice cite than the surrounding nearest neighbours. This is also known as the Wigner-Seitz cell. The inverse of the area contained in a given Wigner-Seitz cell is a good approximation of the local area density in the lattice. Here, we compute the local interparticle spacing a = (3a The next step is to fit equation (6) to the simulation data and compare the resulting coefficients with the known values in the simulation. Rather than using it directly, we integrate equation (6) to obtain where N (r ) is the radial distribution function and k ≡ 2 /(2C 2 s ) is the fit coefficient. The function along with the fit is shown in figure 3. Note that the fit is good in the vicinity of the crystal centre and off towards the edge, as expected for regions where deformations and forces 74.6 are more non-linear, and we obtain k fit = 0.31 cm −2 . To judge the quality of our model, we compute a predicted k sim from the known parameters of the simulation. Using the charge, mass, screening length and measured centre interparticle spacing a 0 in equations (3) and (4), we obtain C s = 2.57 cm s −1 and, with = 2 s −1 , k sim = 0.30 cm −2 . Therefore, k sim ≈ k fit within 3%.
As a further test of the simulation and of our model, a subsonic compression wave was launched through the simulation volume. The disturbance was tracked by computing the average density in the crystal as a function of time and space. The resulting pulse trajectory is shown in figure 4. Fitting the pulse-peak trajectory yields an estimate for the sound speed in the lattice, C fit = 2.50 cm s −1 . Note that this value is very close to that estimated above from the simulation parameters, as expected. Using this value in conjunction with our fitted value for k above yields an estimate for the confinement parameter, = 1.97 s −1 , in good agreement with the simulation value.

Non-linear model
While reasonably good at describing the equilibrium density distribution near the plasma crystal centre, the simple model proposed above fails in two ways: (1) the sound speed of an inhomogeneous crystal is position dependent (C s ≡ C s (κ), κ ≡ κ(x)), and (2) the grain displacement is actually large for large x. Here, we consider a 1D planar non-linear equation to describe the equilibrium density distribution.  Trajectory of a subsonic compression pulse propagated in the simulation crystal shown in figure 2. The dashed line is a fit for the sound speed, C fit = 2.5 cm s −1 . The spatial scale in the plot is normalized to λ D .
First, we rederive equation (5) for the 1D case, assuming uniform density in the direction perpendicular to x. Starting with an infinite homogeneous crystal, we cut a slice of thickness S x from this crystal, perpendicular to the x axis. In the absence of external forces, the crystal will expand in x until it is destroyed. Suppose, however, that we stop the expansion by applying a symmetric external force which vanishes at the line of symmetry (x = 0) but grows with distance. This is the simplest interesting case and it is often realized in experiments. When the confining force is large enough, the slice can remain crystalline, but the particles will change their positions with respect to the initial state. The crystal is deformed.
To describe the deformation more precisely, consider a narrow element dx of the original (infinite) lattice at position x. The number of grains inside this element is dN = n 0 S y dx. Here, n 0 is the grain number density per unit area before deformation, and S y is the (constant) in-plane width perpendicular to x. After deformation, this element occupies a new position X = x + u, where u is the displacement. Assume that the number of grains in this element remains the same (dN = nS y d X ), so that conservation of the grain number takes place. This is equivalent to assuming that no defects are generated through the deformation (the deformation is elastic). We then have equivalent to equation (5). Now, if the forces acting on every element of a crystal layer are in balance, then the layer is in equilibrium. In the case of our deformation, the external forces should compensate a gradient of internal stresses: then the force balance can be written as where F ext is the external force acting on a given test grain at position x, and n −1 d p/dx is the gradient of internal stresses (or 'pressure', by the terminology used in [4]) per grain. Assuming a parabolic confining well in X (the equilibrium coordinates), F ext = −m 2 X and using the relation equation (10) can be rewritten as Equations (3), (4), (9), (12) are the equations of the non-linear model describing the large-scale plasma crystal structure. Assuming a constant screening length λ D , then we additionally have d X/dx = n 0 /n = a/a 0 = κ/κ 0 . Note that equation (11) is valid if the internal stresses depend on the coordinates only indirectly via the crystal density coordinate dependence. It is not valid for the case of a crystal with defects. This has a very simple explanation: if a defect appears, it causes additional stresses in the crystal. The distribution of these stresses depends mainly on the spatial distribution of the defects, not just on the local crystal density. Here, we concentrate on weakly defected crystals, leaving the topic of particle distribution in defected crystals to another paper.
It is not too difficult to find a formally exact solution to the above set of non-linear equations. First, combining equations (9) and (12), we obtain We have, on integrating, This relation gives us the dependence X = X (κ) needed for further integration, yielding Considering equation (4), the dependence of the sound speed on κ is too complicated in general to hope for a simple analytical result in integrating equations (14) and (15) above. Nevertheless, it is interesting to construct simple analytical, albeit approximate, results. To solve this problem, note that the plot in figure 1 is essentially linear for κ < 1, implying that F(κ) has a power law dependence on κ in this regime. Referring again to equation (4), the dependence is approximately F = F 0 κ −2 , with F 0 = constant (see also [10]). Using this approximation, we obtain the following analytic solution: where L = C s (κ 0 )/ . Note that this solution agrees well with equation (7) in the vicinity of x = 0. Equation (16) can be used to fit the observed crystal structure, to find the characteristic length L of the inhomogeneity. Then, the confinement parameter can be computed after obtaining the crystal sound speed in a different experiment.

Experimental results
A particular example of an experimentally observed large-scale crystalline structure is shown in figure 5. It is a single-layer plasma crystal consisting of 6.5 µm grains embedded in the sheath of a low-pressure (1.8 Pa) rf plasma. A horizontal wire lies to the right of the area shown in figure 5. It was used to launch shock waves and solitons into the crystal. These results, as well as a complete description of the experiment, are presented in [11]. Stresses due to the wire induce defect formation in the equilibrium lattice and we suppose that these defects, along with the bulk deformation of the lattice due to the presence of the wire, are responsible for the 1D nature of the lattice deformation. Note that the lattice deformation becomes less 1D near the edges and we select only the central half of the image for further processing.
In order to obtain information about the effective confinement parameter, we start by computing (a) (equation (16)) for the data shown in figure 5. First, we compute a 0 and a(x) following the procedure described above in section 3. Next, we make a two-parameter fit using equation (16) (x → x − x 0 ) to the experimentally determined (a). This is shown in 74.10  [11], we know that the sound speed in the lattice was C s ≡ C s (κ 0 ) = 2.3 ± 0.3 cm s −1 . The corresponding estimate for the confinement parameter is then = 0.7 ± 0.4 s −1 , where most of the error comes from the sound velocity measurement.
Assuming that , Q and λ D are constant in the vicinity of the crystal centre, it is possible to solve equations (14) and (15) numerically for different values of x 0 , L and κ 0 . The hope would be to obtain an additional fit to κ 0 . However, it turns out that the fit is relatively insensitive to the value of κ 0 , at least for the experiment tested here, in the range 0.5 < κ 0 < 2.

Conclusion
Two models, one based on a simple analogy between a plasma crystal and a rotating rigid body that yields an analytic solution and a non-linear model based on a continuum description of the force balance, have been presented. The simple model describes well the equilibrium particle distribution for a simulated cylindrically symmetric 2D plasma crystal. The non-linear model, which yields an approximate analytic solution for κ less than about 2, was solved for the 1D case and used to fit an experimental distribution. In conjunction with data for the sound speed obtained in a separate wave propagation experiment, an estimate for the horizontal confinement parameter was obtained. Both models improve on previous work by taking into account the effect of having screening lengths of the order of the interparticle spacing, in line with previous experimental work in the field.