Effect of Coulomb Screening Length on Nuclear Pasta Simulations
P. N. Alcain, P. A. Giménez Molinelli, J. I. Nichols, C. O. Dorso
EEffect of Coulomb Screening Length on Nuclear Pasta Simulations
P. N. Alcain, P. A. Gim´enez Molinelli, J. I. Nichols and C. O. Dorso
Departamento de F´ısica, FCEN, Universidad de Buenos Aires, N´u˜nez, Argentina andIFIBA-CONICET (Dated: September 19, 2018)We study the role of the effective Coulomb interaction strength and length on the dynamics ofnucleons in conditions according to those in a neutron star’s crust. Calculations were made witha semi-classical molecular dynamics model, studying isospin symmetric matter at sub-saturationdensities and low temperatures. The electrostatic interaction between protons interaction is includedin the form of a screened Coulomb potential in the spirit of the Thomas-Fermi approximation,but the screening length is artificially varied to explore its effect on the formation of the non-homogeneous nuclear structures known as “nuclear pasta”. As the screening length increases, wecan see a transition from a one-per-cell pasta regime (due exclusively to finite size effects) to a moreappealing multiple pasta per simulation box. This shows qualitative difference in the structure ofneutron star matter at low temperatures, and therefore, special caution should be taken when thescreening length is estimated for numerical simulations.
PACS numbers: PACS 24.10.Lx, 02.70.Ns, 26.60.Gj, 21.30.Fe
I. INTRODUCTION
At densities and temperatures expected to exist in neu-tron star crusts ( ρ < ∼ ρ and T < ∼ . M eV , with ρ denot-ing the normal nuclear density), nucleons form structuresthat are substantially different from the “normal”, quasi-spherical nuclei we are familiar with. Such structures,which have been dubbed “nuclear pasta”, have been in-vestigated using various models [1–12] which have shownthem to be the result of the interplay between nuclear andCoulomb forces in an infinite medium. The structure ofthe nuclear pasta is expected to play an important rolein the study of neutrino opacity in neutron stars [13],neutron star quakes and pulsar glitches [14]. In neutronstars, apart from protons and neutrons, there is an elec-tron gas. This electron gas screens the electrostatic longrange proton-proton interaction. This screening effect isoften modeled within the Thomas-Fermi approximation,according to which, the interaction between protons is aYukawa-like potential with a screening length λ : V T F ( r ) = q e − r/λ r According to QFT calculations [15, pp. 175-180], thescreening length is λ ≈
100 fm. For numerical simula-tions, such a long range interaction poses a problem since,to perform correct particle-based simulations, the simu-lation domain (or cell) should be much larger than thelength of the interaction potentials [16]. Using the cor-rect value for λ would then require working with O (10 )particles and it would be computationally very exhaus-tive. Facing this issue, pioneering authors [13, 17] de-cided to work with a much smaller λ = 10 fm, hoping toretain the main qualitative phenomenological aspects ofthe system (competing interactions of different length)but with smaller systems. Even if they were indeed ca-pable of producing “pasta-like” structures, the particular choice of the value for the screening length was arbitraryand based almost solely on computational details. No-tably, this particular value of the screening length wasused by every author using a screened Coulomb poten-tial for particle-based simulations ever since [7, 12, 13].This paper will work on studying to which extent this ar-bitraty choice is physically relevant, expanding previousworks.The role of the length of the screening has been nar-rowly explored by several authors in other models. Forexample, in a 2003 investigation [11], the screening ef-fect of an electron gas on cold nuclear structures wasinvestigated using a static liquid-drop model, and it wasfound that main effect of the gas screening was to ex-tend the range of densities where bubbles and clustersappear and to reduce the range of stability of homoge-neous phases. While the screening was found to be ofminor importance, the study, being static, didn’t includeany spatial or dynamical effect. Another 2005 study [18]used a density functional method to investigate chargescreening on nuclear structures at sub-nuclear densitiesbut still at zero temperature; in particular, cases withand without screening were directly compared. The mainresults of the study were nucleon density profiles usedto quantify the spatial rearrangement of the proton andelectron charge densities. Once again, it was found thatthe density region in which the pasta exists becomesbroader when the Coulomb screening is taken into ac-count, mainly due to the rearrangement of the protons;the authors remark the importance of extending suchstudy to finite temperatures and with dynamical mod-els.More recently, some works [19, 20] began studying theeffect that Coulomb interaction has on the pasta for-mation using dynamical models. The main findings arethat artificial one-per-cell pasta ( pseudo-pasta ) could ex-ist even when Coulomb interaction was removed, andthat they exist due to periodic boundary conditions and a r X i v : . [ nu c l - t h ] N ov finite size [21].In a previous study, we showed that a combination ofclassical molecular dynamics, fragment recognition algo-rithms and a set of topological tools [12] was very effec-tive in the study of the pasta structures. In particular, weshowed that topological observables can be used to clas-sify the different structures into recognizable patterns;this allows for cross-model comparison between struc-tures obtained with different approaches and to a quan-tifiable analysis of the effect the nuclear and Coulombenergy have on the pasta formation and properties. Wewill study the extent to which the screening length λ affects the morphology of the ground states at zero tem-perature of nuclear pasta. To this effect, semi-classicalmolecular dynamic simulations with a screened Coulombpotential and values of λ from 0 to 50 fm were performed.The model used is described in section II, and in sec-tion II A numerical aspects of the Coulomb model usedare discussed. Topological tools were used to quantita-tively analyse the effect the nuclear and Coulomb ener-gies have on the pasta formation and its properties areintroduced in section II C. Results are presented and dis-cussed in section III. II. CLASSICAL MOLECULAR DYNAMICSMODEL
The model used here was and developed to study nu-clear reactions from a semi-classical, particle-based pointof view [22]. The justification for using this model in stel-lar crust environments was presented elsewhere [12], herewe simply mention some basic ingredients of the model.The classical molecular dynamics model
CM D , as in-troduced in [23], is retrofitted with cluster recognition al-gorithms and a plethora of analysis tools. It has been suc-cessfully used in heavy-ion reaction studies to help under-stand experimental data [24], identify phase transitionssignals and other critical phenomena [25–28], explore thecaloric curve [29, 30] and isoscaling [31, 32]. Synopti-cally,
CM D uses two two-body potentials to describe themotion of nucleons by solving their classical equations ofmotion. The potentials,developed phenomenologically byPandharipande [22], are: V np ( r ) = v r exp( − µ r r ) /r − v a exp( − µ a r ) /rV nn ( r ) = v exp( − µ r ) /r where V np is the potential between a neutron and a pro-ton, and V nn is the repulsive interaction between either nn or pp . The cutoff radius is r c = 5 . r > r c both potentials are set to zero. The Yukawaparameters µ r , µ a and µ were determined to yield anequilibrium density of ρ = 0 .
16 fm − , a binding en-ergy E ( ρ ) = 16 MeV/nucleon and a compressibility of250 MeV [22].The main advantage of the CM D model is the pos-sibility of knowing the position and momentum of all particles at all times. This allows the study of the struc-ture of the nuclear medium from a particle-wise point ofview. The output of
CM D , namely, the time evolution ofthe particles in ( r , p ), can be used as input in anyone ofthe several cluster recognition algorithms that some of ushave designed for the study of nuclear reactions [33–35].As explained elsewhere [12, 36, 37], the lack of quan-tum effects, such as Pauli blocking, –perhaps the onlyserious caveat in classical models– ceases to be rele-vant in conditions of high density and temperature (suchas in heavy-ion reactions) or in the low-density andlow-temperature stellar environments, when momentumtrasfer between particles ceases to be important.To simulate an infinite medium, systems with thou-sands of nucleons were constructed using CM D under pe-riodic boundary conditions. Cases symmetric in isospin(i.e. with x = Z/A = 0 .
5, 2500 protons and 2500 neu-trons) were constructed in cubical boxes with sizes ad-justed to have densities between ρ = 0 . f m − ≤ ρ ≤ .
08. Although in the actual neutron stars the protonfraction is low ( x < . A. Coulomb interaction in the Model
To take into account the Coulomb interaction, which isformally of infinite range, in molecular dynamics simula-tion under periodic boundary conditions, it is necessaryto use some approximation. The two most common ap-proaches are the Thomas-Fermi screened Coulomb poten-tial (used with various nuclear models, e.g., in
QM D [7],
CM D [12] and
SSP [13]) and the Ewald summationprocedure [38]. Theoretical estimations for the screen-ing length λ are λ ∼
100 fm, but in the previously men-tioned works, due to computational limitations, a valueof λ = 10 fm was chosen. Our goal on this work is tounderstand the effect this a priori arbitrary choice hason the properties of the ground states of Neutron StarMatter within the framework of CMD.In this work, we used values of λ rangingfrom λ = 0 fm (formally, no Coulomb inter-action) and λ = 20 fm, for densities ρ = { .
005 fm − , .
03 fm − , .
05 fm − , .
08 fm − } , and thecut-off length was chosen at r c = λ . In particular for ρ = 0 .
005 fm − , where “gnocchi” are formed, we ex-tended the analysis to λ = 30 fm and λ = 50 fm to per-form a quantitative analysis on the physical properties ofthe clusters. B. Simulation procedure
The trajectories of the nucleons are then governed bythe Pandharipande and the screened Coulomb potentials.The nuclear system is cooled from T = 1 . T =0 .
001 MeV using isothermal molecular dynamics with theNos´e-Hoover thermostat procedure [40], in the LAMMPSpackage [41]. Systems are cooled in small temperature s(∆ T ≈ . C. Analysis tools
The first of the analysis tools used is the pair corre-lation function, g ( r ), which gives information about thespatial ordering of the nuclear medium. In the previousstudy, g ( r ) showed that nucleons in clusters have an inter-particle distance of about 1 . λ = 10 fm. It is interesting to know if the nearestneighbor distance changes with the screening length.Beyond local measures, the shapes of nuclear struc-tures can be characterized by a set of morphological andtopological observables: their volume, surface area, meancurvature, and Euler characteristic χ . These four objectscomprise the “Minkowski functionals” and completelydescribe all morphological and topological properties ofany three-dimensional object [39]. The computation ofthe mean curvature and χ can be accomplished throughthe Michielsen–De Raed algorithm but requires the map-ping of the nuclear clusters into a polyhedra, proceduredescribed in [12].In [12] it was shown that generic structures, suchas “gnocchi”, “spaghetti”, “lasagna” and “crossed-lasagnas” or “jungle gym” and their inverse structures(with voids replacing particles and vice-versa), all havewell defined and distinct values of the mean curvatureand χ with magnitudes dictated by the overall size ofthe structure, i.e. by the number of particles used. Forstructures with near-zero Euler number, which signalspaghetti, lasagna and their anti-structures, the valuesof the Minkowski functionals are sensitive to the choiceof two parameters: the size given to each particle andthe size of the cells in which we partition the space. Thisclassification is shown in Table I.TABLE I: Classification Curvature - Euler Curvature < ∼ > > ∼ < III. RESULTS AND DISCUSSION
As observed in previous works [19, 20], in absence ofany Coulomb interaction (what would be equivalent to λ = 0), pasta-like structures can be seen, although onlyone per cell. These pseudo-pastas are also shaped inspheres, rods, slabs, anti-rods and anti-spheres, just likethe pasta with Coulomb interaction. The main difference is that, without the Coulomb interaction, we find alwaysone structure per cell, giving the hint that its structureis related to the periodic boundary condition imposed onthe box. The pseudo-pasta exists due to finite size effectsand, if the box was not to exist, the solution would be aninfinite droplet. We notice, however, that when there isCoulomb interaction, the competition between opposinginteractions gives rise to a characteristic length. At sub-saturation densities, this competition is responsible forthe pasta phases, and in the limit of very-low densities(and no screening) it shapes the nuclei we are used to.By increasing the value of λ from 0 to 20 fm we aim toexplore the transition from artificial one-per-cell pasta,to more realistic situations with more than one structureper cell. Moreover, this enables us to assess the physicalimplications of the arbitrary and traditional λ = 10 fmvalue. A. “Critical” Screening Length
A first approach to analyze the nature of the pastaobtained is given by taking a quick glance as the pressureof the different configurations.The pressure is computed by the virial formula P = N k B TV + 13 Σ Ni r i · F i V where N is the number of nucleons in the system and F the force exerted upon each nucleon. The terms in thevirial formula applie only to the interactions specific tothe model, not contemplating the electon gas pressure.This pressure is not to be mistaken with the pressure ex-pected in neutron star crusts, since electrons should beconsidered explicitly in order to calculate it correctly, itis merely a test of the mechanical stability of the config-urations obtained within this model. In figure 1, we seethat for all λ <
10 fm the pressure is negative.The negative pressure is a signal that the non-homogeneous structures found are artificial and that thestructures found can only exist under periodic bound-ary conditions (see [19, 21]). This may be better under-stood by picturing the primitive cell of the simulationas being under the stress caused by its periodic repli-cas. This mean that for such small screening lengths theoverall effective interaction is mostly attractive and pe-riodic boundary conditions are still playing a major rolein shaping the ground state.For λ >
10 fm the pressure becomes positive, meaningthat the structures formed in these configurations are notonly due to periodic boundary conditions, but Coulombinteraction is beginning to play its intended role. Theconfigurations for these values of λ indeed show den-sity fluctuations of length smaller than the size of thecell, which can only be attributed to the coulomb-nuclearcompetition. However, the morphology of the structures,as characterized by the topological measures described insection II C, changes drastically with λ .FIG. 1: Pressure as a function of λ for differentdensities. We see that for λ <
10 fm, the pressure isnegative, implying that periodic boundary conditionsare affecting the morphology of the solution.In order to classify the low temperature ( T =0 .
001 MeV) structures for each value of λ , we study theirmorphology with the analysis tools for the spatial dis-tribution of the particles; namely, the Minkowski func-tionals. In figures 2 we can see the surface, curvatureand Euler number for the ground states, and their de-pendence on λ for different densities.As stated in table I, we expect lasagna and spaghetti to have an Euler number χ = 0. In the gnocchi case,however, each one of them contributes with χ gn = 2.This means that the Euler number of the whole systemwill be χ = 2 · N gn . As the configurations break up intomultiple structures per cell with increasing λ , we expectthe surface to increase as well. As for the curvature, thebehavior described in table I (positive for spaghetti and gnocchi , zero for lasagna and negative for tunnel ) is onlyobserved for λ = 20 fm. Between λ = 7 fm and λ = 10 fmall three of the Minkowski functionals change drasticallybefore reaching well defined values. This indicates thatthere is a transition regime where the structures cannotbe described as any of the traditional pasta. For thismodel of nuclear interaction and Coulomb treatment, itseems the usual λ = 10 fm value is actually too small. B. One vs Many
To better understand how the ground state at low tem-peratures differ on the transition regime from withoutCoulomb interaction to the λ = 20 fm screening length,we show in figure 3 visual representations of the resultsobtained at a set of chosen densities, both with λ = 0and λ = 20 fm.As an example, we show the pair distribution function, g ( r ), for ρ = 0 .
05 fm − in figure 4. In it we see thatthe first peaks of the distribution remain on the samedistances of r = 1 . , . SurfaceCurvatureEuler number
FIG. 2: Minkowski functionals dependence with λ . Wecan see that there is a transition regime between λ = 7 fm and λ = 15 fm, where the Minkowskifunctionals are changing.range structure is governed by the nuclear potential evenat λ = 20 fm, which is evident simply by comparing theorders of magnitude of V n − n and V Coulomb at such short ρ = 0 .
03 fm − , λ = 0 fm ρ = 0 .
03 fm − , λ = 20 fm ρ = 0 .
05 fm − , λ = 0 fm ρ = 0 .
05 fm − , λ = 20 fm ρ = 0 .
08 fm − , λ = 0 fm ρ = 0 .
08 fm − , λ = 20 fm FIG. 3: Difference between pasta with and withoutCoulomb interaction. We can see that the Coulombinteraction splits up the pasta, converting one structureper cell to multiple structures per cell.ranges.When increasing from λ = 15 fm to λ = 20 fm at ρ = 0 .
005 fm − , although qualitatively we see the samebehavior (both show gnocchi ), the average size of clusteris different for these two values of screening length. Thisimplies that the number of clusters is different, hencethe difference observed in the Minkowski functionals. Tostudy this result further, we plot the “gnocchi” size asa function of λ in figure 5. We see that, when we con-sider the standard deviation in the mass distribution, itremains unchanged for λ ≥
20 fm. The average relative FIG. 4: Examples of the radial correlation function for ρ = 0 .
05 fm − and two screening lengths: λ = 20 fm(top panel) and λ = 0 fm (lower panel). Please noticethe difference in the y-scales of the graphs.error on this graph is e ≈ λ (cid:54) = 0) the original λ = 0 pseudo-pasta splits up:from one structure per cell to multiple structures per cell.For intermediate to low values of λ <
20 fm, the effect ofthe periodic boundary conditions is still observable forsome densities and more exotic structures which can beconfused with “true” pasta may exist.
C. The Transition Regime
We now turn to analyze the structures found in thetransition regime.We take, as an example, the lowest density ( ρ =0 .
005 fm − ). As can be seen in figure 6, for λ = 0only one droplet is formed, as expected. For λ = 10 fm,we can see that many gnocchi exist, but some of themstick to their neighbors forming lumps of different sizes.Although Coulomb interaction is now strong enough tobreak the one-pasta found with λ = 0 into many, the re-sulting droplets are not fully fledged gnocchi that can bearranged in a regular lattice such as those for λ = 20 fm. λ = 0 fm λ = 10 fm λ = 20 fm FIG. 6: Different structures got while varying the λ parameter, for ρ = 0 . − . In the transition regime,we find, at λ = 10 fm, that the structure breaks down tomany short-spaghetti -like parts. IV. DISCUSSION AND CONCLUDINGREMARKS
The effect of the screening length of the Coulomb inter-action in simulations of neutron star matter was studiedat densities comparable to that of neutron star crusts.Subsequently along the literature we can find that thevalue of the screening length in the Thomas-Fermi ap-proximation is λ ≈
100 fm. For particle-based simula-tions, due to computational limitations this value washistorically and arbitrarily reduced to λ ≈
10 fm . Thiswas done expecting to mantain the basic phenomenol-ogy when simulating small sytems. We found, though,that there is a critical screening length λ c at which thestructure of the ground state drastically changes. ForPandharipande potential it lies between 10 fm and 15 fm(depending on the density). For λ < λ c , the Coulombinteraction is barely acting and the non-homogeneousstructures emerging from the simulations are due to finitesize effects, as made evident from the negatve pressure ofsuch structures and the fact that there is only one struc-ture per cell. For λ > λ c the pressure becomes positiveand the systems present density fluctuations at a scalesmaller than that of the cell, but not well shaped. Thistransition regime is characterized by large fluctuations inthe surface, curvature and Euler characteristic χ of thestructures. It is only for λ = 20 fm that the morphologyof the structures formed stabilize and cease to deppendon λ . Moreover, the structures in this regime are theusual pasta phases.Because of this, we believe extreme caution shouldbe taken when choosing an arbitrary value for λ , sinceeven though some results at λ = 10 fm can look like theexpected pasta, the results obtained for that particularchoice of λ may be quite different from those in the trueThomas-Fermi approximation. In conclusion, we find thechoice of a good value for λ that is computationally man-ageable and can still adequatelly recover the physics ofthe Thomas-Fermi approximation is no trivial task, anda rigorous study needs to be done prior to the choice ofthe value. A good value for λ must lie in the λ > λ c re-gion which is bound to be dependent on the model usedfor the nuclear interaction. ACKNOWLEDGMENTS
C.O.D. is supported by CONICET Grant PIP0871,P.G.M., J.I.N and P.N.A. by a CONICET scholarship.The three-dimensional figures were prepared using thesoftware
Visual Molecular Dynamics [42]. [1] D. G. Ravenhall, C. J. Pethick and J. R. Wilson, Phys.Rev. Lett. , 2066 (1983). [2] M. Hashimoto, H. Seki and M. Yamada, Prog. Theor. Phys. , 320 (1984).[3] R. D. Williams and S. E. Koonin, Nucl. Phys. A435 , 844(1985).[4] K. Oyamatsu, Nucl. Phys.
A561 , 431 (1993).[5] C. P. Lorenz, D. G. Ravenhall and C. J. Pethick,Phys.Rev. Lett. , 379 (1993).[6] K. S. Cheng, C. C. Yao and Z. G. Dai, Phys. Rev. C55 ,2092 (1997).[7] T. Maruyama, K. Niita, K. Oyamatsu, T. Maruyama, S.Chiba and A. Iwamoto, Phys. Rev.
C57 , 655 (1998).[8] T. Kido, Toshiki Maruyama, K. Niita and S. Chiba, Nucl.Phys.
A663-664 , 877 (2000).[9] G. Watanabe, K. Iida and K. Sato, Nucl. Phys.
A676 ,445 (2000).[10] G. Watanabe, K. Sato, K. Yasuoka and T. Ebisuzaki,Phys. Rev.
C66 , 012801 (2002).[11] G. Watanabe and K. Iida, Phys. Rev.
C68 , 045801(2003).[12] C.O. Dorso, P.A. Gim´enez Molinelli and J.A. L´opez,in “Neutron Star Crust”, Eds. C.A. Bertulani and J.Piekarewicz, in print, Nova Science Publishers, Inc.(2012).[13] C.J. Horowitz, M.A. Perez-Garcia, and J. Piekarewicz,Phys. Rev.
C69 , 045804 (2004).[14] Y. Mochizuki and T. Izuyama, Astrophys. J. , 263(1995).[15] A. L. Fetter and J. D. Wallecka,
Quantum Theory ofMany-Particle Systems , McGraw-Hill (1971)[16] D. Frenkel and B. Smit, “Understanding Molecular Sim-ulations”, 2nd Ed.,Academic Press (2002).[17] Maruyama, Watanabe and Chiba, Prog. Theor. Exp.Phys. , 1, (2012)[18] T. Maruyama, T. Tatsumi, D.N. Voskresensky, T. Tani-gawa and S. Chiba, Phys. Rev.
C72 , 015802 (2005).[19] C. O. Dorso, P. A. Gim´enez Molinelli, J. I. Nichols, J. A.L´opez, http://arxiv.org/abs/1211.5582[20] A. S. Schneider, C. J. Horowitz, J. Hughto, D. K. Berry,http://arxiv.org/abs/1307.1678[21] Binder, K. et al, Am. J. Phys. , 1099[22] A. Vicentini, G. Jacucci and V. R. Pandharipande, Phys.Rev. C31 , 1783 (1985); R. J. Lenk and V. R. Pand-haripande, Phys. Rev.
C34 , 177 (1986); R.J. Lenk, T.J.Schlagel and V. R. Pandharipande, Phys. Rev.
C42 , 372(1990). [23] A. Barra˜n´on, C.O. Dorso, J.A. L´opez and J. Morales,Rev. Mex. F´ıs. , 110 (1999).[24] A. Chernomoretz, L. Gingras, Y. Larochelle, L. Beaulieu,R. Roy, C. St-Pierre and C. O. Dorso, Phys. Rev. C65 ,054613 (2002).[25] A. Barra˜n´on, C.O. Dorso and J.A. L´opez, Rev. Mex. F´ıs. , 93 (2001).[26] A. Barra˜n´on, C.O. Dorso, and J.A. L´opez, Nuclear Phys.
A791 , 222 (2007).[27] A. Barra˜n´on, R. C´ardenas, C.O. Dorso, and J.A. L´opez,Heavy Ion Phys. , 59 (2003).[28] C.O. Dorso and J.A. L´opez, Phys. Rev. C64 , 027602(2001).[29] A. Barra˜n´on, J. Escamilla Roa and J.A. L´opez, Braz. J.Phy., , 904 (2004).[30] A. Barra˜n´on, J. Escamilla Roa and J.A. L´opez, Phys.Rev.
C69 , 014601 (2004).[31] C.O. Dorso, C.R. Escudero, M. Ison and J.A. L´opez,Phys. Rev.
C73 , 044601 (2006).[32] C.O. Dorso, P.A. Gim´enez Molinelli and J.A. L´opez, J.Phys. G: Nucl. Part. Phys. , 115101 (2011); ibid , Rev.Mex. Phys., S 57 (1) , 14 (2011).[33] C.O. Dorso and J. Aichelin, Phys. Lett.
B345 , 197(1995).[34] A. Strachan and C.O. Dorso, Phys. Rev.
C55 , 775(1997); ibid , Phys. Rev.
C56 , 995 (1997).[35] C.O. Dorso and J. Randrup, Phys. Lett.
B301 , 328(1993).[36] J.A. L´opez and C.O. Dorso, Lecture Notes on PhaseTransformations in Nuclear Matter, World Scientific,Hackensack, NJ, USA, ISBN 978-981-02-4007-3 (2000).[37] J. Taruna, “The physics of compact stars”, Ph.D. Thesis,Florida State University, AAT 3321532 (2008).[38] G. Watanabe, K. Sato, K. Yasuoka and T. Ebisuzaki,Phys. Rev.
C68 , 035806 (2003)[39] K. Michielsen and H. De Raedt, Phys. Reports , 461(2001).[40] S. Nose, J. Chem. Phys. , 511 (1984)[41] S. Plimpton, J. Comp. Phys., , 1-19 (1995)[42] Humphrey, W., Dalke, A. and Schulten, K., ‘VMD - Vi-sual Molecular Dynamics’, J. Molec. Graphics14.1