text:
# Untitled Document

Evolution of the polarization of the optical afterglow of the gamma-ray burst GRB 030329

Jochen Greiner, Sylvio Klose, Klaus Reinsch, Hans Martin Schmid, Re’em Sari$, Dieter H. Hartmann, Chryssa Kouveliotou§, Arne Rau, Eliana Palazzi, Christian Straubmeier∗∗, Bringfried Stecklum, Sergej Zharikov††, Gaghik Tovmassian††, Otto Bärnbantner, Christoph Ries, Emmanuel Jehin¶¶, Arne Henden§§, Anlaug A. Kaas‡‡, Tommy Grav$$, Jens Hjorthℓℓ, Holger Pedersenℓℓ, Ralph A.M.J. Wijers#, A. Kaufer¶¶, Hye-Sook Park♢♢, Grant Williams♯♯, Olaf Reimer##

Max-Planck-Institut für extraterrestrische Physik, 85741 Garching, Germany

Thüringer Landessternwarte, 07778 Tautenburg, Germany

Universitäts-Sternwarte Göttingen, 37083 Göttingen, Germany

Institut für Astronomie, ETH Zürich, 8092 Zürich, Switzerland

$California Institute of Technology, Theoretical Astrophysics 130-33, Pasadena, CA 91125, USA

Clemson University, Department of Physics and Astronomy, Clemson, SC 29634, USA

§NSSTC, SD-50, 320 Sparkman Drive, Huntsville, AL 35805, USA

Istituto di Astrofisica Spaziale e Fisica Cosmica, CNR, Sezione di Bologna, 40129 Bologna, Italy

∗∗ I. Physikalisches Institut, Universität Köln, 50937 Köln, Germany

††Instituto de Astronomia, UNAM, 22860 Ensenada, Mexico

Wendelstein-Observatorium, Universitätssternwarte, 81679 München, Germany

¶¶European Southern Observatory, Alonso de Cordova 3107, Vitacura, Casilla 19001, Santiago 19, Chile

§§Universities Space Research Association, U.S. Naval Observatory, P.O. Box 1149, Flagstaff, AZ 86002, USA

‡‡Nordic Optical Telescope, 38700 Santa Cruz de La Palma, Spain

$$University of Oslo, Institute for Theoretical Astrophysics, 0315 Oslo, Norway, and Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA

ℓℓAstronomical Observatory, NBIfAFG, University of Copenhagen, 2100 Copenhagen , Denmark

#Astronomical Institute Anton Pannekoek, Kruislaan 403, 1098 SJ Amsterdam, The Netherlands

♢♢Lawrence Livermore National Laboratory, University of California, P.O. Box 808, Livermore, CA 94551, USA

♯♯MMT Observatory, University of Arizona, Tucson, AZ 85721, USA

##Theoretische Weltraum-und Astrophysik, Ruhr-Universität Bochum, 44780 Bochum, Germany

The association of a supernova with GRB 0303291,2 strongly supports the collapsar model3 of $`\gamma `$-ray bursts (GRBs), where a relativistic jet4 forms after the progenitor star collapses. Such jets cannot be spatially resolved because of their cosmological distances. Their existence is conjectured based on breaks in GRB afterglow light curves and the theoretical desire to reduce the GRB energy requirements. Temporal evolution of polarization5,6,7 may provide independent evidence for the jet structure of the relativistic outflow. Small-level polarization ($``$1-3%)8-17 has been reported for a few bursts, but the temporal evolution of polarization properties could not be established. Here, we report polarimetric observations of the afterglow of GRB 030329 with high signal-to-noise and high sampling frequency. We establish the polarization light curve, detect sustained polarization at the percent level, and find significant variability. The data imply that the afterglow magnetic field has small coherence length and is mostly random, probably generated by turbulence, in contrast with the high polarization detected in the prompt $`\gamma `$-rays from GRB 02120618. Our results suggest a different structure and origin of the magnetic field in the prompt vs. afterglow emission regions.

GRB 030329 triggered the High Energy Transient Explorer, HETE-II, on March 29, 2003 (11:37:14.67 UT)19. The discovery of the burst optical afterglow20,21 was quickly followed by a redshift measurement22 for the burster of z=0.1685 ($``$800 Mpc), thus making GRB 030329 the second closest long-duration23 GRB ever studied, after GRB 9804253. The proximity of GRB 030329 resulted in very bright prompt and afterglow emission, leading to the best-sampled afterglow to date. Detailed optical spectroscopy revealed an underlying supernova (SN 2003dh)1,2 with an astonishing spectral similarity to SN 1998bw (at z=0.0085, associated with GRB 9804253), thus strongly supporting the link of long-duration GRBs with core-collapse supernovae.

The apparent brightness of GRB 030329 also afforded a unique opportunity to study the temporal evolution of polarization in the afterglow phase. Previous single measurements found low-level (1-3 %) polarization8-15 in optical afterglows, and the only reports on variable polarization16,17 were based on few measurements with different instruments and modest signal-to-noise. We have overcome the previous sampling limitations with 31 polarimetric observations of the afterglow of GRB 030329 obtained with the same instrumentation (plus few more with different instruments) over a time period of 38 days for a total of an unprecedented $``$50 hours of observations with an 8 m telescope. For the first time we establish a polarization light curve for an optical gamma-ray burst afterglow.

We performed relative photometry, and derived from each pair of simultaneous measurements at orthogonal angles the Stokes parameters U and Q. In order to obtain the intrinsic polarization of the GRB afterglow, we have to correct for Galactic interstellar polarization (mostly due to dust). We performed imaging polarimetry to derive the polarization parameters of seven stars in the field of GRB 030329, and obtained an interstellar (dust) polarization correction of 0.45% at position angle 155$`\mathrm{deg}`$. Subtraction of the mean foreground polarization was performed in the Q/U plane (Qfp=0.0027$`\pm `$0.0013, Ufp=-0.0033$`\pm `$0.0017). Figure 1 shows an R band image of the GRB 030329 field with the polarization ”vectors” superimposed.

The temporal evolution of the degree and angle of polarization together with the R band photometry is shown in Figure 2. This figure demonstrates the presence of non-zero polarization, $`\mathrm{\Pi }`$0.3-2.5 % throughout a 38-day period, with significant variability in degree and angle on time scales down to hours. Further, the spectropolarimetric data of the first three nights as well as the simultaneous R and K band imaging polarimetry during the second night show that the relative polarization and the position angle are wavelength independent (within the measurement errors of about 0.1 %) over the entire spectral range. These data imply that polarization due to dust in the host galaxy of GRB 030329 does not exceed $``$0.3 %. Dust destruction in the vicinity of the gamma-ray burster due to the strong radiation field would result in a monotonically decreasing polarization degree at constant position angle, which is not observed during the first days. We therefore conclude that the bulk of the observed variability is intrinsic to the afterglow.

Figure 2 shows that while the polarization properties show substantial variability (for which no simple empirical relationship is apparent), the R band flux is a sequence of power laws. During each of the power law decay phases the polarization is of order few percent, different from phase to phase, and variable within the phase, but not in tandem with the ”bumps and wiggles” in the light curve. We observe a decreasing polarization degree shortly after the light curve break at $``$0.4 days (as determined from optical21,24 and X-ray data25,26). Rapid variations of polarization occur $``$1.5 days after the burst, and could be related to the end of the transition period towards a new power law phase starting at $``$1.7 days. Polarization eventually rises to a level of $``$2 %, which remains roughly constant for another two weeks. At late time the underlying Type Ic SN 2003dh1,2 increasingly contributes to the total light and probably also to the observed polarization properties. Asymmetries in this type of supernova can produce polarization of order 1 %27,28, as we observe towards the end of our campaign.

GRB 030329 belongs to a growing group of bursts for which densely monitored afterglow light curves show significant ”bumps” and ”wiggles” relative to a simple power law decay (most notably GRB 021004 and GRB 011211). These bumps and wiggles complicate the interpretation of the polarization properties. While the rapid decrease in polarization degree during the first night is consistent with the model predictions5,6,7, the position angle changes are not, and thus it remains to be proven whether or not these early polarization data support the break at 0.4 days as due to a jet. A connection to theoretical models throughout the entire period covered in Fig. 2 is a daunting task. A detailed comparison of our data with characteristic features of various theoretical models is beyond the scope of this paper, and will be presented elsewhere.

In summary, our data constitute the most complete and dense sampling of the polarization behaviour of a GRB afterglow to date. For GRB 030329 we conclude that the afterglow polarization probably did not rise above $``$2.5 % at any time, and that the polarization did not correlate with the flux. The low level of polarization implies that the components of the magnetic field parallel and perpendicular to the shock do not differ by more than $``$10 %, and suggests an entangled magnetic field, probably amplified by turbulence behind shocks, rather than a pre-existing field. This is in contrast with the high level of polarization detected in the prompt $`\gamma `$-rays from GRB 02120618 and suggests a different structure and origin of the magnetic field in the prompt vs. afterglow emission regions. Evolving polarization properties provide a unique diagnostic tool for GRB studies, and the extremely complex light curve of the optical afterglow of GRB 030329 emphasises that measurements should be carried out with high sampling frequency.

1. Stanek, K.Z., et al. Spectroscopic discovery of the supernova 2003dh associated with GRB 030329. Astrophys. J. 591, L17-L20 (2003)

2. Hjorth, J., et al. A very energetic supernova associated with the $`\gamma `$-ray burst of 29 March 2003. Nature 423, 847-850 (2003)

3. Woosley, S.E., Eastman, R.G., Schmidt, B.P. Gamma-ray bursts and type Ic supernova SN 1998bw. Astrophys. J. 516, 788-796 (1999)

4. M sz ros, P. Theories of gamma-ray bursts. Ann. Rev. Astron. Astrophys. 40, 137-169 (2002)

5. Sari, R. Linear polarization and proper motion in the afterglow of beamed gamma-ray bursts. Astrophys. J. 524, L43-L46 (1999)

6. Gruzinov, A. Strongly polarized optical afterglows of gamma-ray bursts. Astrophys. J. 525, L29-L31 (1999)

7. Ghisellini, G., Lazzati, D. Polarization light curves and position angle variation of beamed gamma-ray bursts. Mon. Not. R. Astron. Soc. 309, L7-L11 (1999)

8. Hjorth, J., et al. Polarimetric constraints on the optical afterglow emission from GRB 990123. Science 283, 2073-2075 (1999)

9. Wijers, R.A.M.J., et al. Detection of polarization in the afterglow of GRB 990510 with the ESO Very Large Telescope. Astrophys. J. 523, L33-L36 (1999)

10. Covino, S., et al. GRB 990510: linearly polarized radiation from a fireball. Astron. Astrophys. 348, L1-L4 (1999)

11. Rol, E., et al. GRB 990712: First indication of polarization variability in a gamma-ray burst afterglow. Astrophys. J. 544, 707-711 (2000)

12. Björnsson, G., Hjorth, J., Pedersen, K., Fynbo, J.U. The afterglow of GRB 010222: A case of continuous energy injection. Astrophys. J. 579, L59-L62 (2002)

13. Masetti, N., et al. Optical and near-infrared observations of the GRB 020405 afterglow. Astron. Astrophys. 404, 465-481 (2003)

14. Covino, S., et al. Polarization evolution of the GRB 020405 afterglow. Astron. Astrophys. 400, L9-L12 (2003)

15. Bersier, D., et al. The strongly polarized afterglow of GRB 020405. Astrophys. J. 583, L63-L66 (2003)

16. Barth, A., et al. Optical spectropolarimetry of the GRB 020813. Astrophys. J. 584, L47-L51 (2003)

17. Rol, E., et al. Variable polarization in the optical afterglow of GRB 021004. Astron. Astrophys. 405, L23-L27 (2003)

18. Coburn, W., and Boggs, S.E. Polarization of the prompt (-ray emission from the (-ray burst of 6 December 2002. Nature 423, 415-417 (2003)

19. Vanderspek, R., Crew, G., Doty, J., Villasenor, J., Monelly, G. GRB 030329 (=H2652): A long, extremely bright GRB localized by HETE WXM and SCX. GCN Circ. 1997 (2003)

20. Price, P.A., et al. The bright optical afterglow of the nearby (-ray burst of 29 March 2003. Nature 423, 844-847 (2003)

21. Uemura, M., et al. Structure in the early afterglow light curve of the (-ray burst of 29 March 2003. Nature 423, 843-844 (2003)

22. Greiner, J., et al. Redshift of GRB 030329. GCN Circ. 2020 (2003)

23. Kouveliotou, C., et al. Identification of two classes of gamma-ray bursts. Astrophys. J. 413, L101-L104 (1993)

24. Burenin, R., et al. First hours of the GRB 030329 optical afterglow. Astron. Lett. 29, 9, 1-6 (2003)

25. Marshall, F.E., Markwardt, C., Swank, J.H. GRB 030329: fading X-ray afterglow with RXTE. GCN Circ. 2052 (2003)

26. Tiengo, A., et al. The X-ray afterglow of GRB 030329. Astron. Astrophys., in press; astro-ph/0305564 (2003)

27. Wang, L., Howell, D. A., H flich, P., and Wheeler, J. C. Bipolar supernova explosions. Astrophys. J. 550, 1030-1035 (2001)

28. Leonard, D. C., Filippenko, A. V., Chornock, R., and Foley, R. J. Photospheric-phase spectropolarimetry and nebular-phase spectroscopy of the peculiar type Ic supernova 2002ap. Publ. Astron. Soc. Pac. 114, 1333-1348 (2002)

Acknowledgements

This work is primarily based on observations collected at ESO, Chile, with additional data obtained at the German-Spanish Astronomical Centre Calar Alto, operated by the Max-Planck-Institute for Astronomy, Heidelberg, jointly with the Spanish National Commission for Astronomy, the NOT on La Palma, Canary Islands, and the Observatorio Astronomico National, San Pedro, Mexico. We are grateful to the staff at the Paranal, Calar Alto and NOT observatories, in particular A. Aguirre, M. Alises, S. Hubrig, A.O. Jaunsen, C. Ledoux, S. Pedraz, T. Szeifert, L. Vanzi and P. Vreeswijk for obtaining the service mode data reported here.

Figure 1: R band image of the field centred on the optical afterglow of GRB 030329. The 27 FORS1 (Focal Reducer and Low Dispersion Spectrograph at VLT/Antu) measurements are shown as ”vectors” where the length is a measure of the polarization degree, and the orientation indicates the position angle. While the afterglow varies in degree as well as angle, the polarization of seven field stars is constant within 0.1% in polarization degree and 1.5 degree in position angle throughout the 38 days. Linear polarization was measured from sets of exposures with different retarder plate position angles. Imaging polarimetry was obtained during the first four nights (when the afterglow was brighter than 17th mag) from sixteen different retarder plate angles, and from eight angles thereafter. Since the FORS1 polarization optics allows determination of the degree of polarization to an accuracy of ¡3x10-4 and of the polarization angle to $``$ 0.2 deg, we consider the above variance of the field stars to represent the systematic error over the 38-day time period. Observations of polarimetric standard stars reproduced their tabulated values within 5%. The FORS1 retarder plate zero point angle of -1.2 degrees was subtracted from the polarization angle. The position angle has a systematic uncertainty of 1.5 deg, which was added in quadrature to the statistical errors.

Figure 2: Evolution of the polarization during the first 38 days. Top and middle panels show the polarization degree in percent and the position angle in degrees. The red data points are from imaging polarimetry with FORS1/VLT. Spectropolarimetry (blue symbols, 700-800 nm range) was performed during the first three nights. The green points were obtained with CAFOS at the 2.2m telescope at Calar Alto on March 30/31. The magenta data point was obtained with AFOSC at the 2.56m NOT telescope on 2003 April 2. The bottom panel shows the residual R band light curve after subtraction of the contribution of a power law t-1.64 describing the undisturbed decay during the time interval 0.5-1.2 days after the GRB (i.e., after the early break at 0.4 days), thus leading to a horizontal curve. The symbols correspond to data obtained through: the literature (black), the 1m USNO telescope at Flagstaff (blue), the OAN Mexico (cyan), and FORS1/VLT (red). Lines indicate phases of power law decay, with the first one from early data21 (not shown). Yellow bars mark re-brightening transitions. Contributions from an underlying supernova (solid curved line) do not become significant until $``$10 days after the GRB.


format: text

text:
# Decoherence from Internal Dynamics in Macroscopic Matter-Wave Interferometry

## I Introduction

Superposition of states lies at the heart of quantum mechanics and gives rise to many of its paradoxes. Not only can a particle go through two paths simultaneously, but the wavefunction of a pair of particles flying apart from each other is also entangled into a non-separable superposition of states. However, such strange phenomena have never been observed in our macroscopic world. It has been an important question why and how quantum weirdness disappears in large systems Zurek91 .

Environment, usually described by a huge number of variables, can destroy coherence among the states of a quantum system. This is decoherence. The environment is watching the path followed by the system, and thus suppressing interference effects and quantum weirdness. In macroscopic systems, such process is so efficient that we see only its final result: the classical world around us Joos00 ; Zurek03 . For truly macroscopic superpositions, decoherence occurs on a very short time-scale that it is almost impossible to observe quantum coherences. However, mesoscopic systems present the possibility of investigating the process of decoherence and the transition from quantum to classical behavior Haroche98 . So far many experiments have been realized to generate mesoscopic superpositions Brune96 ; Monroe96  and to decohere them in a controlled way Myatt00 . Recently considerably large molecules have been used to investigate the decoherence, the transition from quantum to classical. For example, the researchers in Wien have observed interference of de Broglie waves of fullerenes (C60 or C70 molecules) and even bigger ones Arndt99 ; Nariz01 ; Brezger02 ; Hornberger03 ; Hackermueller03 ; Hackermueller04 . In this experiment the fullerenes are quite hot as well as big, which means they contain complicated dynamics of their internal degrees of freedom. The internal thermal energy is almost one order of magnitude larger than the kinetic energy of its center of mass (CM) motion. A question naturally arises: is the complex dynamics of the internal degrees of freedom harmful for the interference of the CM motion?

Usually, decoherence has been generated from the coupling with an outer environment such as other particles or fluctuating electromagnetic fields. However, a macroscopic object generically possesses its own environment in itself, namely the internal dynamics (ID) Kolovsky94a  when only small part of the total system, e.g. its CM, is under consideration. In this paper, we would like to address a question: when and how the ID decoheres interference of a macroscopic object. We also show nontrivial expectation that the weak localization of large molecules in a disordered potential can be destroyed by the decoherence generated from the ID without any external perturbation breaking the time reversal symmetry.

Let us consider a macroscopic object consisting of $`N`$ particles exposed to the external potential $`V_{ex}`$. The total Hamiltonian can be written as

$`H`$ $`=`$ $`{\displaystyle \underset{i}{\overset{N}{}}}{\displaystyle \frac{p_i^2}{2m_i}}+V_{in}(\{x_i\})+{\displaystyle \underset{i}{\overset{N}{}}}V_{ex}(x_i)`$ (1)
$`=`$ $`{\displaystyle \frac{P^2}{2M}}+V_{ex}(X)+{\displaystyle \underset{\alpha }{\overset{N1}{}}}{\displaystyle \frac{\mu _\alpha ^2}{2m_\alpha }}+V_{in}(\{\xi _\alpha \})+\mathrm{\Gamma }(X,\{\xi _\alpha \})`$
$``$ $`H_{CM}+H_{in}+\mathrm{\Gamma },`$
where $`V_{in}`$ is an internal (or confinement) potential. $`P`$ and $`X`$ are a momentum and a coordinate of the CM with a mass $`M(=m_i)`$ respectively, while $`\mu _\alpha `$ and $`\xi _\alpha `$ the same quantities of the internal degrees of freedom with the reduced mass $`m_\alpha `$. $`\mathrm{\Gamma }`$ represents the coupling between the CM and the ID. Since one finds $`V_{in}`$ is determined only from the relative coordinates $`\xi _\alpha `$, i.e. $`V_{in}(\{x_i\})=V_{in}(\{\xi _\alpha \})`$, the coupling term $`\mathrm{\Gamma }`$ depends only on the external potential $`V_{ex}`$. It is easy to show that if $`V_{ex}`$ does not correspond to a simple form such as constant, linear, and harmonic, the non-zero coupling $`\mathrm{\Gamma }`$ always exists Park03 . The CM motion can be entangled with its ID when the anharmonic external potential is applied. We call such a non-trivial external potential “nonlinear” since the corresponding force is nonlinear. It is noted that existence of the external potential has nothing to do with generation of the decoherence.

In a usual two slit experiment, a macroscopic object freely flies to a screen. Therefore, no entanglement between CM and ID arises, neither does the decoherence from ID. However, one can ask what happens if the repulsive potential produced by the slits is considered. For example, the van der Waals interaction between the C70 and the grating was assumed to correctly explain experimental results Brezger02 . Since C70 molecule is too small to see the anharmonic shape of the external potential the coupling between the CM and the ID is hardly expected. One can still insist, however, that not only if the slit wall contains more rapidly varying repulsive potential but also if the size of the object becomes larger, for example by using an insulin, then the coupling term might manifest itself. In this paper we consider such situation that the coupling between CM and ID is non-negligible. It must be mentioned that we ignore all the other sources of decoherences except those originating from the ID for our discussion, but will briefly comment it. We also note that neither we derive any new formula nor estimate any values of this effect. We would like to show basic physical mechanism and possibility of the phenomena.

## II Review of decoherence

First, we briefly review the decoherence from an external environment. Figure 1(a) shows usual two path interferometry, where one of the paths interacts with a “bath” of particles, i.e. an environment. One can write an initial total wavefunction in the following way

$$\mathrm{\Psi }(0)=\frac{1}{\sqrt{2}}\left[\psi _1(0,X)+\psi _2(0,X)\right]\chi (0,\xi ),$$
(2)
where $`\psi _{1(2)}`$ and $`\chi `$ denote a wavefunction of a system moving along the path 1(2) and that of a bath described by variables $`\xi `$, respectively. Under the dynamics, which includes interaction with the bath only along the path 2, the wavefunction becomes

$$\mathrm{\Psi }(t)=\frac{1}{\sqrt{2}}\left[\psi _1(t,X)\chi _1(t,\xi )+\psi _2(t,X)\chi _2(t,\xi )\right],$$
(3)
where $`\chi _{1(2)}`$ is a state of the bath with the system going around the path 1(2). Here, we assume the coupling with the bath is small enough to have no influence on the system, but changes only a state of the bath Stern90 ; Fiete03 .

From the reduced density matrix the interference term is given as

$$\mathrm{Re}\left[\psi _1^{}(t,X)\psi _2(t,X)𝑑\xi \chi _1^{}(t,\xi )\chi _2(t,\xi )\right].$$
(4)
The physical meaning of Eq. (4) is obvious. The first term contains usual information of the interference pattern of the system going through the two paths. The second term corresponds to the visibility which measures the decoherence. Eq. (4) allows one to interpret the reduction of the interference in terms of a reduction in the overlap of the bath states for the two paths. Stern, Aharony and Imry Stern90  argue that one can make the identification

$$e^{i\widehat{\varphi }}=𝑑\xi \chi _1^{}(t,\xi )\chi _2(t,\xi ),$$
(5)
where $`e^{i\widehat{\varphi }}𝑑\xi \chi ^{}(0,\xi )e^{i\varphi (\xi )}\chi (0,\xi )`$, and $`\varphi (\xi )=𝑑tV_I(t,X_2(t),\xi )`$. Here $`V_I`$ denotes the interaction potential between the system and the environment in the interaction picture, i.e. $`V_I\mathrm{exp}(iH_{env}t)V(X_2,\xi )\mathrm{exp}(iH_{env}t)`$ in which $`H_{env}`$ and $`V`$ represent the Hamiltonian of the environment and its interaction with the system, respectively. Eq. (5) implies that the reduction of the interference can also be ascribed to the accumulating phase uncertainty of the system on the interacting path being subject to an uncertain potential. In this sense decoherence is often referred to as “dephasing”.

## III Decoherence from internal dynamics

### III.1 Two slit interferometry

Now let us consider two path interferometry of a macroscopic object with its complicated ID. The CM and the ID now play roles of a system and an environment, respectively. We take into account the case that the object is moving freely as shown in Fig. 1(b), so that there is no entanglement between the CM and the ID. The motion of the CM is easily described by a plane wave, $`\psi _0(t,X)`$. The CM of the object then shows perfect coherence since the CM dynamics is completely isolated. When one increases the path length of one of the two paths by amount of $`l`$, clear interference pattern is expected from the term $`\psi _1^{}(t,X)\psi _2(t,X)=\psi _0^{}(t+\mathrm{\Delta }t,X)\psi _0(t,X)`$ between two CM states, where $`\mathrm{\Delta }t=l/v`$ ($`v`$ is the velocity of the CM). The final wavefunction is given as

$$\mathrm{\Psi }(t)=\frac{1}{\sqrt{2}}\left[\psi _0(t+\mathrm{\Delta }t,X)+\psi _0(t,X)\right]\chi (t,\xi ).$$
(6)
The time delay $`\mathrm{\Delta }t`$ is only applicable to the CM since the ID independently evolves in time. The overlap of the ID always yields that of the same states, i.e. complete coherence. Without entanglement with the ID the decoherence of the CM cannot be generated.

Let us consider now the case an object goes though nonlinear external potentials along the paths as shown in Fig. 1(c). We assume these nonlinear potentials are equivalent for the two paths. In a two slit interferometry, in general, the two slits are made to have the same geometry. Even though the CM motion is entangled with the ID during the passage through such nonlinear region, this does not generate any decoherence. The reason is that the evolution of the ID is always equal for the two paths, so that $`\chi _1(t,\xi )=\chi _2(t,\xi )`$, i.e. $`𝑑\xi \chi _1^{}(t,\xi )\chi _2(t,\xi )=1`$. The entanglement with internal environment makes the phases of each of the partial waves of the CM uncertain in viewpoint of Eq. (5), but does not alter the relative phase. In Fig. 1(d), we introduce additional delay for one of the paths. When the delay is given after the passage through the nonlinear region, the situation exactly corresponds to the case described in Fig. 1(b). The only difference is that one starts not with an initial ID state $`\chi _1(0,\xi )`$ \[$`=\chi _2(0,\xi )`$\] but with $`\chi _1(t_0,\xi )`$ \[$`=\chi _2(t_0,\xi )`$\], where $`t_0`$ is the time when the object departures from the nonlinear region. When the delay is given before, one can see that the situation is also the same as the case shown in Fig. 1(b) by considering the argument related to Eq. (6).

It should be noted that in usual decoherence generated from an outer environment it is not easy to find the case that the two paths have the same environment. If the decoherence occurs mainly by interacting with other particles, the system going through two paths accumulates different random phases from scattering with different particles of different states. The system going through the two paths thus see different environments. This is the reason why the decoherence from an outer environment has been dealt with the setup shown in Fig. 1(a). One example that the same environment is applied to the two interfering waves is the interaction of an interfering electron with zero-point (or vacuum) fluctuation, where the electron does not decohere Kumar87 ; Rammer88 ; Imry02 .

### III.2 Quantification of decoherence from internal dynamics

In the above discussion it has been shown that it is not easy to see the decoherence generating from the ID with usual simple geometry of interferometry. The only way to observe the decoherence from the coupling with ID is that the ID’s should see different nonlinear interactions for each path. Without loss of generality this situation can be represented as the case that only one path contains the nonlinear region. The overlap of the ID can then be written as

$$\chi (0,\xi )\left|e^{iH_{in}\mathrm{\Delta }t}\widehat{T}\mathrm{exp}\left[i_0^{\mathrm{\Delta }t}𝑑t\left(H_{in}+\mathrm{\Gamma }\right)\right]\right|\chi (0,\xi ),$$
(7)
where $`\mathrm{\Delta }t`$ and $`\widehat{T}`$ denote the interaction time within the nonlinear potential and the time ordering operator, respectively. Asymmetric geometry of interferometry is sometimes useful, for example, measurement of phase of the transmission coefficient through a quantum dot Schuster97 , where a quantum dot is plugged into one of the arms of an Aharonov-Bohm ring. From the Aharonov-Bohm oscillation one can determine the phase shift of electron passing through the quantum dot.

The quantity given in Eq. (7) is known as so called fidelity Nielsen00 . The decay of such a quantity determines the decoherence rate. One important remark is that even for the ID with a few degrees of freedom the fidelity decays exponentially if its dynamics is chaotic Peres84 ; Jalabert01 . It opens possibility that the decoherence can occur in molecules consisting of even several atoms from entanglement with its ID Park03  (See also Adachi88 ; Kolovsky94b ; Kubotani95 ; Nakazato96 ; Tanaka02 ; Cohen04  for decoherence generated from chaos) in a certain special condition. When a single coherent state, a minimal wavepacket, is chosen as an initial state of ID, one can expect golden rule or Lyapunov decay depending on the strength of the coupling Jacquod01 , in which the fidelity decay does not much depend on the initial condition for a given energy as far as the completely chaotic dynamics is concerned for the ID. This situation corresponds to the ID governed by rather small number of degrees of freedom at low temperature. As we mentioned in the beginning, however, for small molecules it is hard to expect the coupling between the CM and the ID. To see this effect in the systems with small degrees of freedom we need something different from usual molecules.

For rather bigger systems, which we are interested in, it is not easy to directly calculate how much the interaction between the CM and the ID make an influence on the state of the ID because the internal degrees of freedom consist of many particles. First let us consider the zero temperature case. If the interaction is strong enough to generate any kind of elementary excitations such as phonons, charge density waves, magnons for magnetic systems, and so on, then the CM will lose his coherence completely. Nothing happens for the system going through the path 1 in Fig. 1, while the state of the ID through the path 2 in Fig. 1 is changed from the ground state to the excited state of the elementary excitation. By checking the state of the ID one can see which path the system go through. It is nothing but a complete decoherence. In this case it is crucial to know whether the excitation is gapless or not. At finite temperature the situation is much more complicated. In the beginning we assumed that there is no other sources of decoherences except those originating from the ID. Unfortunately this is no longer true since an object with finite temperature is always coupled with outer environment by emitting black body radiation or possibly other radiation from thermal vibration. It is another issue how much this effect degrades the coherence of the motion of the CM. Now one needs to compare the decoherence from the coupling with the outer environment and that with the ID. If the radiation is not so harmful for the coherence, e.g. the case that the wavelength of the radiation is much larger than the difference of two paths, it is meanigful to calculate Eq. (7). This calculation can be done by using the field theoretical technique once the time dependent interaction $`\mathrm{\Gamma }`$ is known. Surely it is still hard task. We do not want to calculate it in this paper, but only point out the possible existence of the decoherence generating from the coupling with the ID.

Another remark on Eq. (7) is its physical interpretation in terms of quantum measurement. Here, the total system including CM and ID is isolated from the external world. The decoherence is generated from entanglement only within a macroscopic object itself. In this sense no external observer exists. No information is transferred from the object to outside. Since a system (the CM) does not know the difference between an internal and an outer environment, the ID plays a role of an observer. Information thus flows from the CM to the ID. Following Zurek’s argument Zurek03 , one can say the ID is watching the CM.

### III.3 Weak localization in disordered systems

Finally, let us consider one interesting example, weak localization in disordered systems Bergman84 . Two time reversal paths multiply reflecting from the random scatterers leads to the localization of a wavefunction due to their constructive interference as shown in Fig. 2(a). Such localization is fragile for both the decoherence and the perturbation breaking time reversal symmetry. A random potential generates complicated and chaotic dynamics, which can give rise to the entanglement between CM and ID. At first glance the decoherence from the ID is hardly expected to arise since the total system has time reversal symmetry. One might think that the state of the ID of the CM rotating clockwise must be the same as that of counter-clockwise. It can be shown, however, that the coupling to the ID can destroy the weak localization by generating the decoherence from the ID.

To prove the appearance of the decoherence let us consider the overlap between two ID, namely $`\chi _1`$ and $`\chi _2`$, of the time reversal paths of CM after a round trip along a closed loop. Since we are interested only in the ID, the influence from its coupling to the CM can be regarded as a time dependent external perturbation, i.e. $`\mathrm{\Gamma }(X_{1(2)}(t),\xi )\mathrm{\Gamma }_{1(2)}(t,\xi )`$. Note $`\mathrm{\Gamma }_{1(2)}(t+\tau ,\xi )=\mathrm{\Gamma }_{1(2)}(t,\xi )`$, where $`\tau `$ denotes the duration time taken for a round trip around the closed loop. The final state $`\chi _{1(2)}(\tau )`$ is then given as

$$\chi _{1(2)}(\tau ,\xi )=\widehat{T}\mathrm{exp}\left\{i_0^\tau 𝑑t\left[H_{in}(\xi )+\mathrm{\Gamma }(X_{1(2)}(t),\xi )\right]\right\}\chi _0.$$
(8)
Even though $`X_1(t)=X_2(\tau t)`$ holds, i.e. $`\mathrm{\Gamma }_1(t,\xi )=\mathrm{\Gamma }_2(\tau t,\xi )`$, one finds $`\chi _1(\tau )\chi _2(\tau )`$ due to existence of the time ordering operator $`\widehat{T}`$. To make it more clear, let us consider a simple example: three scatterers well localized in space as shown in Fig. 2(b). During free propagation between two scatterers the CM is decoupled from the ID. We assume the process of collision with the scatterers is short enough to be described by delta-function in time. The interaction term $`\mathrm{\Gamma }_1`$ for the clockwise propagation can then be given as

$$\mathrm{\Gamma }_1(t,\xi )=f_1(\xi )\delta (t\tau _1)+f_2(\xi )\delta (t\tau _2),$$
(9)
where $`\tau _1`$ and $`\tau _2`$ are the collision times upon the first and the second scatterer, respectively, and $`0<\tau _1<\tau _2<\tau `$. After one round trip, the states of the ID for the clockwise and the counter-clockwise $`\chi _1`$ and $`\chi _2`$ are respectively given as

$`\chi _1(\tau )=U_3K_2U_2K_1U_1\chi _0`$
$`\chi _2(\tau )=U_1K_1U_2K_2U_3\chi _0,`$ (10)
where by using the eigenstates of $`H_0`$ ($`H_0|i=E_i|i`$) one obtains $`U_{n,ij}=\mathrm{exp}[iE_i\mathrm{\Delta }t_n]\delta _{ij}`$, and $`K_{n,ij}=i|\mathrm{exp}[if_n(\xi )]|j`$. Here, $`\mathrm{\Delta }t_1=\tau _1`$, $`\mathrm{\Delta }t_2=\tau _2\tau _1`$, and $`\mathrm{\Delta }t_3=\tau \tau _2`$. It is obvious that in general $`U_3K_2U_2K_1U_1U_1K_1U_2K_2U_3`$ since $`K_1`$ and $`K_2`$ are not diagonal. Consequently one can find that in general $`\chi _1(\tau )|\chi _2(\tau )<1`$. The weak localization of the CM of a macroscopic object in disordered potentials can be destroyed due to the coupling to the ID.

## IV Summary

In summary, we have investigated the decoherence generated from the internal dynamics of a macroscopic object. In a usual setup of two slit interferometry, it is hard to expect the appearance of such decoherence. Only asymmetric geometry of the interfering paths containing anharmonic external potential allows one to observe the decoherence from the internal dynamics. Such decoherence can then be measured by the fidelity given in Eq. (7). In this case, the internal degrees of freedom of a macroscopic object are watching its center of mass motion. The weak localization of the center of mass motion of a macroscopic object in disordered potentials can also be destroyed by such decoherence without any external perturbation breaking time reversal symmetry.

## Acknowledgments

SWK would like to thank Fritz Haake and Klaus Hornberger for helpful discussions. The part of this work has been done when SWK has participated in the focus program of Asia Pacific Center for Theoretical Physics (APCTP) in Pohang, Korea.


text:
# Improved neighbor list algorithm in molecular simulations using cell decomposition and data sorting method

## 1 Introduction

Some molecular simulation techniques such as molecular dynamics and Monte Carlo method are widely used to study the physical properties and chemical processes which contain a large number of atoms in statistical physics, computational chemistry, and molecular biology . All these methods involve evaluation of the total interatomic potential energy $`V_{\text{tot}}`$ of $`N`$ atoms and its gradients. The potential energy contains various interatomic interactions in the physical system, and is usually the function of internal coordinates of atoms. For example the potential energy of liquids and gases is often described as a sum of two-body (or pairwise) interactions over all atom pairs. A common choice of two-body interatomic interaction expression is Lennard-Jones potential function, which is a simple function of the distance $`r_{ij}`$ between atom $`i`$ and $`j`$, and is shown as follows,

$$V_{\text{LJ}}(r_{ij})=4ϵ\left[\left(\frac{\sigma }{r_{ij}}\right)^{12}\left(\frac{\sigma }{r_{ij}}\right)^6\right].$$
(1)
The total potential energy is a sum of two-body interactions over all atom pairs,

$$V_{\text{tot}}=\frac{1}{2}\underset{i=1}{\overset{N}{}}\underset{\begin{array}{c}j=1\\ ji\end{array}}{\overset{N}{}}V_{\text{LJ}}(r_{ij}).$$
(2)
In molecular dynamics simulation, evaluation of Eq.(2) and its gradient usually costs most of CPU time. Apparently, direct calculation of Eq.(2) requires $`N^2`$ steps. If we use Newton’s third law, the total calculation steps can be decreased to $`N(N1)/2`$. Obviously it is formidable to carry out such a calculation when there are many atoms in the system, and some methods are strongly needed to reduce the redundancy in evaluation of Eq.(2).

Firstly a cutoff distance $`r_{\text{cut}}`$ is introduced in potential functions, and both potential functions and their gradients beyond the cutoff distance are assumed to be zero. This treatment can reduce the computing time greatly by neglecting all atoms beyond the cutoff distance, since interactions between these atoms are zero and needn’t to be considered. However, straightforward determination of which atoms are within cutoff distance needs to evaluate all interatomic distances over all atom pairs, and this procedure scales $`O(N^2)`$ as the system size.

Effective reduction of unnecessary interatomic distance evaluation can be accomplished by Verlet table algorithm  and cell linked list algorithm . However, there is a tradeoff between overhead for maintaining neighbor list table and reduction of unnecessary interatomic distance calculation.

Verlet table algorithm and cell linked list algorithm have been intensively studied and have shown significant reduction in total computing time . Glikman et al proposed an alternative Verlet table method called *relational method* which updates the neighbor list recursively and is called in every timestep , Eisenhauer considered the locality of the neighbor list and used a simple slab-based decomposition in his work . In 1994, Hansen and Evans employed cell method in conventional Verlet table algorithm to accelerate the neighbor list construction speed , and Frenkel also considered combination of two conventional neighbor list algorithms in his book . In our work improved neighbor list algorithm is proposed and the overhead to maintain the neighbor list table has been reduced to order $`O(N)`$. The details and benchmark of this algorithm are given in the paper.

## 2 Conventional neighbor list algorithms and the improved algorithm

For the details of conventional Verlet table and cell linked list algorithms one can refer to Allen and Tildesley’s book , and here we only mention some fundamental and important facts. Some graphs in this section are drawn in 2D for the convenience of illustration, however, all related discussions can be easily generalized to 3D systems.

### 2.1 Conventional Verlet table algorithm

The basic idea of Verlet table algorithm is to construct and maintain a list of neighboring atoms for each atom in the system . During the simulation, the neighbor list will be updated periodically in a fixed interval, or automatically when the displacements of some atoms are larger than a certain value . In this algorithm the potential cutoff sphere with radius $`r_{\text{cut}}`$ is surrounded by a “skin” $`r_s`$, to give a larger sphere with radius $`r_\text{s}+r_s`$ . When constructing a neighbor list for an atom $`i`$, another atom $`j`$ is considered as a “neighbor” if the distance between them $`r_{ij}r_{\text{cut}}+r_s`$. It should be noticed that the “skin” $`r_s`$ should be large enough, so that between two times of neighbor list reconstruction no atom can penetrate through the skin into the cutoff sphere of another atom if it is not in the neighbor list of that atom.

Conventionally in this algorithm, interatomic distances of all atom pairs need to be calculated, so the total CPU time scales $`O(N^2)`$ as the number of atoms. However, evaluation of the total potential energy and its gradient using the neighbor list is efficient because only atoms appearing in the list, i.e., those in the sphere of $`r_{\text{cut}}+r_s`$, are checked, thus the overall procedure scales $`O(NN_{\text{NL}})`$, where $`N_{\text{NL}}`$ is the average number of neighbors of an atom in the material and is independent to the system size $`N`$.

Verlet table algorithm has been proven to be efficient when the number of atoms is relatively small and the atoms move slowly. Its main advantage is the high efficiency of potential/force evaluation using the neighbor list, as the atoms in the list are only a few more than actual needed (which $`r_{\text{cut}}<r_{ij}<r_{\text{cut}}+r_s`$). On the other hand, its main drawback is that construction of the neighbor list scales $`O(N^2)`$ (needs $`N(N1)/2`$ steps to build the list for all atoms). Moreover, with the increase of atoms’ mobility, either the “skin” or the frequency of reconstructing the neighbor list must increase. Both of them make the overall simulation time increases dramatically.

### 2.2 Conventional cell linked list algorithm

The cell linked list algorithm is effective when the number of atoms is large. In this algorithm, the simulation domain is partitioned into several cells, and the edge of cells is equal to or larger than $`r_{\text{cut}}`$. All atoms are assigned to these cells by their positions, at the same time a linked list of the atom indices is created. At the beginning of a simulation, an array containing a list of cell neighbors for each cell is created, and this list remains fixed unless the simulation domain changes during the simulation .

A cell $`m`$ is considered as a neighbor of another cell $`n`$, if there are a point at $`𝐫_i`$ in cell $`m`$ and another point at $`𝐫_j`$ in cell $`n`$ such that $`|𝐫_i𝐫_j|r_{\text{cut}}`$. Since in the conventional method the edge of the cell is equal to or larger than $`r_{\text{cut}}`$, and considering the periodic boundary condition, there are 8 (for 2D systems) or 26 (for 3D systems) neighbors for each cell. Thus the neighbors of each atom can be found in the cell where the atom located and the neighboring cells.

The construction of “neighbor list”, i.e., assigning each atom to appropriate cells, scales $`O(N)`$, but we can see that a big number of atoms need to be checked in the potential/force evaluation procedure, and this is rather inefficient compared to Verlet table algorithm. A common choice of the cell edge is the potential cutoff distance $`r_{\text{cut}}`$, thus for each atom, all atoms in 27 cells, or in the volume of $`27r_{\text{cut}}^3`$, will be checked in the potential/force evaluation procedure111If Newton’s 3rd law is used and the interatomic interaction is described by a two-body potential function, only half of neighboring cells needs to be searched, in other words, only 14 cells are needed for 3D systems. However, if multiple-body potential functions are adopted, especially those using bond order terms, we cannot safely do this. For example in Brenner potential, calculation of $`N_{ij}^{conj}=1+[_{ki,j}^{carbon}f_{ij}^c(r_{ik})F(X_{ik})]^2+[_{li,j}^{carbon}f_{jl}^c(r_{jl})F(X_{jl})]^2`$, i.e., determination of whether a bond is part of a conjugated system or not, needs to reference all neighbors of every atom. In this case, searching only half of the neighboring cells is not correct. Therefore, in general cases we say all 27 cells need to be searched.. Ideally, only atoms in the volume of $`\frac{4}{3}\pi r_{\text{cut}}^34.189r_{\text{cut}}^3`$ fall in the cutoff distance, and accordingly need to be checked.

However, if a small cell edge is used, volume containing the atoms need to be checked will be dramatically reduced. For example if the cell edge is $`\frac{1}{2}r_{\text{cut}}`$, the volume will be $`(2.5)^3r_{\text{cut}}^3`$, only 57.87% of the previous volume, as Fig. 2 shows for 2D systems.

Furthermore, the cell edge can be very small so that only one atom can be contained in the cell, as described in Ref. . This method is reported to has higher performance than conventional algorithm. However, according to our tests, it is still not as fast as Verlet table algorithm when the number of atoms is relatively small.

From above discussions, we can see that the advantage of the cell linked list algorithm is the fast and efficient construction of the “neighbor list”, but its disadvantage is that too many atoms need to be checked in the potential/force evaluation procedure, and improvements to this algorithm are not competitive to Verlet table algorithm with the increasing of implementation complexity when the number of atoms is relatively small.

### 2.3 The improved neighbor list algorithm

Verlet table algorithm will scale linearly as the number of atoms if cell decomposition approach is adopted. Together with data sorting and partial updating method, we get an improved neighbor list algorithm with high performance.

#### 2.3.1 Cell decomposition approach

As cell decomposition approach has been introduced and used previously , here we only give its main features. In this approach two conventional algorithms are combined together. Firstly the simulation domain is partitioned into several cells, and then each atom is assigned to these cells by their positions, in the following the neighbor list is constructed by searching only in neighboring cells, instead of checking all atom pairs in the system. It can be seen that the construction of the neighbor list scales $`O(cN)`$, where $`c`$ is a constant and independent to $`N`$. For a system containing more than 1000 atoms, this approach improves the overall performance significantly.

For systems with high atomic mobility, we find that the cell edge of $`\frac{1}{2}r_{\text{cut}}`$ gives best performance after carrying out many test runs on some different kinds of computers.

#### 2.3.2 Partial updating of the neighbor list

In certain systems where a small number of atoms have high mobility and most others only oscillate at equilibrium positions, the choice of skin in Verlet table algorithm can be very tricky. Small value ensures that atoms which are stored in the list and beyond the cutoff distance ($`r_{\text{cut}}<rr_{\text{cut}}+r_s`$) are as less as possible, however, updating of neighbor list will become frequent. On the other hand, the potential/force evaluation efficiency is lowered when using large skin value, it can be tought to find a good value.

Fortunately, this problem can be solved in the framework of cell decomposition approach. As searching for neighbors is confined to several cells, after the neighbor list is out-of-date and needs to be updated, reconstruction procedure can be carried out for atoms in relevant cells only, instead of for every atom in the system.

In this method, each cell is accompanied with a “dirty flag”, which denotes if the neighbor list corresponding to this cell needs to be updated. At every timestep all atoms’ accumulative displacements will be checked in each cell, if the sum of the largest two displacements is greater than $`r_s`$, then the “dirty flag” of this cell and neighboring cells will be turned on, and the neighbor list construction procedure is carried out for these cells. After the construction procedure these flags are turned off.

For aforementioned system this method effectively improves the overall performance because the CPU time on updating the neighbor list is minimized. In addition, the skin can be small enough so that the potential/force evaluation efficiency is maximized.

It needs to be noticed that in this method the neighbor lists of different atoms should be stored separately (e.g., the neighbor list array can be defined in the form of nlist\[N\]\[MAXNL\], where N is the number of atoms and MAXNL is the maximum number of neighbors of an atom). This, however, may slightly increase the memory requirement although this is not serious in modern computers where the internal memory is at least hundreds of MB.

#### 2.3.3 Acceleration of data access by data sorting

Considering pipeline architectures of modern CPUs, further effort can be taken to maximize the computing performance: sorting the storage sequences of atoms in the memory, thus the memory locations of atoms which in the same cell or neighboring cells are also close to each other, then the data can be loaded and cached in the CPU more efficiently.

For a better understanding of this method, we can consider the simulation of gas and liquid materials which have high atomic mobility. In the beginning, after the molecular structure is initialized, the data of positions, velocities and accelerations are well ordered and stored in the memory continuously. Sometimes data of all neighboring atoms can be loaded into CPU data cache if original data is well organized. But when the simulation is going on and the atoms are moved here and there, memory locations of these data become more and more disordered, and the neighboring atoms can be seldom loaded in the same cache line. Then the CPU is hard to find the data of neighboring atoms in the data cache, and has to stall the execution and fetch required data into the data cache. However, irrelevant data along with fresh loaded data pollute the data cache, therefore the data cache miss rate can be high. This phenomena can be detected in a long time simulation, as an example shown in Fig. 4, where a very significant performance degradation can be seen.

The solution of this problem is rather straightforward, i. e., sorting the data of atoms by their positions and making the memory locations of the neighboring atoms as close as possible. By adopting this method, no explicit performance degradation can be seen in long time simulations. The procedure to carry out the data sorting is shown in Algorithm 1. In this procedure, firstly we find the longest direction of the domain, and partition the domain into several layers along this direction, lastly the data of all atoms are stored layer by layer. Please note that a “weak” sorting is used, it means that data of atoms inside a layer may be disordered. It can be seen that “weak” data sorting procedure scales $`O(N)`$ as the number of atoms.

Together with the cell decomposition approach in section 2.3.1, the data locality is enhanced. Thus when parallelizing the simulation program on SMP platforms, the data can be easily and well partitioned, and the neighbor list construction can be carried out by each CPU in the computer simultaneously, thus a high parallelization execution can be achieved.

In this work, overall procedure for constructing the neighbor list222Full neighbor list construction is carried out at the first timestep and at every timestep after the data sorting, and partial updating is performed between two times of data sorting. is shown in Algorithm 2. Here it should be mentioned that data sorting procedure may not be carried out for low atomic mobility systems.

## 3 Results and discussions

A molecular dynamics simulation program using Lennard-Jones two-body potential function is developed to compare the performance of three different neighbor list algorithms. The benchmarks are carried out on a Compaq Alpha Server DS20 with two EV67/667 MHz processors and Tru64 5.1A operating system installed, and the program is written in Fortran 90 and compiled by Compaq Fortran compiler V5.5–1877. We also run the same benchmarks on a PC with one Intel Pentium III 866 MHz CPU and Red Hat Linux 8.0 installed, and a HP RX2600 workstation with two Intel Itanium2 900 MHz processors and Red Hat Linux Advanced Workstation installed. Three CPUs have different architectures: Alpha EV67 is a typical RISC (Reduced Instruction Set Computing) CPU, Intel Pentium III is a CISC (Complex Instruction Set Computing) CPU, while Itanium2 is declared as a brand new architecture named EPIC (Explicitly Parallelized Instruction Computing). The performance on three platforms are different, however the difference among three algorithms are qualitatively similar. All results presented in this section are based on the benchmark on Alpha Server DS20.

In order to measure the performance of an algorithm quantitatively, a new unit named *atom$``$step/second* is defined. It can be simply calculated by multiplying the number of atoms and number of timesteps divided by the number of CPU time elapsed in second. Larger value of this quantity stands for better performance. If computer hardware architectures are identical, and the simulation program scales $`O(N)`$ perfectly as the system size, then this quantity is proportional to the CPU performance.

In our benchmark, firstly some Argon atoms are randomly placed in the domain according to the predetermined density. Then the molecular dynamics simulation in canonical ensemble is performed, and the number of steps is $`10^2`$ for $`10^4`$ atoms and more, or $`10^3`$ for $`10^310^4`$ atoms, or $`10^4`$ for 999 atoms or less. The temperature of the system is 300 K.

The time integration algorithm is implementated by velocity Verlet scheme, and the canonical ensemble simulation is implemented by using Nośe–Hoover thermostat . Real physical units, instead of reduced units, are used. We collect the parameters in the simulation and show in Table 1.

The simulation program is parallized using OpenMP technique . Serial and parallel version of the program are compiled by turning on or off the relevant compiling options of OpenMP333In Compaq Fortran compiler on Tru64/Alpha option -omp enables OpenMP, in HP Fortran compiler on HP-UX/Itanium option +Oopenmp/+Onoopenmp enables/disables OpenMP, and in Intel Fortran compiler option -openmp enables OpenMP..

In order to verify the improved neighbor list algorithm, all neighbor lists are dumped to disk files and compared with those in the execution with Verlet table algorithm. In the verification stage, some statistical quantities, such as total potential energy, total kinetic energy, transient temperature of system and the trajectory of all atoms have been recorded in the interval of 10 timesteps, and data generated from three algorithms are compared and ensure they differ in round-off errors only. Intensive tests show that the neighbor lists from the improved algorithm are exactly the same as those from the other two algorithms, and three different algorithms give same results precisely.

Several simulations with different number of atoms and the same density are carried out, and the performance is calculated and recorded.

Firstly we measure the neighbor list constructing time with different number of atoms for three algorithms, and results are shown in Fig. 5. In order to measure the CPU time accurately, a very large number of atoms are used. From the results we can see that the neighbor list constructing time of Verlet table algorithm scales $`O(N^2)`$ as the number of atoms, while the other two algorithms scale $`O(N)`$. We can also see that the neighbor list constructing time of the improved algorithm is slightly larger than that of the cell linked list algorithm. However, the neighbor list needs to be constructed in every timestep in the cell linked list algorithm, but it needs not in our improved algorithm, thus the CPU time used to maintain the neighbor list in the improved algorithm is much less than those in the other two algorithms.

Then we carry out the real benchmark and calculate the overall performance of the simulation with different algorithms. Results are shown in Fig. 6 (single processor results) and Fig. 7 (dual-processor results)

From the results of Fig. 6, we can know that the proposed improved algorithm improves the overall performance of the molecular dynamics simulation significantly. When the system is small, the performance of the improved algorithm is as high as that of Verlet table algorithm. As the system size increases and there are more and more atoms, the performance of Verlet table algorithm decreases rapidly, but the performance of the improved algorithm becomes even better. As the system is as big as $`7\times 10^5`$ atoms, the performance of Verlet table algorithm becomes very bad, and is much lower than that of the improved algorithm. On the other hand, although the cell linked list algorithm can handle a very large system due to its small memory requirement, its performance is much lower than that of the improved algorithm. For large systems, the improved algorithm in our work is about 2 $``$ 3 times faster than conventional neighbor list algorithms.

From Fig. 7 we can see that the improved algorithm in SMP platforms exhibits higher performance than the other two algorithms. The conclusions in above discussions for single-processor are also valid for dual-processor.

## 4 Conclusions

Nowadays as the continuously increasing of computing power, huge and complicated molecular simulations will be attempted and involve a larger number of atoms and more complex potential functions. The expectation of running molecular simulations faster and easier for larger systems on existing platforms makes it important to improve the neighbor list algorithm in order to reduce the unnecessary interatomic distance calculations. In this paper, we proposed an improved order $`O(N)`$ neighbor list algorithm which incorporates the advantages of conventional Verlet table and cell linked list algorithms. In the new algorithm, cell decomposition approach is adopted to accelerate the neighbor list construction speed, and data sorting method is used to maximize the CPU data cache hit rate. In addition, partial updating of the neighbor list is introduced to minimize unnecessary neighbor list reconstructions. We carried out the benchmarks to these algorithms using molecular dynamics simulations with Lennard-Jones potential function, and compared the performance of the improved algorithm with that of the other two algorithms. The results show that the improved algorithm outperforms conventional Verlet table and the cell linked list algorithms by a factor of $`23`$ in single-processor and SMP platforms.

## Acknowledgments

This work is supported by the Singapore-MIT Alliance.