Continuum theory for two-dimensional complex plasma clusters

In a continuum approximation for the particle density ρ(r), the total potential over a particle located at r is given by

$$\begin{equation} \phi(\mathbf{r})=\int_{r'\leqslant R_m} \rho(\mathbf{r'})\frac{q^2\,{\rm e} ^{-\kappa|\mathbf{r}-\mathbf{r'}|}}{|\mathbf{r}-\mathbf{r'}|}\,{\rm d}^2r', \end{equation} \tag{ 4 }$$

where R_m is the system's radius. When the range of the potential, 1/κ, is small compared to the system size (i.e. when κR_m ≫ 1), a small region around r gives the main contribution to the integral of equation (4). In this case, the effective region of integration is approximately a circle with radius of the order of 1/κ and centered at the position vector r. For such a circular symmetry in the region of integration, we can expand ρ(r') around r (i.e. $\rho(\mathbf{r'})=\rho(\mathbf{r})+(\mathbf{r'}-\mathbf{r})\cdot\boldsymbol{\nabla}\rho(\mathbf{r})+|\mathbf{r'}-\mathbf{r}|^2\nabla^2\rho(\mathbf{r})/4+\cdots$) and one can see that the second term vanishes in the integral of ρ(r') e^{−κ|r−r'|}/|r − r'|. Therefore, the use of just the first term, i.e. of the local density ρ(r), may give a good approximation for ϕ(r). This is called the LDA.

On the other hand, when R_m is equal to a few units of 1/κ or less, FSE must be considered. These effects were taken into account only in the 3D case by Henning et al [9] and we will use a similar FSE consideration in 2D. Even with a constant density, the integral of equation (4) cannot be given in terms of elementary functions. In this case, we consider a position independent FSE by simply integrating over the region |r' − r| ⩽ R_m, i.e. a disc of radius R_m centered at r, which gives

$$\begin{equation} \fl \phi(\mathbf{r})=\phi(\rho(\mathbf{r}))=\rho(\mathbf{r})\int_{|\mathbf{r'}-\mathbf{r}|\leqslant R_m}\frac{q^2\,{\rm e}^{-\kappa|\mathbf{r}-\mathbf{r'}|}}{|\mathbf{r}-\mathbf{r'}|}\,{\rm d}^2r'=\frac{2\pi q^2\rho(\mathbf{r})}{\kappa}(1-{\rm e}^{-\kappa R_m}). \end{equation} \tag{ 5 }$$

As the system has locally a predominant triangular lattice symmetry, the density is related to the nearest particles' distance x by $\rho =2/(\sqrt {3}x^2)$. In the limit of strong screening, i.e. a big value of κ and therefore a short range of the pair potential, each particle effectively interacts only with its nearest-neighbors. In this case, the total potential ϕ(r) is approximately equal to 6 (the number of nearest-neighbors of the triangular lattice) times V_pair(x(r)).

Now we make the consideration that each particle in the cluster always interacts with only six effective particles at a distance x (see figure 1(a)) through an effective potential given by V_eff(x) = ϕ(x)/6, even when the interaction has a long-range character. In fact, each effective particle is in a nearest-neighbor position but represents the interaction with all the particles in one sixth of the system as is shown in figure 1(a). In the short-range case, the effective potential becomes the pair potential and the effective particles become the nearest-neighbor particles. We considered that the total potential over a particle can be given explicitly in terms of x, i.e. ϕ(r) = ϕ(ρ(r)) = ϕ(x(r)), which is always true for the LDA without a position dependent FSE.

Zoom In Zoom Out Reset image size

Equation (5) is a continuum approximation and therefore it is accurate only when the interaction range 1/κ is much bigger than the particles' spacing x (κx ≪ 1). In fact, for Yukawa systems, this approximation fails when κx ≫ 1, or even κx ∼ 1 [18], and one should consider correlation (discretization) effects. In this case, each particle interacts, effectively, only with few particles in its neighborhood. As those particles have a triangular lattice-like arrangement, we approximate their interaction energy by that given by a perfect triangular lattice. A good approximation for this energy was obtained in [18] by summing up a set of well-defined rings. To include the FSE, we must perform the sum up to the ring of radius ≈R_m, which yields

$$\begin{equation} \phi(x)=\frac{6q^2\,{\rm e}^{-\kappa x}}{x\left(1-{\rm e}^{-\frac{3\sqrt{3}}{2\pi}\kappa x}\right)}\left(1-{\rm e}^{-\kappa R_m}\right). \end{equation} \tag{ 6 }$$

Notice that in the regime of κx ≪ 1 (or ρ/κ² ≫ 1) equation (5) is recovered. In the rest of this article, we will use equation (5) and concentrate on 10⁻⁴ ⩽ α ⩽ 10¹. Nevertheless, the already good results obtained for average and great values of α (∝κ³) [11, 12] can still be further improved by using equation (6).

In order to compensate the external force, $\mathbf {F}_{\mathrm { ext}}=-V_0r\hat {\mathbf {r}}$, one must consider a deformation of the original hexagon formed by the effective particles. To do so, we choose the dislocation shown in figure 1(b). In this case, the distance between the central particle and the three effective particles in the positive (negative) direction of $\hat {\mathbf {r}}$ was increased (decreased) by Δx/2. The angle between two neighbor particles remains π/3 and the total potential remains the same up to the first order in Δx.

The magnitude of the force due to an effective particle is $F(x)=-\frac {\partial V_{\mathrm { eff}}}{\partial x }=-\frac {1}{6}\frac {\partial \phi }{\partial x}$. Therefore, from figure 1(b) and from the equilibrium of forces, we have

$$\begin{equation} \left(1+2\cos (\pi/3)\right)\left[F(x-\Delta x/2)-F(x+\Delta x/2)\right]=V_0r \end{equation} \tag{ 7 }$$

which, for small Δx, becomes

$$\begin{equation} -2\frac{\partial F}{\partial x}\Delta x=V_0 r. \end{equation} \tag{ 8 }$$

For a smooth dependence of x with r we associate Δx/x to the derivative of x with respect to r. This is based on the fact that the variation of the nearest-neighbor distance from the bottom rhombus to the top one, i.e. the difference between their sides, in figure 1(b), is Δx and the radial distance between their centers is x. By substituting Δx = x(dx/dr) in equation (8), we obtain the following differential equation:

$$\begin{equation} \frac{x}{3}\frac{\partial^2 \phi}{\partial x^2}\frac{{\rm d}x}{{\rm d}r}=V_0r. \end{equation} \tag{ 9 }$$

This last equation written in terms of the density becomes

$$\begin{equation} \left(\frac{\partial\phi}{\partial\rho}+\frac{2}{3}\rho\frac{\partial^2\phi}{\partial\rho^2}\right)\boldsymbol{\nabla}\rho+\boldsymbol{\nabla}V_{\rm ext}=0, \end{equation} \tag{ 10 }$$

where V_ext(r) represents a general external potential.

The obtained differential equation has a general form different from other well-known approaches. In a variational approach [11], the minimization of the total potential energy functional

$$\begin{equation} U[\rho]=\frac{1}{2}\int\limits_{r\leqslant R_m} \rho(\mathbf{r})\phi(\rho(\mathbf{r}), \mathbf{r})\,{\rm d}^2r+\int\limits_{r\leqslant R_m} \rho(\mathbf{r})V_{\rm ext}(\mathbf{r})\,{\rm d}^2r \end{equation} \tag{ 11 }$$

under the constraint $\int \rho (\mathbf {r}) \,{\mathrm { d}}^2r=N$, gives

$$\begin{equation} \frac{1}{2}\boldsymbol{\nabla}\left(\phi+\rho\frac{\partial\phi}{\partial\rho}\right)+\boldsymbol{\nabla}V_{\rm ext}=0. \end{equation} \tag{ 12 }$$

Whilst in a hydrostatic approach [12], the Euler equation

$$\begin{equation} \frac{\boldsymbol{\nabla}P}{\rho}+\boldsymbol{\nabla}V_{\rm ext}=0, \end{equation} \tag{ 13 }$$

where P = −∂(ϕ/2)/∂(1/ρ) = (ρ²/2)∂ϕ/∂ρ is the pressure, gives

$$\begin{equation} \frac{1}{\rho}\boldsymbol{\nabla}\left(\frac{\rho^2}{2}\frac{\partial\phi}{\partial\rho}\right)+\boldsymbol{\nabla}V_{\rm ext}=0. \end{equation} \tag{ 14 }$$

In spite of the differences between equations (10), (12) and (14), by using an approximation for ϕ of the form ϕ(r) = Bρ(r), where B is a constant, they provide the same result for the density, i.e. ρ(r) = ρ₀ − V_ext(r)/B where ρ₀ is a constant of integration. The differences between these equations appear when ϕ has other dependences with ρ and r. These differences will be investigated in future work. Now we use the result ρ(r) = ρ₀ − V_ext(r)/B together with the approximation of equation (5) and the external potential V_ext(r) = V₀r²/2 to obtain the following equation for the density:

$$\begin{equation} \frac{\rho(r)}{\kappa^2}=\frac{\rho(0)}{\kappa^2}-\frac{\kappa^2r^2}{4\pi\alpha\left(1-{\rm e}^{-\kappa R_m}\right)}. \end{equation} \tag{ 15 }$$

The values of ρ(0)/κ² and κR_m are obtained from the boundary and normalization conditions. As remarked in section 2, they must depend only on the parameters α and N.

A continuum approximation for the density must satisfy the normalization condition $\int \rho (\mathbf {r})\,{\mathrm { d}}^2r=N$, which in our case is written as

$$\begin{equation} 2\pi\int_0^{R_m}\rho(r)r\,{\rm d}r=N. \end{equation} \tag{ 16 }$$

This equation can be used to eliminate one unknown term of equation (15), while the second one can be eliminated from a boundary condition.

The most simple boundary condition is to say that ρ(R_m) = 0. However it does not indeed happen, although becomes a good approximation in the large N limit where, as observed in experiments and simulations, lim_N→∞ρ(R_m)/ρ(0) = 0 ∀α. The use of this condition and equation (16) in (15) gives

$$\begin{equation} \frac{\rho(\kappa r)}{\kappa^2}=\frac{\kappa^2 R_m^2-\kappa^2r^2}{4\pi\alpha\left(1-{\rm e}^{-\kappa R_m}\right)}, \end{equation} \tag{ 17 }$$

where

$$\begin{equation} \frac{(\kappa R_m)^4}{\alpha(1-{\rm e}^{-\kappa R_m})}=8N. \end{equation} \tag{ 18 }$$

Except by the factor (1 − e^−κR_m) descendant from FSE, this result was obtained by Totsuji et al [11] by using a variational approach. There, the result ρ(R_m) = 0 comes out from the minimization of the total energy obtained from the LDA without FSE and correlation effects. When the latter effects were considered [11], a non-zero density at the edge came out naturally, satisfying lim_N→∞ρ(R_m)/ρ(0) = 0.

Notice that, from the effective particles approach of figure 1(b), a straightforward boundary condition appears. When the central particle of figure 1(b) is at the boundary of the cluster, the three effective particles at the top do not exist. Therefore, the equilibrium of forces considered in equation (7) is now written as (1 + 2 cos(π/3))F(x_m) = V₀R_m, where x_m is the nearest-neighbor distance at the boundary, or equivalently

$$\begin{equation} -\frac{1}{3}\frac{\partial \phi}{\partial x}\Big|_{x_m}=V_0R_m. \end{equation} \tag{ 19 }$$

Despite the consideration of just three effective particles in this case, the effective potential continues to be given by ϕ(x)/6. This is in fact a correction for our choice of the term (1 − e^−κR_m) for the FSE. This choice overestimates the potential at the edge by a factor that approximates the number 2 as κR_m increases.

By placing the potential ϕ of equation (5) into equation (19) we obtain an expression relating x_m and R_m

$$\begin{equation} -\frac{1}{3}\left[\frac{-8\pi q^2}{\sqrt{3}\kappa x_m^3}\right]\left(1-{\rm e}^{-\kappa R_m}\right)=V_0R_m. \end{equation} \tag{ 20 }$$

An important and non-trivial consideration must be made here: for a particle at the edge, just a fraction of its Wigner–Seitz cell is inside the circle of radius R_m (see figure 2). This fraction is approximately one half for R_m/x_m > 1. Due to this fact, the density (given by the inverse of Wigner–Seitz cell area) at the edge is related to x_m by

$$\begin{equation} \rho(R_m)=\frac{4}{\sqrt{3}x_m^2}. \end{equation} \tag{ 21 }$$

This correction in the density should be made always that the superior limit of integration of the normalization condition shown in equation (16) is R_m. In this case, the total area occupied by all Wigner–Seitz cells must be equal to πR²_m but it is overestimated if we always have $A_{\mathrm { WS}}=1/\rho =\sqrt {3}x^2/2$ (see figure 2). This problem has not been considered previously in the literature but it becomes important in small systems. Equation (21) is not exact but gives a good approximation to correct this problem.

Zoom In Zoom Out Reset image size

Using equations (21) and (20), we obtain

$$\begin{equation} \frac{\rho(R_m)}{\kappa^2}=\sqrt{3}\left(\frac{\kappa R_m}{\pi \alpha \left(1-{\rm e}^{-\kappa R_m}\right)}\right)^{2/3}. \end{equation} \tag{ 22 }$$

The above equation is one of the main theoretical results obtained in this paper. It will be determinant for the good accuracy of the theory at small number of particles. By using equation (22) in equation (15), we can find ρ(0) in terms of R_m. Finally, the density at any point is then given by

$$\begin{equation} \frac{\rho(\kappa r)}{\kappa^2}=\sqrt{3}\left(\frac{\kappa R_m}{\pi \alpha \left(1-{\rm e}^{-\kappa R_m}\right)}\right)^{2/3}+\frac{\kappa^2 (R_m^2-r^2)}{4\pi\alpha\left(1-{\rm e}^{-\kappa R_m}\right)}, \end{equation} \tag{ 23 }$$

where κR_m can be found by using equation (16), which gives

$$\begin{equation} \sqrt{3}\left(\frac{\sqrt{\pi}(\kappa R_m)^4}{\alpha \left(1-{\rm e}^{-\kappa R_m}\right)}\right)^{2/3}+\frac{(\kappa R_m)^4}{8\alpha(1-{\rm e}^{-\kappa R_m})}=N. \end{equation} \tag{ 24 }$$

For a given number of particles, equations (23) and (24) inform that the distances in the system scale with R_m, i.e.

$$\begin{equation} \rho(r)=R_m^{-2}f(N,r/R_m), \end{equation} \tag{ 25 }$$

where f(N,r/R_m) is independent on the potentials' parameters. This would imply that the values of κ and α cannot be obtained separately from the final configuration of an experiment, but only from a relation between them. Another experiment with a different N must be done in order to get a new relation.

A density profile obeying equation (25) can be obtained when the potential ϕ is, up to a constant, approximated by ϕ(r) = Bρ(r). By using this approximation, the method of the previous two subsections (i.e. any of equations (10), (12) and (14) together with $\rho (R_m)=4/(\sqrt {3}x_m^2)$ and equation (19)) gives

$$\begin{equation} R_m^2\rho(r)=f(N,r/R_m)=\sqrt{3}c^2+\frac{c^3}{4}\left(1-\frac{r^2}{R_m^2}\right), \end{equation} \tag{ 26 }$$

where c = c(N) is obtained from

$$\begin{equation} c^3+8\sqrt{3}c^2=8N/\pi \end{equation} \tag{ 27 }$$

and R_m is related to B by

$$\begin{equation} 2V_0R_m^4=Bc^3. \end{equation} \tag{ 28 }$$

The factor B can be a function of N, R_m and of the potentials' parameters; and may be obtained, for example, by: (a) an integral of the pairwise potential $B=\int V_p(r)\,{\mathrm { d}}^2r$, as it was done in equation (5), or (b) by a derivative of the Madelung energy (equation (6)), B = ∂ϕ/∂ρ, evaluated at the mean density $\skew3\bar {\rho }=N/\pi R_m^2$. Indeed, the approximation ϕ = Bρ is very general and can be used in many systems beyond Yukawa [16, 17], and so does equations (26) and (27).

The scale invariance of the theory is broken when position dependent FSE (small κR_m) or correlation effects (small ρ/κ²) are taken into account. Using the latter, Totsuji et al [19] were able to obtain approximate values of κ and q², for a given V₀, by considering the α-dependence of R²_mρ(0). Their results are applicable with the assumption that α ≫ 1. But, as we will see in section 4, if all values of α are possible a priori, α can be a bi-valued function of R²_mρ(0).