{'fraction_non_alphanumeric': 0.039510075069142635, 'fraction_numerical': 0.01281255291527911, 'mean_word_length': 4.7228682170542635, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Strongly coupled Yukawa plasma layer in a harmonic trap

Observations made in dusty plasma experiments suggest that an ensemble of electrically charged solid particles, confined in an elongated trap, develops structural inhomogeneities. With narrowing the trap the particles tend to form layers oriented parallel with the trap walls. In this work we present theoretical and numerical results on the structure of three-dimensional many-particle systems with screened Coulomb (Yukawa) inter-particle interaction in the strongly coupled liquid phase, confined in one-dimensional harmonic trap, forming quasi-2D configurations. Particle density profiles are calculated by means of the hypernetted chain approximation (HNC), showing clear signs of layer formation. The mechanism behind the formation of layer structure is discussed and a method to predict the number of layers is presented. Molecular dynamics (MD) simulations provide validation of the theoretical results and detailed microscopic insights.

Strongly coupled charged particle systems (e.g. plasmas), where the inter-particle interaction energy exceeds the thermal kinetic energy of the constituent particles, can often be approximated by the one component plasma (OCP) model [1]. In cases when the dynamics of one dominant particle species decouples from that of the other components, it might be sufficient to provide a detailed description only for the former, while approximating the contribution of rest with a continuous neutralizing background. In the case of warm dense matter, for instance, the ions can follow classical particle trajectories, while the electrons experience quantum degeneracy and realize a homogeneous background [2,3]. In the case of dusty plasmas it is the very different time scales and charge states that decouple the dust dynamics (with characteristic times in the order of 10 ms and 10 4 elementary charges per dust) from the microscopic interactions with the electrons and ions present in the gas discharge (with typical times in the ns to µs range and unit charges). In both cases, only one of the plasma components (the ions in warm dens matter, and the dust particles in dusty plasmas) needs to be traced explicitly, the contribution of the background is reflected in the particular shape of the inter-particle interaction. If the polarizability of the isotropic neutralizing background is taken into account, the system can be approximated by the Yukawa one component plasma (YOCP) model. In this case the electrostatic pairwise interaction between the particles becomes screened, which screening can be approximated by the Dedye-Hückel mechanism resulting in an exponential decay superimposed to the bare Coulomb potential.
Since the pioneering experiments published in 1994 [4][5][6] laboratory dusty plasmas have been extensively used to gain insight into the microscopic details of macroscopic and collective phenomena. In large area radio frequency (RF) discharges, using micrometer sized monodisperse spherical dust particles it is easy to form a single layer of highly charged dust grains. This ensemble can form large (consisting of tens of thousands of dust grains) structures with the particles ordered in triangular lattices. Simply by changing discharge conditions (e.g. the RF power) or by other forms of external energy coupling (e.g. laser heating [7] or DC pulsing [8]) the system can undergo a solid to liquid transition and stabilize in the strongly coupled liquid state. However, the apparent single layer structure does not mean that the system is strictly two-dimensional. The vertical confinement is defined by the interplay of gravity and the electric field in the RF sheath, ultimately forming a potential well experienced by each dust particle. Both the vertical equilibrium position and the effective trap frequency depend on the charge-to-mass (q/m) ratio of the dust grains and the discharge conditions. At finite temperature the dust particles experience small amplitude vertical oscillations, forming quasi-2D configurations.
Early experimental investigations on quasi-2D system have been published by Lin I [9,10] observing anomalous diffusion, as well as layering and slow dynamics. Theoretical work on quasi-2D systems include studies of in-plane and out-of-plane polarized wave propagation and structural transitions into multi-layered configurations [11][12][13][14][15] building heavily on numerical simulations. An analytical approach targeting the particle density distribution in the trap would provide deeper insight into the double or triple layer formation, as observed in the experiments. A similar problem has been studied by Wrighton [16][17][18] for spherical confinement comparing mean field and hypernetted chain calculations to numerical results. In Klumov's work [19], crystallization of quasi-2D system under both harmonic and hard wall confinement was studied, in which the wetting phenomenon was observed in the hard wall case. Further theoretical studies on radially, as well as linearly confined liquids [20][21][22][23][24][25] include the pre-diction of the radial density distribution, mechanisms of oscillatory density profiles and ensuing solvation forces, in both the weak and strong coupling regimes.
In this work, we study the evolution of the density profile of quasi-2D Yukawa systems in a one-dimensional harmonic trap in the strongly coupled liquid state. The system is infinite in the x and y directions, where the particles can move freely. A harmonic trapping potential is applied in the z direction. Let n s denote the surface density of the quasi-2D layer projected onto the x-y plane. In this case we define a as the Wigner-Seitz radius in the projected plane as a 2 = 1/(πn s ).
The Yukawa interaction potential energy and the harmonic trap potential are, respectively where r is the three-dimensional inter-particle distance, κ is the dimensionless Yukawa screening parameter, ω t is the trap frequency, q and m are the electric charge and mass, equal for all particles.
The strength of the electrostatic coupling can be characterized with the Coulomb coupling parameter, nominally expressing the ratio of the potential to kinetic energies per particle, and is defined as where the thermodynamic β = 1/(k B T ). For Yukawa systems an effective coupling Γ eff (κ) < Γ [12,26] can be introduced, that depends on the Yukawa screening parameter a does more accurately characterize the state of a particular system. The strongly coupled domain is defined by Γ 1. Using the nominal 2D plasma frequency ω 2 p = 2πq 2 n s /(ma) to define the frequency and time units we introduce the dimensionless parameter t characterizing the strength of the trapping potential, as Within the frame of the quasi-2D YOCP model the behavior of a system is fully determined by the three dimensionless parameters Γ, κ and t.
Both the equilibrium properties and the dynamics of the system have been studied by molecular dynamics (MD) simulation and by theoretical analysis. Here we report on the equilibrium studies. Our MD simulations trance the trajectories of 10 000 particles in an external trapping potential defined by V (z) and periodic boundary conditions in x and y in a cubic simulation domain. Initial positions are assigned to the particles based on a simple initial barometric estimate. Inter-particle forces are summed for all particles within a radius of R ≈ 44a for each particle. Before conducting any measurements the system is given enough time (approximately of 4 000 plasma oscillation cycles) to reach equilibrium during the initial thermalization phase using the velocity backscaling thermostat. This is verified by observing the temperature stability after the thermaization is turned off. During the measurements data is collected and averaged over approx. 2 000 plasma oscillation cycles (100 000 time-steps).
In order to obtain a theoretical description of the density distribution within the trap we invoke the density functional theory (DFT) [27,28]. By minimizing the grand potential, a general, but quite complex expression for n(r), the 3D number density has been reported by Evans [28] (eq. 26 in chapter 3), see also the Appendix of [16] in terms of c (r, r ; [n(r)]), the direct correlation function (DCF) of the system. The notation emphasizes that the DCF is a unique, albeit unknown functional of the density profile. c(r, r ) is also connected with h(r, r ), the pair correlation function through the Ornstein-Zernike equation (OZ).
The latter is related to the pair distribution function through g(r, r ) = h(r, r ) + 1.
The density profile is determined by the selfconsistency relation [27] n(r) ∝ exp −βV (r) + n(r )c 0 (r − r ; n) dr , (5) where c 0 is the DCF of a reference system with uniform density. The structure of eq. (5) tells us that c(r, r ) plays the role of the effective interaction potential in the system (note that ϕ(r) does not appear explicitly in eq. (5)), In order to solve c and h simultaneously, the OZ relation provides the first equation, the second relation is derived by applying the hypernetted chain (HNC) approximation [27,29]. HNC has been successfully used in various problems relating to strongly coupled Coulomb and Yukawa systems [30]. The HNC approximation is based on the neglect of the so-called "bridge" (or irreducible, i.e. not derivable from the combined operations of "parallel connection" and "series connection" of Mayer diagrams) diagrams: their contribution is assumed to be negligible in the case of long range potentials. As a result, one obtains the general basic relationship g(r, r ) = exp [−βϕ(r − r ) + h(r, r ) − c(r, r )] , (7) which in combination with the OZ equation (4), provides a solution for h(r, r ) and c(r, r ). Restricting the solutions to homogeneous and isotropic functions, DCF satisfies the simpler variant of the OZ equation In equation (7), when r is large enough, g(r) → 1, h(r) → 0, the expression reduces to −βϕ ( r − r ) = c(r − r ). Using this asymptotic formula for the whole range of r, one arrives to the mean field (MF) approximation. Another self-consistent approach to calculating DCF for homogeneous fluids was described in [31].
Since the system is uniform in the x-y plane, the density is non-uniform only along the z direction, the parametrization can be simplified to n(r) = n(z), with the normalization condition n s = n(z) dz. The integral in equation (5) can be split intoz part and radial part separately. We can write eq. (5) as n(z) ∝ exp −βV (z) + 2π n(z )c 0 (z, z , ρ ; n)ρ dρ dz .
(9) We can rewrite eq. (9) in terms of the dimensionless quantitiesñ = na 2 ,z = z/a,Ũ (z) = β Γ U (z), etc. as ρ being the 2D distance between two projected particle positions in the x-y plane. Then, provided that c 0 (r − r ; n) is known one can obtain the density profile by the iterative solution of eq. (10). In view of eq. (6) the simplest approximation for c 0 (z − z , ρ) is to ignore correlations and set ϕ eff (r, r ) = ϕ(r − r ), i.e. c 0 (r − r ) = −βϕ(r − r ) (11) This is tantamount to a mean field approximation. Substituting (11) into (10), (10) reduces to In the sequel we drop the symbol for the dimensionless quantities. It should be kept in mind that all the length variables are in the unit of the Wigner-Seitz radius a. In the following, for simplicity, we only discuss the results for κ = 0.4. The results of the MF calculation and their comparison with the results of the MD simulation are given in Fig. 1. If there is no particle-particle interaction, then Γ → 0, the profile is Gaussian n(z) ∝ exp[−Γt 2 z 2 ], mapping the Maxwell distribution of non-interacting particles. Proceeding to higher Γ values, we expect the MF method to gradually fail to reproduce the numerical data as it is only valid for weak coupling (low Γ, i.e. low density or high temperature). Inspecting the MF profiles in Fig. 1, covering moderate and strong Γ values, we observe that those reasonably match the MD results at low Γ values, while fail spectacularly at high Γ-s. In particular MF does not provide the non-monotonic behavior related to layer formation. In fact, it has been proved analytically that the MF density profiles have to be monotonic on both sides of the maximum [32], and thus MF is unable to predict the formation of multiple layers. Next we calculate the DCF from the HNC approximation, then derive the density profile. Then in eq. (8) is chosen as the average number density of the plasma between the two layer boundaries. The boundary is defined arbitrarily as the points where the density value is one percent of the maximum density in the layer. This, of course requires the knowledge of the density, whose determination is the purpose of the calculation. All this lends itself to an iteration scheme. The resulting protocol for the numerical calculation is portrayed in Fig. 2. The algorithm is based on work published in [33]. The major improvement of the HNC over the MF calculation is that it correctly reproduces the splitting of the system into multiple layers. Figs. 3 and 4 show com- parisons with MD data at moderately high coupling values for cases when multi-peak profiles form. The similar multi-peak profile was reported in previous works [20,23,24,[34][35][36][37].
Generally, both t and Γ can affects the density profile formation, but the acting mechanisms are different. From Figs. 3 and 4 we can see that when t decreases, the layer becomes wider, developing more peaks (layers) in the density profile. The coupling parameter Γ only affect the density modulation amplitude. The relation between the trapping strength and the density profile was also studied in [38,39].
MD results in Figs. 5 and 6 show density profiles for a set of Γ values. With increasing Γ, the system develops very sharp density peaks with deep minima between the peaks, a sign of the formation of crystal-like ordering (no exchange of particle between layers). As Γ decreases, the particles experience larger vertical oscillation amplitudes, and the sides of neighboring peaks in the density profile fuse. The position and the number of the density peaks is mostly independent of the coupling, the distribution in the liquid state resembles that of the solid with lower amplitude of the density modulation. Based on the idea that a multiple layers may be regarded as a slab extracted from a 3D lattice [40][41][42], one may determine the number of layers by a simple algorithm. Assuming there are N layers, for simplicity a planar square lattice is associated to each layer. In this case the unit square side length is where w is the distance between the two outermost peaks. Comparing x N and d N , the lowest value of N at which x N > d N is the prediction for the number of layers (see Table I). In [40] the total energy of the system with different number of layers, including the particle interaction and the trapping energy is calculated. Despite operating with different physical quantities, the principles of the two methods are equivalent. In [43], the shell structure of the density profile in a confined colloidal particle system was studied by Monte Carlo simulation. A similar configuration was studied in [22] applying both the MF and HNC formalism, however, as those results are restricted to the weak coupling regime no layer formation was reported. In that work the correlation energy, based on local density approximation (LDA), and with this the excess chemical potential was calculated. In anther study [44] the density profile of a polyelectrolyte system in a harmonic trap was calculated. A similar shell  structure formation was observed. In [45] a similar DFT method was applied on hard-core Yukawa dusty plasma in a spherically harmonic trap. In [46] a study on van der Waals fluids with different confining potentials was presented including the modulation of the density profiles.
In our calculations, relying on the HNC approximation, the stability of the solution scheme appears to be limited to Γ ≤ 300 (for κ = 0.4), at higher couplings no converged solution could be found. For this reason comparing the HNC and MD result at the strongest coupling is performed at Γ = 300. Fig. 7 shows that even though the HNC result underestimates the density modulation amplitude, it predicts correctly the number and position of the peaks.
In conclusion, in this paper, we have analyzed the density profile of quasi-2D Yukawa plasma in harmonic trap in the strongly coupled liquid regime. We utilize the HNC approximation to calculate the density profile, which matches MD simulation result very well at low to moderate couplings. At stronger coupling, near the solidification the HNC calculations underestimate the density modulation amplitudes. We have presented a method to predict the number of layers in the solid phase based on the comparison between the in-plane inter-particle distance and the inter-layer distance. We confirm the validity of this prediction in the liquid state as well.
An outlook into future research directions include: (1) theoretical models beyond the standard HNC approximation could be applied to derive the direct correlation function c(r); (2) a method for the more accurate incorporation of correlation effect could be worked out to extend applicability to higher Γ-s, where the system is in an intermediate state between liquid and solid phases; (3) application of the current model to anharmonic confinement potentials, like hard-wall traps.
The authors thank Jeff Wrighton's help on numerical calculation in this paper. The authors are grateful for financial support from the Hungarian Office for Research, Development, and Innovation NKFIH grants K-134462 and K-132158.