document: the debate about the properties of the melting phase transition of two - dimensional ( 2d ) systems did not lose its intensity over the past several decades . 
   recent developments in the fabrication of 2d materials @xcite simultaneously seek for , and may provide clarification of the details of the transition . 
   a milestone , and still the most widely accepted theory available , is the kosterlitz - thouless - halperin - nelson - young ( kthny ) picture @xcite . in the underlying physical process 
   two separate , continuous transitions can be distinguished , as the solid transforms into a liquid in quasi - equilibrium steps by slow heating . during the first stage the translational ( positional ) order vanishes , while in the second stage the orientational order decays . 
   all this is mediated by the unbinding of ( i ) dislocation pairs into individual dislocations , and ( ii ) dislocations into point defects @xcite . the strength of this theory consists in its compatibility with the mermin - wagner theorem that forbids the existence of exact long range positional order in 2d for a wide range of pair potentials , at finite temperatures @xcite . 
   the most criticized weakness of it , however , is that it assumes a dilute , unstructured distribution of the lattice defects , which is in contradiction with observations , where the alignment and accumulation of dislocations into small angle domain walls was found @xcite .    since the birth of the kthny theory , the examination of its validity for systems with different pair interactions has been in focus . 
   investigations started with hard - sphere ( disk ) , lennard - jones , and coulomb systems . 
   more recently , systems characterized by dipole - dipole and debye - hckel ( screened coulomb or yukawa ) inter - particle interactions became important due to the significant advances achieved in the field of colloid suspensions @xcite and dusty plasmas @xcite . 
   to illustrate the incongruity of both experimental and simulation results that had accumulated over the last three decades on investigations of classical single - layer ( 2d ) many - body systems , we list a few examples :    * first order phase transition to exist was reported for lennard - jones systems @xcite and hard - disk systems @xcite , for the phase - field - crystal ( pfc ) model @xcite , as well as for coulomb and dipole systems @xcite in simulations , and in experiments with halomethanes and haloethanes physisorbed on exfoliated graphite @xcite , as well as in experiments on a quasi - two - dimensional suspension of uncharged silica spheres @xcite . * 
   second order ( or single step continuous ) transition was found in dusty plasma experiments @xcite , for a hard - disk system @xcite , electro - hydrodynamicly excited colloidal suspensions @xcite , as well as for coulomb @xcite and yukawa @xcite systems . 
   * kthny - like transition was reported in a dusty plasma experiment @xcite and related numerical simulations @xcite , for the harmonic lattice model @xcite , in experiments and simulations of colloidal suspensions @xcite , for lennard - jones @xcite , yukawa @xcite , hard disk @xcite , dipole - dipole @xcite , gaussian - core @xcite , and electron systems @xcite , as well as for a system with @xmath1 repulsive pair potential @xcite , weakly softened core @xcite , and for vortices in a w - based superconducting thin film @xcite . 
   the effect of the range of the potential on two - dimensional melting was studied in @xcite for a wide range of morse potentials . 
   it has been shown , that extended - ranged interatomic potentials are important for the formation of a `` stable '' hexatic phase . 
   similar conclusion was drawn in @xcite for modified hard - disk potentials . 
   the effect of the dimensionality ( deviation from the mathematically perfect 2d plane ) on the hexatic phase was discussed for lennard - jones systems in @xcite . 
   it was found , that an intermediate hexatic phase could only be observed in a monolayer of particles confined such that the fluctuations in the positions perpendicular to the particle layer was less than 0.15 particle diameters . 
   the timeline of the results listed above shows a general trend : in earlier studies , first or second order phase transitions were identified in particle simulations , but subsequently , as the computational power increased with time , since approximately the year of 2000 , particle based numerical studies became in favor of the kthny theory . 
   a possible resolution of the ongoing debate is given in @xcite , where extensive monte carlo simulations of 2d lennard - jones systems have revealed the metastable nature of the hexatic phase . 
   this seems to support pfc simulations @xcite operating on the diffusive time - scale ( averaging out single particle oscillations ) , which is significantly longer , than what monte carlo ( mc ) or molecular dynamics ( md ) methods can cover .    in this paper we will show that the observation of the hexatic phase is strongly linked with the thermodynamic equilibration of the systems . 
   the necessary equilibration time , in turn , strongly depends on the measured quantity of interest . 
   local , or single particle properties can equilibrate very rapidly , while long - range , or collective relaxations usually take significantly longer . 
   we find , consequently , that monitoring the velocity distribution function alone to verify the equilibration of the system is insufficient . 
   the idea , that numerical simulations may have related equilibration issues ( called as kinetic bottlenecks ) was raised already in 1993 in @xcite . 
   we have performed extensive microcanonical md simulations @xcite in the close vicinity of the expected solid - liquid phase transition temperature , @xmath2 , for repulsive screened coulomb ( also called yukawa or debye - hckel ) pair - potential with the potential energy in form of @xmath3 where @xmath4 is the debye screening length , @xmath5 is the electric charge of the particles , and @xmath6 is the vacuum permittivity . to characterize the screening we use the dimensionless screening parameter @xmath7 , where @xmath8 is the wigner - seitz radius , and @xmath9 is the particle density 
   this model potential was chosen because of its relevance to several experimental systems consisting of electrically charged particles , like dusty plasmas , charged colloidal suspensions , and electrolytes . here 
   we show results obtained for @xmath10 . 
   our earlier studies @xcite identified the melting transition ( without clarifying its nature ) to take place around the coulomb coupling parameter @xmath11 for this strength of screening . 
   time is measured in units of the nominal 2d plasma oscillation period with @xmath12 where @xmath13 is the mass of a particle . 
   our simulations are initialized by placing @xmath14 particles ( in the range of 1,920 to 740,000 ) into a rectangular simulation cell that has periodic boundary conditions . 
   the particles are released from hexagonal lattice positions , with initial velocities randomly sampled from a predefined distribution . at the initial stage , which has a duration @xmath15 ( thermalization time ) , 
   the system is thermostated by applying the velocity back - scaling method ( to follow the usual approach used in many previous studies ) to reach near - equilibrium state at the desired ( kinetic ) temperature . 
   data collection starts only after this initial stage and runs for a time period @xmath16 ( measurement time ) without any additional thermostation .    to characterize the level of equilibration we study the time and system size dependence of the following quantities :    * momenta of the velocity distribution function , @xmath17 , * the configurational temperature , @xmath18 @xcite , and * the long - range decay of the @xmath19 pair - correlation and @xmath20 bond - angle correlation functions @xcite .    while in the case of the first two quantities @xmath21 , in the simulations targeting the correlation functions , @xmath15 is varied over a wide range and the measurement time 
   is chosen to be @xmath22 to avoid significant changes ( due to ongoing equilibration ) during the measurement .      using the maxwell - boltzmann assumption for the velocity distribution in thermal equilibrium in the form @xmath23 where @xmath24 , in two - dimensions the first four velocity moments are : @xmath25    to measure the relaxation time of the velocity distribution function we have performed md simulations with particle numbers @xmath26 and @xmath27 , with initial velocity components ( @xmath28 and @xmath29 ) sampled from a uniform distribution between @xmath30 and @xmath31 , in order to start with the desired average kinetic energy , but being far from equilibrium . 
   figure [ fig : vmom ] shows the time evolution of the first eight velocity moments normalized with their theoretical equilibrium values . 
   as already mentioned , the initial conditions are far from the equilibrium configuration ( perfect lattice position and non - thermal velocity distribution ) . 
   ( color online ) moments of the computed velocity distribution functions relative to the theoretical equilibrium values vs. simulation time at temperatures slightly above ( full lines ) and below ( dashed lines ) the melting point , @xmath2 . 
   the dashed lines are mostly hidden behind the full lines , indicating a low sensitivity on the temperature . 
   the dark red curve shows functional fit in the form @xmath32 to @xmath33 . 
   @xmath26.,width=302 ]    we can observe , that the velocity momenta have initial values very different from the expected maxwell - boltzmann equilibrium distribution . 
   the values approach the equilibrium value asymptotically with regular oscillations . these oscillations ( or fluctuations ) are typical for microcanonical md simulations , where the total energy of the system is constant , while there is a permanent exchange of potential and kinetic energies . 
   the relaxation time can be found by fitting the curves with an exponential asymptotic formula in the form @xmath32 . 
   we find , that the relaxation of the velocity distribution can be characterized by a short relaxation time of @xmath34 , and this is independent of system size and temperature in the vicinity of the melting point .      in 1997 , rugh @xcite pointed out that the temperature can also be expressed as ensemble average over geometrical and dynamical quantities and derived the formula for the configurational temperature : @xmath35 where @xmath36 . 
   as the central quantity in this expression is the inter - particle force acting on each particle , in case of finite range interactions ( like the yukawa potential ) , the configurational temperature is sensitive on the local environment within this range . 
   simulations were performed for a series of particle numbers between @xmath37 and @xmath38 with initial velocities sampled from maxwellian distribution . 
   figure  [ fig : tconf](a ) shows examples from runs with @xmath26 for the time evolutions , while fig . 
   [ fig : tconf](b ) presents relaxation time data computed ( similarly as above ) for different kinetic temperatures . 
   ( color online ) ( a ) time evolution of the configurational temperature @xmath18 for different kinetic temperatures @xmath39 . 
   ( b ) relaxation time vs. kinetic temperature . 
   ( @xmath26).,width=302 ]    we observe relaxation times about an order of magnitude longer ( @xmath40 ) compared to the velocity distribution , and a strong temperature dependence in the vicinity of the melting point . 
   no significant system size dependence was found . 
   the central property used to identify the hexatic phase is traditionally the long - range behavior of the pair- , and bond - order correlation functions , @xmath19 and @xmath20 , respectively @xcite . to be able to compute correlations at large distances , one naturally has to use large particle numbers , 
   otherwise the periodic boundary conditions introduce artificial correlation peaks . 
   this trivial constraint led to investigations of larger and larger systems by different groups . 
   figures  [ fig : corr ] and [ fig : corr2 ] show correlation functions for systems consisting of @xmath41 particles , for a set of increasing equilibration times provided to the systems before performing the data collection . 
   ( color online ) log  log plots of ( a ) an example of @xmath42 pair correlation function with its upper envelope , ( b ) a series of envelope curves of pair correlation functions , ( c ) @xmath20 bond - order correlation functions measured after letting the systems equilibrate for various times indicated , @xmath15 , at a temperature 1 percent above the melting point . 
   the systems consisted of @xmath41 particles , the data acquisition took @xmath43 and started after @xmath15 has elapsed.,width=302 ]     ( color online ) same as fig  [ fig : corr ] with semi - logarithmic scales.,width=302 ]    we can observe a clear long - time evolution of the correlation functions . on the double - logarithmic plot 
   the @xmath19 pair - correlation functions show already at early times a long - range decay , which is faster than power - law [ fig . 
   [ fig : corr](a , b ) ] , while the @xmath20 orientational correlations smooth out to near perfect straight lines [ fig . 
   [ fig : corr](c ) ] , representing power - law type decay for relatively long times . on the semi - logarithmic graphs all the @xmath19 functions have almost straight upper envelopes [ fig . 
   [ fig : corr2](a , b ) ] in the intermediate distance range @xmath44 , where the statistical noise is still negligible . 
   this indicates almost pure exponential decay , although the characteristic decay distance does decrease with increasing simulation time . on the other hand , it is only the last @xmath20 orientational correlation function , belonging to the longest simulation , that shows linear apparent asymptote on the semi - logarithmic scale [ fig . 
   [ fig : corr2](c ) ] , representing a clear exponential decay , meaning the lack of long range order . to conclude these observations : in short simulations we observe short - range positional and quasi - long - range orientational order , signatures of the hexatic phase , which , however vanish if we provide the system longer time for equilibration . as a consequence , in the case we would stop the simulation at , e.g. , @xmath45 ( which already means simulation time - steps in the order of @xmath46 , as plasma oscillations have to be resolved smoothly ) we may identify the system to be in the hexatic phase , exactly as shown in @xcite , which , however is not the true equilibrium configuration . 
   in addition , as the accessible length scale strongly depends on the system size ( typically less than 1/3 of the side length of the simulation box ) , smaller systems apparently equilibrate faster . 
   we have found @xmath47 to be sufficient to reach equilibrium for a system of @xmath37 particles , while @xmath48 was needed for @xmath41 . 
   to verify , that the observed slowdown of relaxation is not an artifact of the applied microcanonical ( constant @xmath49 ) simulation , we have implemented the computationally much more demanding , but in principle for phase transition studies better suited isothermal - isobaric ( constant @xmath50 ) molecular dynamics scheme @xcite . 
   although the @xmath50 simulations were performed for much smaller systems ( @xmath51 ) , limiting the calculation of the correlation functions to a shorter range and resulting in higher noise levels , the same long - time tendency of decaying long - range correlations could be identified as already shown with the computationally much more efficient microcanonical simulations . 
   during the equilibration of an interacting charged many - particle system we have identified three different stages of relaxation :    * the velocity distribution does approach the maxwellian distribution within a few plasma oscillation cycles . in the close vicinity of the melting transition the speed of this process is found to be independent of temperature and system size . * 
   compared to the velocity distribution function , the configurational temperature ( determined by the local neighborhood within the range of the inter - particle interaction potential ) relaxes at time scales about an order of magnitude longer for our systems . 
   the relaxation time scale is not sensitive to the system size , but has a strong dependence on the temperature . * 
   the equilibration of the long range correlations is significantly slower compared to the above quantities , and depends strongly on the systems size ( larger systems need longer time to equilibrate ) .    from this study 
   we can conclude , that increasing the system size in particle simulations alone can be insufficient and can result in misleading conclusions , as the length of the equilibration period also plays a crucial role in building up or destroying correlations .    in the vast majority of the earlier numerical studies on charged particle ensembles ( as listed in the introduction ) no simulation time is specified , given to the system to equilibrate before the actual measurement were performed , neither is the method of characterizing the quality of the equilibrium described . 
   based on these results , we suspect , that the rapidly increasing computational resources in the first decade of the 21st century beguiled increasing the system sizes in particle simulations without increasing the length of the simulated time intervals . in the majority of these studies the systems may got stuck in the metastable hexatic phase , instead of settling in the true equilibrium configuration . 
   67ifxundefined [ 1 ] ifx#1 ifnum [ 1 ] # 1firstoftwo secondoftwo ifx [ 1 ] # 1firstoftwo secondoftwo `` `` # 1''''@noop [ 0]secondoftwosanitize@url [ 0 ] 
   + 12$12  & 12#1212_12%12@startlink[1]@endlink[0]@bib@innerbibempty @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop `` , ''  ( ,  ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop _ _ ,  ed . 
   ( ,  ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( )
summary: for two - dimensional many - particle systems first - order , second - order , single step continuous , as well as two - step continuous ( kthny - like ) melting transitions have been found in previous studies . 
   recent computer simulations , using particle numbers in the @xmath0 range , as well as a few experimental studies , tend to support the two - step scenario , where the solid and liquid phases are separated by a third , so called hexatic phase . 
   we have performed molecular dynamics simulations on yukawa ( debye - hckel ) systems at conditions earlier predicted to belong to the hexatic phase . 
   our simulation studies on the time needed for the equilibration of the systems conclude that the hexatic phase is metastable and disappears in the limit of long times . 
   we also show that simply increasing the particle number in particle simulations does not necessarily result in more accurate conclusions regarding the existence of the hexatic phase . 
   the increase of the system size has to be accompanied with the increase of the simulation time to ensure properly thermalized conditions .