Strongly coupled Yukawa plasma layer in a harmonic trap
SStrongly coupled Yukawa plasma layer in a harmonic trap
Hong Pan, Gabor J. Kalman, and Peter Hartmann Department of Physics, Boston College, Chestnut Hill, Massachusetts, 02467, USA Institute for Solid State Physics and Optics, Wigner Research Centre for Physics, P.O.Box. 49, H-1525 Budapest, Hungary (Dated: January 25, 2021)Observations made in dusty plasma experiments suggest that an ensemble of electrically chargedsolid particles, confined in an elongated trap, develops structural inhomogeneities. With narrowingthe trap the particles tend to form layers oriented parallel with the trap walls. In this work we presenttheoretical and numerical results on the structure of three-dimensional many-particle systems withscreened Coulomb (Yukawa) inter-particle interaction in the strongly coupled liquid phase, confinedin one-dimensional harmonic trap, forming quasi-2D configurations. Particle density profiles arecalculated by means of the hypernetted chain approximation (HNC), showing clear signs of layerformation. The mechanism behind the formation of layer structure is discussed and a method topredict the number of layers is presented. Molecular dynamics (MD) simulations provide validationof the theoretical results and detailed microscopic insights.
PACS numbers: 52.27.Gr, 52.27.Lw
Strongly coupled charged particle systems (e.g. plas-mas), where the inter-particle interaction energy exceedsthe 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 onedominant particle species decouples from that of theother components, it might be sufficient to provide adetailed description only for the former, while approx-imating the contribution of rest with a continuous neu-tralizing background. In the case of warm dense matter,for instance, the ions can follow classical particle trajec-tories, while the electrons experience quantum degener-acy and realize a homogeneous background [2, 3]. Inthe case of dusty plasmas it is the very different timescales and charge states that decouple the dust dynamics(with characteristic times in the order of 10 ms and 10 elementary charges per dust) from the microscopic in-teractions with the electrons and ions present in the gasdischarge (with typical times in the ns to µ s range andunit charges). In both cases, only one of the plasma com-ponents (the ions in warm dens matter, and the dust par-ticles in dusty plasmas) needs to be traced explicitly, thecontribution of the background is reflected in the partic-ular shape of the inter-particle interaction. If the polariz-ability of the isotropic neutralizing background is takeninto account, the system can be approximated by theYukawa one component plasma (YOCP) model. In thiscase the electrostatic pairwise interaction between theparticles becomes screened, which screening can be ap-proximated by the Dedye-H¨uckel mechanism resulting inan exponential decay superimposed to the bare Coulombpotential.Since the pioneering experiments published in 1994 [4–6] laboratory dusty plasmas have been extensively usedto gain insight into the microscopic details of macroscopicand collective phenomena. In large area radio frequency(RF) discharges, using micrometer sized monodisperse spherical dust particles it is easy to form a single layerof highly charged dust grains. This ensemble can formlarge (consisting of tens of thousands of dust grains)structures with the particles ordered in triangular lat-tices. Simply by changing discharge conditions (e.g. theRF power) or by other forms of external energy cou-pling (e.g. laser heating [7] or DC pulsing [8]) the systemcan undergo a solid to liquid transition and stabilize inthe strongly coupled liquid state. However, the appar-ent single layer structure does not mean that the systemis strictly two-dimensional. The vertical confinement isdefined by the interplay of gravity and the electric fieldin the RF sheath, ultimately forming a potential well ex-perienced by each dust particle. Both the vertical equi-librium position and the effective trap frequency dependon the charge-to-mass ( q/m ) ratio of the dust grains andthe discharge conditions. At finite temperature the dustparticles experience small amplitude vertical oscillations,forming quasi-2D configurations.Early experimental investigations on quasi-2D systemhave been published by Lin I [9, 10] observing anomalousdiffusion, as well as layering and slow dynamics. Theoret-ical work on quasi-2D systems include studies of in-planeand out-of-plane polarized wave propagation and struc-tural transitions into multi-layered configurations [11–15]building heavily on numerical simulations. An analyticalapproach targeting the particle density distribution inthe trap would provide deeper insight into the double ortriple layer formation, as observed in the experiments.A similar problem has been studied by Wrighton [16–18]for spherical confinement comparing mean field and hy-pernetted chain calculations to numerical results. In Klu-mov’s work [19], crystallization of quasi-2D system underboth harmonic and hard wall confinement was studied,in which the wetting phenomenon was observed in thehard wall case. Further theoretical studies on radially, aswell as linearly confined liquids [20–25] include the pre- a r X i v : . [ phy s i c s . p l a s m - ph ] J a n diction of the radial density distribution, mechanisms ofoscillatory density profiles and ensuing solvation forces,in both the weak and strong coupling regimes.In this work, we study the evolution of the densityprofile of quasi-2D Yukawa systems in a one-dimensionalharmonic trap in the strongly coupled liquid state. Thesystem is infinite in the x and y directions, where the par-ticles can move freely. A harmonic trapping potential isapplied in the z direction. Let n s denote the surface den-sity of the quasi-2D layer projected onto the x - y plane.In this case we define a as the Wigner-Seitz radius in theprojected plane as a = 1 / ( πn s ).The Yukawa interaction potential energy and the har-monic trap potential are, respectively ϕ ( r ) = q e − κ r/a r , and V ( z ) = mω z , (1)where r is the three-dimensional inter-particle distance, κ is the dimensionless Yukawa screening parameter, ω t isthe trap frequency, q and m are the electric charge andmass, equal for all particles.The strength of the electrostatic coupling can be char-acterized with the Coulomb coupling parameter, nomi-nally expressing the ratio of the potential to kinetic en-ergies per particle, and is defined asΓ = βq /a, (2)where the thermodynamic β = 1 / ( k B T ). For Yukawasystems an effective coupling Γ eff ( κ ) < Γ [12, 26] canbe introduced, that depends on the Yukawa screeningparameter a does more accurately characterize the stateof a particular system. The strongly coupled domain isdefined by Γ (cid:29) ω =2 πq n s / ( ma ) to define the frequency and time units weintroduce the dimensionless parameter t characterizingthe strength of the trapping potential, as t = ω ω = mω a πq n s = mω a q . (3)Within the frame of the quasi-2D YOCP model thebehavior of a system is fully determined by the threedimensionless parameters Γ, κ and t .Both the equilibrium properties and the dynamics ofthe system have been studied by molecular dynamics(MD) simulation and by theoretical analysis. Here wereport on the equilibrium studies. Our MD simulationstrance the trajectories of 10 000 particles in an externaltrapping potential defined by V ( z ) and periodic bound-ary conditions in x and y in a cubic simulation domain.Initial positions are assigned to the particles based on asimple initial barometric estimate. Inter-particle forcesare summed for all particles within a radius of R ≈ a for each particle. Before conducting any measurements the system is given enough time (approximately of 4 000plasma oscillation cycles) to reach equilibrium duringthe initial thermalization phase using the velocity back-scaling thermostat. This is verified by observing thetemperature stability after the thermaization is turnedoff. During the measurements data is collected and aver-aged over approx. 2 000 plasma oscillation cycles (100 000time-steps).In order to obtain a theoretical description of the den-sity distribution within the trap we invoke the densityfunctional theory (DFT) [27, 28]. By minimizing thegrand potential, a general, but quite complex expressionfor n ( r ), the 3D number density has been reported byEvans [28] (eq. 26 in chapter 3), see also the Appendixof [16] in terms of c ( r, r (cid:48) ; [ n ( r )]), the direct correlationfunction (DCF) of the system. The notation empha-sizes that the DCF is a unique, albeit unknown func-tional of the density profile. c ( r, r (cid:48) ) is also connectedwith h ( r, r (cid:48) ), the pair correlation function through theOrnstein-Zernike equation (OZ). h ( r, r (cid:48) ) = c ( r, r (cid:48) ) + (cid:90) c ( r, r (cid:48)(cid:48) ) n ( r (cid:48)(cid:48) ) h ( r (cid:48)(cid:48) , r (cid:48) ) d r (cid:48)(cid:48) (4)The latter is related to the pair distribution functionthrough g ( r, r (cid:48) ) = h ( r, r (cid:48) ) + 1.The density profile is determined by the self-consistency relation [27] n ( r ) ∝ exp (cid:20) − βV ( r ) + (cid:90) n ( r (cid:48) ) c ( r − r (cid:48) ; n ) d r (cid:48) (cid:21) , (5)where c is the DCF of a reference system with uniformdensity. The structure of eq. (5) tells us that c ( r, r (cid:48) ) playsthe role of the effective interaction potential in the system(note that ϕ ( r ) does not appear explicitly in eq. (5)), − βϕ eff ( r − r (cid:48) ) = c ( r − r (cid:48) ) (6)In order to solve c and h simultaneously, the OZ re-lation provides the first equation, the second relation isderived by applying the hypernetted chain (HNC) ap-proximation [27, 29]. HNC has been successfully used invarious problems relating to strongly coupled Coulomband Yukawa systems [30]. The HNC approximation isbased on the neglect of the so-called “bridge” (or irre-ducible, i.e. not derivable from the combined operationsof “parallel connection” and “series connection” of Mayerdiagrams) diagrams: their contribution is assumed to benegligible in the case of long range potentials. As a result,one obtains the general basic relationship g ( r, r (cid:48) ) = exp [ − βϕ ( r − r (cid:48) ) + h ( r, r (cid:48) ) − c ( r, r (cid:48) )] , (7)which in combination with the OZ equation (4), providesa solution for h ( r, r (cid:48) ) and c ( r, r (cid:48) ). Restricting the solu-tions to homogeneous and isotropic functions, DCF sat-isfies the simpler variant of the OZ equation h ( r ) = c ( r ) + ¯ n (cid:90) c ( | r − r (cid:48) | ) h ( r (cid:48) ) d r (cid:48) . (8)In equation (7), when r is large enough, g ( r ) → h ( r ) →
0, the expression reduces to − βϕ ( r − r (cid:48) ) = c ( r − r (cid:48) ). Us-ing this asymptotic formula for the whole range of r , onearrives to the mean field (MF) approximation. Anotherself-consistent approach to calculating DCF for homoge-neous fluids was described in [31].Since the system is uniform in the x - y plane, thedensity is non-uniform only along the z direction, theparametrization can be simplified to n ( r ) = n ( z ), withthe normalization condition n s = (cid:82) n ( z ) d z . The inte-gral in equation (5) can be split into z part and radialpart separately. We can write eq. (5) as n ( z ) ∝ exp (cid:20) − βV ( z ) + 2 π (cid:90) (cid:90) n ( z (cid:48) ) c ( z, z (cid:48) , ρ (cid:48) ; n ) ρ (cid:48) d ρ (cid:48) d z (cid:48) (cid:21) . (9)We can rewrite eq. (9) in terms of the dimensionlessquantities ˜ n = na , ˜ z = z/a , ˜ U ( z ) = β Γ U ( z ), etc. as˜ n (˜ z ) = ˜ n s exp[ − Γ ˜ U (˜ z )] (cid:82) ∞−∞ exp[ − Γ ˜ U (˜ z )] d˜ z (10)˜ U (˜ z ) = t ˜ z − ˜ W (˜ z )˜ W (˜ z ) = 2 π Γ (cid:90) (cid:90) ˜ n (˜ z (cid:48) ) c (˜ z, ˜ z (cid:48) , ˜ ρ (cid:48) ; n )˜ ρ (cid:48) d˜ ρ (cid:48) d˜ z (cid:48) ρ being the 2D distance between two projected par-ticle positions in the x - y plane. Then, provided that c ( r − r (cid:48) ; n ) is known one can obtain the density profileby the iterative solution of eq. (10). In view of eq. (6)the simplest approximation for c ( z − z (cid:48) , ρ ) is to ignorecorrelations and set ϕ eff ( r, r (cid:48) ) = ϕ ( r − r (cid:48) ), i.e. c ( r − r (cid:48) ) = − βϕ ( r − r (cid:48) ) (11)This is tantamount to a mean field approximation. Sub-stituting (11) into (10), (10) reduces to U ( z ) = t z + 2 (cid:82)(cid:82) n ( z (cid:48) ) exp[ − κd ( z,z (cid:48) ,ρ (cid:48) )] d ( z,z (cid:48) ,ρ (cid:48) ) ρ (cid:48) d ρ (cid:48) d z (cid:48) d ( z, z (cid:48) , ρ (cid:48) ) = (cid:112) ( z − z (cid:48) ) + ρ (cid:48) . (12)In the sequel we drop the symbol for the dimensionlessquantities. It should be kept in mind that all the lengthvariables are in the unit of the Wigner-Seitz radius a . Inthe following, for simplicity, we only discuss the resultsfor κ = 0 .
4. The results of the MF calculation and theircomparison with the results of the MD simulation aregiven in Fig. 1. If there is no particle-particle interaction,then Γ →
0, the profile is Gaussian n ( z ) ∝ exp[ − Γ t z ],mapping the Maxwell distribution of non-interacting par-ticles. Proceeding to higher Γ values, we expect the MFmethod to gradually fail to reproduce the numerical dataas it is only valid for weak coupling (low Γ, i.e. low den-sity or high temperature). Inspecting the MF profiles inFig. 1, covering moderate and strong Γ values, we observethat those reasonably match the MD results at low Γ val-ues, while fail spectacularly at high Γ-s. In particular MFdoes not provide the non-monotonic behavior related to layer formation. In fact, it has been proved analyticallythat the MF density profiles have to be monotonic onboth sides of the maximum [32], and thus MF is unableto predict the formation of multiple layers. FIG. 1. Density profile comparison between the MD simula-tion and the MF approximation (continuous line) for differentΓ and t values. (a) Γ = 8, t = 0 .
2, (b) Γ = 16, t = 0 .
1, (c)Γ = 64, t = 0 .
2, (d) Γ = 64, t = 0 . Next we calculate the DCF from the HNC approxima-tion, then derive the density profile. The ¯ n in eq. (8)is chosen as the average number density of the plasmabetween the two layer boundaries. The boundary is de-fined arbitrarily as the points where the density value isone percent of the maximum density in the layer. This,of course requires the knowledge of the density, whosedetermination is the purpose of the calculation. All thislends itself to an iteration scheme. The resulting proto-col for the numerical calculation is portrayed in Fig. 2.The algorithm is based on work published in [33]. FIG. 2. Numerical iteration loop for calculating density pro-file.
The major improvement of the HNC over the MF cal-culation is that it correctly reproduces the splitting ofthe system into multiple layers. Figs. 3 and 4 show com-
FIG. 3. Density profile comparison between MD simulation(continuous line) and the HNC approximation, for differenttrapping strengths. Γ = 32FIG. 4. Density profile comparison between MD simulation(continuous line) and the HNC approximation, for differenttrapping strengths. Γ = 64 parisons with MD data at moderately high coupling val-ues for cases when multi-peak profiles form. The sim-ilar multi-peak profile was reported in previous works[20, 23, 24, 34–37].Generally, both t and Γ can affects the density pro-file formation, but the acting mechanisms are different.From Figs. 3 and 4 we can see that when t decreases, thelayer becomes wider, developing more peaks (layers) inthe density profile. The coupling parameter Γ only affectthe density modulation amplitude. The relation betweenthe trapping strength and the density profile was alsostudied in [38, 39].MD results in Figs. 5 and 6 show density profiles for a set of Γ values. With increasing Γ, the system developsvery sharp density peaks with deep minima between thepeaks, a sign of the formation of crystal-like ordering (noexchange of particle between layers). As Γ decreases, theparticles experience larger vertical oscillation amplitudes,and the sides of neighboring peaks in the density profilefuse. The position and the number of the density peaksis mostly independent of the coupling, the distributionin the liquid state resembles that of the solid with loweramplitude of the density modulation.Based on the idea that a multiple layers may be re-garded as a slab extracted from a 3D lattice [40–42], onemay determine the number of layers by a simple algo-rithm. Assuming there are N layers, for simplicity a pla-nar square lattice is associated to each layer. In this casethe unit square side length is x N = √ πN . The inter-layerdistance is d N = wN − , where w is the distance betweenthe two outermost peaks. Comparing x N and d N , thelowest value of N at which x N > d N is the predictionfor the number of layers (see Table I). In [40] the totalenergy of the system with different number of layers, in-cluding the particle interaction and the trapping energyis calculated. Despite operating with different physicalquantities, the principles of the two methods are equiva-lent. In [43], the shell structure of the density profile in aconfined colloidal particle system was studied by MonteCarlo simulation. FIG. 5. Evolution of the density profiles with increasing Γvalues as determined by MD simulations. t = 0 . A similar configuration was studied in [22] applyingboth the MF and HNC formalism, however, as those re-sults are restricted to the weak coupling regime no layerformation was reported. In that work the correlation en-ergy, based on local density approximation (LDA), andwith this the excess chemical potential was calculated. Inanther study [44] the density profile of a polyelectrolytesystem in a harmonic trap was calculated. A similar shell √ Nπ N − κ = 0 .
4, the row with green color corresponds to thepredicted N value.FIG. 6. Evolution of the density profiles with increasing Γvalues as determined by MD simulations. t = 0 . t = 0 .
05, (b) t = 0 . structure formation was observed. In [45] a similar DFTmethod was applied on hard-core Yukawa dusty plasmain a spherically harmonic trap. In [46] a study on van derWaals fluids with different confining potentials was pre-sented including the modulation of the density profiles.In our calculations, relying on the HNC approxima-tion, the stability of the solution scheme appears to belimited to Γ ≤
300 (for κ = 0 . c ( r ); (2) a method for the more accurate in-corporation of correlation effect could be worked out toextend applicability to higher Γ-s, where the system isin an intermediate state between liquid and solid phases;(3) application of the current model to anharmonic con-finement potentials, like hard-wall traps.The authors thank Jeff Wrighton’s help on numericalcalculation in this paper. The authors are grateful forfinancial support from the Hungarian Office for Research,Development, and Innovation NKFIH grants K-134462and K-132158. [1] M. Baus and J.-P. Hansen, Statistical mechanics of sim-ple coulomb systems, Physics Reports , 1 (1980).[2] J. Cl´erouin, P. Arnault, C. Ticknor, J. D. Kress, andL. A. Collins, Unified concept of effective one compo-nent plasma for hot dense plasmas, Phys. Rev. Lett. ,115003 (2016).[3] Z.-Q. Wang, J. Tang, Y. Hou, Q.-F. Chen, X.-R. Chen,J.-Y. Dai, X.-J. Meng, Y.-J. Gu, L. Liu, G.-J. Li, Y.-S. Lan, and Z.-G. Li, Benchmarking the effective one-component plasma model for warm dense neon andkrypton within quantum molecular dynamics simulation,Phys. Rev. E , 023302 (2020).[4] J. H. Chu and L. I, Direct observation of coulomb crystalsand liquids in strongly coupled rf dusty plasmas, Phys.Rev. Lett. , 4009 (1994).[5] H. Thomas, G. E. Morfill, V. Demmel, J. Goree, B. Feuer-bacher, and D. M¨ohlmann, Plasma crystal: Coulombcrystallization in a dusty plasma, Phys. Rev. Lett. ,652 (1994).[6] A. Melzer, T. Trottenberg, and A. Piel, Experimentaldetermination of the charge on dust particles formingcoulomb lattices, Physics Letters A , 301 (1994).[7] V. Nosenko, J. Goree, and A. Piel, Laser method of heat-ing monolayer dusty plasmas, Physics of Plasmas ,032106 (2006).[8] Z. Donk´o, P. Hartmann, P. Magyar, G. J. Kalman, andK. I. Golden, Higher order structure in a complex plasma,Physics of Plasmas , 103701 (2017).[9] W. Juan and L. I, Anomalous diffusion in strongly cou-pled quasi-2d dusty plasmas, Phys. Rev. Lett. , 3073(1998).[10] L. Teng, P. Tu, and L. I, Microscopic observationof confinement-induced layering and slow dynamics ofdusty-plasma liquids in narrow channels, Phys. Rev. Lett. , 245004 (2003).[11] O. Bystrenko, Structural transitions in one-dimensionallyconfined one-component plasmas, Phys. Rev. E ,025401 (2003).[12] P. Hartmann, G. J. Kalman, Z. Donk´o, and K. Ku-tasi, Equilibrium properties and phase diagram of two-dimensional yukawa systems, Phys. Rev. E , 026409(2005).[13] K. Qiao and T. W. Hyde, Structural phase transi-tions and out-of-plane dust lattice instabilities in verti-cally confined plasma crystals, Phys. Rev. E , 026406(2005).[14] C. Henning, P. Ludwig, A. Filinov, A. Piel, andM. Bonitz, Ground state of a confined yukawa plasmaincluding correlation effects, Phys. Rev. E , 036404(2007).[15] Z. Donk´o, P. Hartmann, and G. J. Kalman, Two-dimensional dusty plasma crystals and liquids, Journalof Physics: Conference Series , 012016 (2009).[16] J. Wrighton, J. W. Dufty, H. K¨ahlert, and M. Bonitz,Theoretical description of coulomb balls: Fluid phase,Phys. Rev. E , 066405 (2009).[17] H. Bruhn, H. K¨ahlert, T. Ott, M. Bonitz, J. Wrighton,and J. W. Dufty, Theoretical description of sphericallyconfined, strongly correlated yukawa plasmas, Phys. Rev.E , 046407 (2011).[18] J. Wrighton, H. K¨ahlert, T. Ott, P. Ludwig, H. Thom- sen, J. Dufty, and M. Bonitz, Charge correlations in aharmonic trap, Contributions to Plasma Physics , 45(2012).[19] B. A. Klumov and G. Morfill, Effect of confinement onthe crystallization of a dusty plasma in narrow channels,JETP Letters , 409 (2008).[20] T. K. Vanderlick, L. E. Scriven, and H. T. Davis, Molec-ular theories of confined fluids, The Journal of ChemicalPhysics , 2422 (1989).[21] Y.-X. Yu, A novel weighted density functional theory foradsorption, fluid-solid interfacial tension, and disjoiningproperties of simple liquid films on planar solid surfaces,The Journal of Chemical Physics , 024704 (2009).[22] M. Girotto, A. P. dos Santos, T. Colla, and Y. Levin,Yukawa particles in a confining potential, The Journal ofChemical Physics , 014106 (2014).[23] K. Nyg˚ard, Local structure and density fluctuations inconfined fluids, Current Opinion in Colloid & InterfaceScience , 30 (2016).[24] G. A. Mansoori and S. A. Rice, Confined fluids: Struc-ture, properties and phase behavior, in Advances inChemical Physics (John Wiley & Sons, Ltd, 2014)Chap. 5, pp. 197–294.[25] C. Henning, H. Baumgartner, A. Piel, P. Ludwig, V. Gol-ubnichiy, M. Bonitz, and D. Block, Ground state of a con-fined yukawa plasma, Phys. Rev. E , 056403 (2006).[26] T. Ott, M. Bonitz, L. G. Stanton, and M. S. Murillo, Cou-pling strength in coulomb and yukawa one-componentplasmas, Physics of Plasmas , 113704 (2014).[27] Theory of simple liquids, in Theory of Simple Liquids(Fourth Edition) , edited by J. P. Hansen and I. R. Mc-Donald (Academic Press, Oxford, 2013) fourth editioned.[28] Fundamentals of inhomogenous fluids, in
Fundamentalsof Inhomogenous Fluids , edited by D. Henderson (MarcelDekker, NY, 1992).[29] J. S. Rowlinson, The equation of state of dense systems,Reports on Progress in Physics , 169 (1965).[30] K. Ng, Hypernetted chain solutions for the classical one-component plasma up to Γ = 7000, The Journal of Chem-ical Physics , 2680 (1974).[31] G. Rickayzen, P. Kalpaxis, and E. Chacon, A self-consistent approach to a density functional for homoge-neous fluids, The Journal of Chemical Physics , 7963(1994).[32] H. Pan, Static and dynamic properties of strongly coupledquasi-2D Yukawa plasma layers , Ph.D. thesis, MorrisseyCollege of Arts and Sciences, Boston College (2019).[33] J. F. Springer, M. A. Pokrant, and F. A. Stevens, Inte-gral equation solutions for the classical electron gas, TheJournal of Chemical Physics , 4863 (1973).[34] I. K. Snook and W. van Megen, Solvation forces in simpledense fluids. i, The Journal of Chemical Physics , 2907(1980).[35] T. S. Ingebrigtsen and J. C. Dyre, The impact range forsmooth wall–liquid interactions in nanoconfined liquids,Soft Matter , 4324 (2014).[36] C. Ebner, M. A. Lee, and W. F. Saam, Nonuniform one-dimensional classical fluids: Theory versus monte carloexperiments, Phys. Rev. A , 959 (1980).[37] E. Kierlik and M. L. Rosinberg, Density-functional the-ory for inhomogeneous fluids: Adsorption of binary mix-tures, Phys. Rev. A , 5025 (1991).[38] T. Ott, M. Bonitz, Z. Donk´o, and P. Hartmann, Superdif- fusion in quasi-two-dimensional yukawa liquids, Phys.Rev. E , 026409 (2008).[39] N. G. Almarza, C. Mart´ın, and E. Lomba, Phase be-havior of the hard-sphere maier-saupe fluid under spatialconfinement, Phys. Rev. E , 031501 (2009).[40] H. Totsuji, T. Kishimoto, and C. Totsuji, Structure ofconfined yukawa system (dusty plasma), Phys. Rev. Lett. , 3113 (1997).[41] H. Totsuji, T. Kishimoto, Y. Inoue, C. Totsuji,and S. Nara, Yukawa system (dusty plasma) in one-dimensional external fields, Physics Letters A , 215(1996).[42] E. C. Oguz, R. Messina, and H. L¨owen, Multilayeredcrystals of macroions under slit confinement, J. Phys.: Condens. Matter , 424110 (2009).[43] A. Ricci, P. Nielaba, S. Sengupta, and K. Binder, Order-ing of two-dimensional crystals confined in strips of finitewidth, Phys. Rev. E , 011405 (2007).[44] S. Dutta and Y. S. Jho, Shell formation in short like-charged polyelectrolytes in a harmonic trap, Phys. Rev.E , 012503 (2016).[45] F. Gu, H.-J. Wang, and J.-T. Li, Density functional the-ory for the ground state of spherically confined dustyplasma, Phys. Rev. E , 056402 (2012).[46] Y. Tang and J. Wu, Modeling inhomogeneous van derwaals fluids using an analytical direct correlation func-tion, Phys. Rev. E70