Confinement Induced Ordering in Fluid of Hard Ellipsoids

Han Miao，Hong-ru Maa.No.18 High School of Chongqing，Chongqing 400020，China .School of Mechanical Engineering，and Key Laoratory of Artificial Structures and Quantum Control（Ministry of Education），Shanghai Jiao Tong University，Shanghai 200240，China（Dated:Received on June 7，2015；Accepted on Octoer 18，2015）



Confinement Induced Ordering in Fluid of Hard Ellipsoids

Han Miaoa∗，Hong-ru Mab
a.No.18 High School of Chongqing，Chongqing 400020，China b.School of Mechanical Engineering，and Key Laboratory of Artificial Structures and Quantum Control（Ministry of Education），Shanghai Jiao Tong University，Shanghai 200240，China
（Dated:Received on June 7，2015；Accepted on October 18，2015）

The ordering configurations of a fluid of anisotropic ellipsoids under the confinement of two apposing impenetrable walls are studied by Monte Carlo simulations.The excess adsorption of the fluid on the walls with respect to the aspect ratio has a maximum at the critical aspect ratio of 2.9 in high-density ellipsoid fluids，indicating an orientational ordering in the adjacent region of the walls，which is confirmed by probing into the density configurations and the orientational order parameter in the adjacent region of the walls for varying aspect ratios.In addition，the orientational order parameter in the bulk fluid at the same density is calculated，and it indicates an isotropic state as the bulk density is still below the bulk isotropic-to-nematic transition.Therefore，it can be concluded that the anisotropic ordering near the walls in the ellipsoid fluid that exhibits isotropic in the bulk is induced by the confinement effect of the walls.

Confinement，Ellipsoid，Monte Carlo，Isotropic-to-nematic

I.INTRODUCTION

Hard-sphere（HS）fluid［1-7］is an ideal theoretical model that can be used to describe the fluid phase of matter，and it is also the simplest one that predicts a phase transition between a fluid phase and a solid phase.The HS model has been studied extensively for more than four decades，and thus deep understanding the various aspects of the model has been well established.Among these，one of the most important aspects is about the impact of the hard-wall confinement on the properties of the HS fluid.For example，G¨otzelmann et al.statistically computed the direct correlation function and the surface tension of uniform hard-sphere mixtures within the frame of the weighted density functional theory［8］.Later on，Roth et al.［9］studied the structures of polydisperse hard-sphere systems by using modified fundamental measure theory（FMT）originally proposed by Rosenfeld［10］，and predicted the existence of the oscillatory size segregation.On the other hand，Lin et al.experimentally investigated the hindered diffusion of an isolated Brownian particle system confined within two flat walls and implied that ignoring the pressure and velocity fields led to underestimation of the influence of the walls on the dragging force［11］.

Since 1960s an increasing interest has been focused on the study of isotropically interacting hard-body systems，which has led to a lot of achievements such as applications of the scaling particle theory［12］and a vast knowledge of correlation function［13］，whereas these studies are still based on single direct assumptions.In many real systems the effects caused by the anisotropic shape of particles cannot be ignored［14-19］.As a matter of fact，many anisotropic systems，such as rods［20］，spherocylinders［21］，and plates［22］have been proposed and examined in great detail.For example，Pfleiderer et al.［23］performed a numerical study on the stretchedfcc phase of hard ellipsoids of revolution and identified the SM2 phase for the range of aspect ratio 3.0-6.0. Recently，Cohen et al.［24］identified the structure of a simple dense fluid of ellipsoids using the direct confocal imaging technique that was dictated by the interplay between rotations and translations caused by the anisotropic shape.Cleaver et al.［25］simulated thin lyotropic liquid crystal films for which the anchoring strength and orientation were varied.Teixeira et al.［26］carried out a simulation study of the effects of nanoconfinement on a system of hard Gaussian overlap particles.These studies lead to the conclusion［27-32］that various particle shape and confinement have remarkable effects on the configurations of fluid.

To date，however，systems composed of monodisperse ellipsoids confined between two parallel hard walls have not been fully studied，and hence deserve further investigations［33-35］.Specifically，when the system is strongly confined in at least one direction，such as confined between two closely apposing walls or in a narrow cylindrical pore［36，37］，the confining perturbations may have significant effects on the structures of the anisotropic systems.Furthermore，hard ellipsoid as a typical anisotropic geometry can represent solidshapes from a hard disk to a needle by simply changing the aspect ratio，and thus the fluid of ellipsoids serves as an ideal model systems of anisotropic rigid particles.In order to deeply understand the effects of the rigid confinements on the various properties of an anisotropic ellipsoid system，we focus on the fluid of ellipsoids confined in one of the most representative confinements，i.e.two apposing flat walls that are readily accessed in experiments.

In this work，we systematically study the fluid of 3D hard prolate ellipsoids with various densities and aspect ratios confined between two impenetrable walls with varying distances using Monte Carlo simulations.

II.THEORY AND MODEL

We consider a model fluid consisting of N identical hard ellipsoids confined between two parallel impenetrable walls with the normal direction along the z axis located at the positions of z=0 and z=Lz，respectively. Reflective boundary condition is utilized at the z direction for the surfaces of the walls，while periodic boundary conditions are imposed on the x and y directions that have lengths of Lxand Ly，respectively.Accordingly，the bulk density is determined as ρb=N/LxLyLz. The geometry of the ellipsoid is specified by the lengths of three principal semi axes a，b，and c，where b is set to be equal to c and a is chosen as the longest semi axis in our simulations，and thereby the aspect ratio of the ellipsoid is simply characterized by a/b.Three unit vectors l，m，and n along the principal axes are implemented to describe the orientation of an ellipsoid，and hence the vector r indicating the position of an arbitrary point on the surface of the ellipsoid satisfies

For a fluid system consisting of anisotropic particles，there are two interchangeable types of entropy in the simulation system［38］，i.e.the translational entropy Etwhich denotes how many states the structure possesses，and the rotational entropy Erwhich is related to the number of distinct states that the particles can be arranged in the same structure.The free energy F in low-density systems is dominated by the translational entropy Et（F→Et），and thus in order to minimize F，the translational entropy Etshould be maximized.In other words，the value of Etis much larger than that of Erin a low density system，and therefore a change in Ercan rarely influence Et.In a marked contrast，a change in Ercan drastically affect Etin a high density system，in order to maximize Etthe anisotropic particles are aligned along a predominant direction facilitating the transformation of Erinto Et.Moreover，the conversion of Erinto Etresults in the loss of itself，which is proceeded via one of the following two possible mechanism∶（i）Lpis realized by the formation of an ordered nematic phase with the long axis of ellipsoids aligning parallel to the walls，the depletion potential of the walls induces more ellipsoids to be aligned parallelly to the walls［39］，（ii）Lvis proceeded by the formation of another ordered nematic phase with the long axis of ellipsoids aligning vertical to the walls［40］.The depletion potential of ellipsoids that results in the alignment of more ellipsoids to be vertical to the walls.So for the sake of maximizing Etwe need to compare Lpwith Lvto determine which mechanism is dominant.

If Lpis more dominant than Lv，the long axes of ellipsoids prefer to be ordered parallelly to the wall，indicating a stronger directing effect from the depletion potential of the walls.Otherwise the long axes of ellipsoids prefer to be aligned normal to the walls，indicating a stronger directing effect from the depletion potential of the ellipsoids.If Lpand Lvare comparable in magnitude，the alignments of the ellipsoids can be either parallel or vertical to the walls.

In this anisotropic system，one interesting quantity is the reduced excess adsorption［41，42］，which is defined as

where the absolute value of Γ represents the adsorption capability of the walls to the ellipsoids，z is the centerof-mass of the ellipsoid，ρ∗（z）is the local reduced number density，and ρ∗bis the reduced number density in the bulk.The reduced number density is defined as the number density multiplied by the volume V∗of the cube that circumscribes an ellipsoid，i.e.ρ∗=ρV∗.V∗is set to be the unit of volume in the simulations to guarantee that the unit of length in every simulation is adjusted with the aspect ratio to conserve V∗.The local orientation order parameter S，defined as the average of the largest eigenvalue of the Q tensor［43］is

which is another critical quantity in the anisotropic system that quantifies the orientational order.Here（α=x，y，z）is the Cartesian component of the longest axial unit vector of ellipsoid i.The value of the order parameter is close to zero for an isotropic phase，or close to 1 for a nematic phase in a finite system［35］.

The NV T Monte Carlo simulation is performed to obtain the equilibrium properties［44］.Hard-core interactions are implemented between ellipsoids as well as between ellipsoids and the walls，i.e.the interaction is zero when there is no any overlap between ellipsoids，or between the ellipsoids and the walls，otherwise the interaction is set to be infinite.

In the simulations，different sampling algorithms introduce errors to different degrees.In molecular dynamics simulations，a basic task is to estimate whetherthe results can be validated or not［45，46］.Instead，in Monte carlo simulations［47］，an importance-sampling algorithm introduced by Metropolis et al.is often adopted，in which every accessible point in configurational space is considered to be reached in a finite number of Monte Carlo steps from any given point.Moreover，many efficient Monte Carlo simulations have not been proven to be ergodic［44］.In order to lower the systematic errors in the biased sampling algorithms，we carry out the simulations with a large number of ellipsoidal particles，N=1000.All results are obtained for prolate ellipsoids with aspect ratios a/b of 1.0-3.5.Each simulation consists of a succession of trial moves where the average acceptance of translational/orientational moves is about 25%-30%.Each trial move of a chosen ellipsoid consists of a small random translational displacement and a small angle of random rotation.If the trial move results in any overlap between the moved ellipsoid with others，it is rejected，and otherwise，it is accepted.It takes 106MC steps per particle for equilibration and 2×106MC steps for data collection.

FIG.1 Plots of local density distributions with respect to z for various aspect ratios.The bulk density is fixed as ρb=0.3.

III.RESULTS AND DISCUSSION

Figure 1 shows the local density profiles of the ellipsoids with respect to the distance z to the wall at z=0 for four different aspect ratios a/b=1.5，2.4，2.7，and 3.0，and a fixed bulk density of ρb=0.3.There are two obvious features in the four density profiles.The first one is that the density vanishes when approaching z=b，which is dictated by the excluded-volume effect of the hard ellipsoid.The other is that there is a remarkable peak between z=b and z=a in each density profile，but the peak width and height are sensitive to the aspect ratio，i.e.that the peak is widened while the peak height is lowered as the aspect ratio increases.The presence of the peak indicates the layering arrangement of particles induced by the adsorption of the walls，and the variation tendency is dictated by the orientational configurations of the ellipsoids from parallel to normal to the walls that weakening the layering effect.Accordingly，the magnitude of oscillation is also reduced by increasing the anisotropy of ellipsoids.For this low density system of ρb=0.3，the bulk state is an isotropic fluid.Limited local densities when z is close to b for all four aspect ratios suggest that few ellipsoids are aligned with their long axes parallel to the walls.

Based on these density distributions，the surface excess adsorption is calculated for various aspect ratios and plotted in Fig.2.We find that the reduced excess adsorption Γ increases（or its magnitude decreases）with the aspect ratio in Fig.2，and we ascribe this to the orientational ordering effect of the two walls to the ellipsoids.Specifically，when the aspect ratio a/b=1.0，which leads to the isotropic hard sphere system，our result is similar to that reported by Roth et al.［34］.In general，the reduced surface excess adsorption is rather weak at the considered range of aspect ratios in Fig.2，which simply indicates that the effect of the wall con-finement is weak in the anisotropic fluid of ellipsoids at the low density.The weak excess adsorption is attributed to the factor that the translational entropy Etis expected to be more dominant over the rotational entropy Erat ρb=0.3［48］，and thereby the ellipsoids prefer to be more uniformly distributed in the entire available space.

FIG.2 Reduced surface excess adsorption as a function of the aspect ratio for the fixed bulk density of ρb=0.3.

FIG.3 Reduced excess adsorption with respect to the aspect ratios for the bulk densities ρb=0.4 and ρb=0.6，respectively.

The excess adsorption curves with respect to the aspect ratio of ellipsoids for two higher bulk density of ρb=0.4 and ρb=0.6 are present in Fig.3.For ρb=0.4，similar to that of ρb=0.3，the excess adsorption monotonically increases as the aspect ratio is increased.Surprisingly，when the bulk density is increased to ρb=0.6，the excess adsorption does exhibit a monotonic behavior as the aspect ratio varies from that at two lower bulk densities of ρb=0.3 and ρb=0.4.Indeed，it increases at the range of small aspect ratios until it reaches a maximum point at a/b≈2.9，and then it abruptly drops off. At the aspect ratio of a/b=3.0，the magnitude of the excess adsorption Γ has dropped down to Γ≈-0.076.

In order to understand this underlying mechanism in the non-monotonic variation behavior of the excess adsorption with respect to the aspect ratio，we probe into the internal configurations of the fluid at the high bulk density of ρb=0.6 for various aspect ratios，i.e.the inhomogeneous density profiles induced by the perturbation of the walls.Figure 4 shows four typical density distributions obtained for distinct aspect ratios located around the maximal point of the excess adsorption，a/b=2.7，2.8，2.9，and 3.0.Generally，each density distribution exhibits two notable peaks.One peak is located at the position very close to z=b，which is contributed by the ellipsoids aligned parallelly to the wall. The other peak is located at the position slightly further than z=a，and hence it is resulted in by the aggregation of the ellipsoids with the along axes normal to the wall.Importantly，when the aspect ratio increases from a/b=2.7 smaller than the critical value of 2.9 to be 3.0，the magnitudes of the two peaks change drastically.On the one side，the two peaks are maintained until they merges into a unique peak in the limit case of isotropic hard spheres as the aspect ratio is decreased.On the other side，the first peak is reduced indicating that less ellipsoids are aligned with the long axes parallel to the wall，whereas the second peak is strengthened indicating that more ellipsoids with the long axes normal to the wall are aggregated to the surface with a distance of around the length of the semi-major axis to the wall. In particular，for a/b=2.8，the first peak is only slightly weaker than the second one，while for a/b=2.9，the first peak has changed to be significantly weaker than the second one.For further increased a/b=3.0，the first peak has become very weak compared with the second peak.The varying tendency of the magnitudes of the two peaks when the aspect ratio is increased from a/b=2.7 going through the critical value of 2.9 to 3.0 leads to an obvious conclusion that the aggregation of ellipsoids near the wall is predominantly governed by the ellipsoids aligned with the long axes normal to the walls when the aspect ratio exceeds 2.9 in the highdensity ellipsoid fluid confined between two walls.

In order to understand the variation of the peaks further，we estimate the probability of the ellipsoids aligned in the angle interval P（θ＜θ0），where θ is the angle between the long axis of the ellipsoid and the z direction，i.e.，θ=0 means that the long axis of the ellipsoid is along the z axis while θ=π/2 means that the ellipsoid is aligned in the plane of the wall.Here we choose θ0=π/5 as a specific value to represent the case of small angles，and the data are collected in the slit area with the distance to the wall from b to 1.4a.The probability of P（θ＜θ0），which measures the aligning degree of the ellipsoids with the long axes normal to the walls，is plotted in Fig.5.We find that the probability increases slowly as the aspect ratio is till α=2.9，and then the probability increases rapidly.This indicates that more and more ellipsoids tend to be aligned perpendicularly to the wall when the aspect ratio exceeds 2.9，which results in the abrupt decrease（or increase in the magnitude）of the excess adsorption observed in Fig.3.

At a dense density ρb=0.6，the rotational entropy Eris not negligible［34］，and can be transformed into Et. Thus，the directed ordering of the ellipsoids near thewalls is dictated by the delicate balance between Lpand Lv.For aspect ratios smaller than 2.9，Fig.4 implies that the competition between Lpand Lvis too fierce and Lpis comparable to Lp，and thereby the ellipsoids are aligned either parallelly or vertically to the walls.However，for aspect ratios larger than 2.9，there exists a single dominant peak because Lvis much larger than Lp，indicating that the vertical configurations are dominant.In fact，Fig.4 also suggests that the depletion potential between the ellipsoids is more pronounced that the depletion potential between the ellipsoids and the walls，and as a result the ellipsoids prefer to stand on the walls.

FIG.4 Typical density distributions with respect to z for four aspect ratios a/b=2.7，2.8，2.9，and 3.0.The bulk density is fixed as ρb=0.6.

FIG.5 Probability of the ellipsoids aligned normal to the walls within a small angle interval θ＜π/5 as a function of the aspect ratio for the bulk density of ρb=0.6.The data are collected in the slit area with the distance to the wall from b to 1.4a.

Inspired by the variation of the ordering configurations near two walls sensitive to the aspect ratio，we examine the order parameter S of the hard ellipsoids within the wall confinement for two aspect ratios a/b=2.8 and 3.0，which are smaller and larger than the critical value of 2.9，respectively.The results of orientational order parameters with respect to the bulk density in Fig.6 reveal that the variation of orientational order parameters are also sensitive to the aspect ratio near the critical value of a/b=2.9.For a/b=2.8 smaller than the critical value，the orientational order parameter maintains roughly a constant low value of 0.1 indicating a near isotropic state.However，for a/b=3.0 larger than the critical value，S increases almost linearly as the bulk density increasing.For example，for ρb=0.7，S is increased to be 0.55，which indicates an isotropic-nematic transition near the walls.This observation supports the results further in Fig.4.For instance，for a/b=2.8 at ρb=0.7，the low orientational order suggests the the ellipsoids do not have a preferred alignment between parallel and vertical to the walls，indicating that Lvand Lpare comparable.Whereas，when the aspect ratio becomes a/b=3.0，the ellipsoids are directed to be ordered by the walls，and accordingly，Lvbecomes more dominant than Lp.

In order to distinguish the spontaneous ordering in the bulk anisotropic system from the induced ordering by the wall confinement，we compare their orientational order parameters of varying bulk densities for a specificaspect ratio of a/b=3.0.In Fig.7，we can see that S in the bulk systems only slightly increases as the bulk density increases and sustains a small value indicating that the bulk system of the ellipsoids with a/b=3.0 is truly isotropic.Furthermore，it has been known that the transition from isotropic to nematic in the bulk system occurs above the bulk density of ρb=0.7［49-53］. Therefore we are able to conclude that the anisotropic ordering behavior in the confined system arises from the confinement effect of the walls.

FIG.6 Orientational order parameter S with respect to the bulk density ρbfor two typical aspect ratios a/b=2.8 and 3.0，respectively.

FIG.7 Orientational order parameter S of the ellipsoids confined between two walls or in the bulk with respect to the bulk density for a fixed aspect ratio of a/b=3.0.

IV.CONCLUSION

We study the structure of a hard ellipsoid fluid confined between two hard walls by using the NV T Monte Carlo simulation.Our results suggest that the reduced surface excess adsorption is a sensitive parameter to identify the transition between an isotropic state and a nematic phase.Specifically，for the high bulk density ρb=0.6，but below the bulk transition density from isotropic to nematic，the reduced excess adsorption surprisingly exhibits a non-monotonic change as the aspect ratio varies，i.e.that it increases（or decreases in magnitude）at small aspect ratios until reaches a maximum at a critical aspect ratio of a/b=2.9，and then decreases abruptly.This finding is confirmed by the analysis of different relevant quantities，the density configurations and the orientational order parameter for various aspect ratios.Based on these consistent results we conclude that when the aspect ratio a/b≥2.9 the confinement effect of the walls induces an ellipsoid fluid to be ordered with the the long axes of ellipsoids aligned normally to the walls，while the walls have minor directing effect on the ordering of the ellipsoids when a/b＜2.9.This will place an useful，tight theoretical constraint on the investigation of an isotropic-nematic phase transition in the future.

V.ACKNOWLEDGEMENTS

This work is supported by the National Natural Science Foundation of China（No.10874111，No.11304169，and No.11174196）.

［1］D.E.Sullivan and G.Stell，J.Chem.Phys.69，5450（1978）.

［2］J.C.Neu，Phys.Rev.Lett.82，1072（1999）.

［3］C.N.Patra and S.K.Ghost，J.Chem.Phys.106，2752（1997）.

［4］R.D.Groot，N.M.Faber，and J.P.van der Eerdenn，Mol.Phys.62，861（1987）.

［5］V.M.Bedanov and F.M.Peeters，Phys.Rev.B 49，2667（1994）.

［6］V.A.Schweigert and F.M.Peeters，Phys.Rev.B 51，7700（1995）.

［7］R.Zangi，J.Phys.:Condens.Matter 16，S5371（2004）.

［8］B.G¨otzelmann，A.Haase，and S.Dietrich，Phys.Rev. E.53，3456（1996）.

［9］R.Roth，R.Evans，A.Lang，and G.Kahl，J.Phys.: Condens.Matter 14，12063（2002）.

［10］Y.Rosenfeld，Phys.Rev.Lett.63，980（1989）.

［11］B.H.Lin，J.Yu，and S.A.Rice，Phys.Rev.E 62，3909（2000）.

［12］H.Reiss，H.L.Frisch，E.Helfand，and J.L.Lebowitz，J.Chem.Phys.32，119（1960）.

［13］P.Tarazona，Phys.Rev.A 31，2672（1985）.

［14］A.Donev，J.Burton，F.H.Stillinger，and S.Torquato，Phys.Rev.B 73，054109（2006）.

［15］T.Schilling，S.Pronk，B.M.Mulder，and D.Frenkel，Phys.Rev.E 71，036138（2005）.

［16］R.Blaak，B.M.Mulder，and D.Frenkel，J.Chem.Phys. 120，5486（2004）.

［17］Y.Li，H.Miao，H.R.Ma，and J.Z.Y.Chen，Soft Matter 9，11461（2013）.

［18］C.Avendano and F.A.Escobedo，Soft Matter 8，4675（2012）.

［19］K.Wojciechowski and D.Frenkel，Comp.Met.Sci. Technol.10，235（2004）.

［20］M.A.Bates and D.Frenkel，J.Chem.Phys.112，10034（2000）.

［21］P.G.Bolhuis and D.Frenkel，J.Chem.Phys.106，666（1997）.

［22］A.G.Moreira and R.R.Netz，Euro.Phys.J.E 8，33（2002）.

［23］P.Pfleiderer and T.Schilling，Phys.Rev.E 75，020402（R）（2007）.

［24］A.P.Cohen，E.Janai，E.Mogilko，A.B.Schofield，and E.Sloutskin，Phys.Rev.Lett.107，238301（2011）.

［25］D.J.Cleaver and P.I.C.Teixeira，Chem.Phys.Lett. 338，1（2001）.

［26］P.I.C.Teixeira，F.Barmes，F.C.Anquetil-Deck，and D.J.Cleaver，Phys.Rev.E 79，011709（2009）.

［27］F.Barmes and D.J.Cleaver，Phys.Rev.E 71，021705（2005）.

［28］F.Barmes and D.J.Cleaver，Phys.Rev.E 69，061705（2004）.

［29］M.Moradi，R.J.Wheatley，and A.Avazpour，Phys. Rev.E 72，061706（2005）.

［30］P.I.C.Teixeira，A.Chrzanowska，G.D.Wall，and D. J.Cleaver，Mol.Phys.99，889（2001）.

［31］G.W.Wall and D.J.Cleaver，Mol.Phys.101，1105（2003）.

［32］G.W.Wall and D.J.Cleaver，Phys.Rev.E 56，4306（1997）.

［33］W.H.Li and H.R.Ma，J.Chem.Phys.119，585（2003）.

［34］D.Frenkel，B.M.Mulder，and J.P.McTague，Phys. Rev.Lett.52，287（1984）.

［35］J.Talbot，D.Kivelson，M.P.Allen，G.T.Evans，and D.Frenkel，J.Chem.Phys.92，3048（1990）.

［36］P.A.Kralchevsky and K.Nagayama，Adv.Colloid Interface Sci.85，145（2000）.

［37］M.Breuer，J.Bernsdorf，T.Zeiser，and F.Durst，Internat.J.Heta Fluid 21，186（2000）.

［38］M.S.Searle and D.H.Williams，J.Am.Chem.Soc. 114，10690（1992）.

［39］H.Miao，Y.Li，and H.R.Ma，J.Chem.Phys.140，154904（2014）.

［40］B.G¨otzelmann，R.Roth，S.Dietrich，M.Dijkstra，and R.Evans，Europhys.Lett.47，398（1999）.

［41］R.Roth and S.Dietrich，Phys.Rev.E 62，6926（2000）.

［42］J.R.Henderson and F.van Swol，Mol.Phys.51，991（1984）.

［43］U.Fabbri and C.Zannoni，Mol.Phys.58，763（1996）.

［44］D.Frenkel and B.Smit，Understanding Molecular Simulation，New York:Academic Press，A Division of Harcourt Inc.，（2002）.

［45］W.F.van Gunsteren and A.E.Mark，J.Chem.Phys. 108，6109（1998）.

［46］W.F.van Gunsteren，J.Dolenc，and A.E.Mark，Curr. Opin.Struct.Biol.18，149（2008）.

［47］X.G.Hong and K.Tamura，Chin.Phys.Lett.20，1315（2003）.

［48］D.Frenkel and B.M.Mulder，Mol.Phys.55，1171（1985）.

［49］X.Zheng，W.Iglesias，and P.Palffy-Muhoray，Phys. Rev.E 79，057702（2009）.

［50］M.Baus，J.L.Colot，X.G.Wu，and H.Xu，Phys.Rev. Lett.59，2184（1987）.

［51］R.Holyst and A.Poniewierski，Mol.Phys.68，381（1989）.

［52］U.P.Singh and Y.Singh，Phys.Rev.A 33，2725（1986）.

［53］J.F.Marko，Phys.Rev.Lett.60，325（1988）.

∗Author to whom correspondence should be addressed.E-mail：miaohan1105@126.com