text:
# Size Effects in Heavy Ions Fragmentation

## 1  Introduction

M. Rousseau et. al. have studied light charged particle emission from a compound nucleus of $`{}_{}{}^{40}Ca`$ obtained from $`{}_{}{}^{28}Si`$+$`{}_{}{}^{12}C`$ reaction, providing evidence of binary emission of $`{}_{}{}^{8}Be`$, in a data fitting scenario formed by a sequential evaporation process, a subsequent fragment emission from the $`{}_{}{}^{40}Ca`$ compound nucleus and finally the residual nucleus freeze-out .

In the case of inverse kinetic reactions, fragmentation occurs in two stages, namely a first pre-equilibrium stage where participants interact and deposit excitation energy in the spectators, and a second stage where the excited spectator residual decays .

Considering small excitations, there is a delay in fission fragment emission in the range of 1000 fm/c and for larger excitations the charged particles are emitted at the same time as heavy fragments ($`Z_{imf}`$= 8- 16), with an integral rupture time in the range of 200 fm/c after thermalization . Kolomietz et. al. have obtained “rise-plateau-rise isobaric caloric curves in the Thomas-Fermi approximation with an effective force SkM, showing that in the plateau there is liquid-gas coexistence .

Gourio et. al. have obtained light particle emission times around 200 fm/c, using two-particle correlation functions in order to analyze INDRA and GANIL data . D’Enterría et. al. have obtained similar caloric curves studying Brehmstahlung photon emission in Heavy Ion Collisions using a thermal model . D’Agostino et. al. have shown the usefulness of an abnormal change in kinetic energy variance as a first order phase transition signature . Nuclear collisions for energies in tenths of $`meVs/nucleon`$ produce excited systems that break-up into several intermediate size fragments . These multifragmentation phenomena, similar to those occurring in inertial confinement reactions  and in the nanostructures surface synthesis , could be the adaptation of a macroscopic phase transition to a finite and transient system.

In the realm of nuclear physics, the possibility of reaching a critical behavior in heavy ions collisions has motivated several studies on critical exponents. The first one focused on proton-Xe and proton-Ar collisions by Purdue group  and others afterwards . More recently, modern detection technology allowed caloric curve experimental computation  , i.e. the relation between system temperature and excitation energy at fragmentation time. Nevertheless, the influence of residual finite size, a transient evolution and experimental limitations due to observing only final fragments, turn out to be problematic.

As explained by Pochodzalla y Trautmann  , there are difficulties in the reconstruction of laid down energy via exit channels and therefore in caloric curve computation. Other problems arise from the fluctuations in system size , “side feeding” effects on final mass distribution, and in the use of final spectra modified by collision evolution .

As a matter of fact, experiments using different “thermometers” have reached contradictory caloric curves such as the typical rise-plateau-rise plot , a peculiar “rise-plateau” caloric curve , and even a “rise-rise” caloric curve without the plateau characteristic of a first-order phase transition .

The suppression of the final temperature growth could be related either to a low density rupture , or to a high fragmentation energy destroying hot intermediate fragments .

With so many variables affecting the nuclear caloric curve determination, the question to answer is whether it is possible to extract this caloric curve in nuclear collisions. More recently, this was dealt with a Molecular Dynamics study in nuclear systems inside of a container  and for a freely expanding classical system .

Even when those computations can not be compared, neither between them nor with experiments, those exercises were able to obtain a caloric curve. Sugawa and Horiuchi , used Antisymmetrized Molecular Dynamics, in order to study a uniformly excited, contained system with constant pressure, and obtained a “rise-plateau-rise” caloric curve. On the other side, Strachan and Dorso  used a uniformly excited Lennard-Jones system with free expansion and obtained a caloric curve with a “rise-plateau” shape. The difference between these results arises from the collective expansion manifest in finite systems (but absent in infinite and contained systems), acting as an energy reservoir and avoiding a final temperature rise .

And Gobet et. al. have obtained event by event “rise-plateau-rise” caloric curves for collisions of hydrogen clusters ions with Helium targets, considering this caloric curve as a first-order phase transition signature in a finite system. The plateau is interpreted as latent heat associated with this phase transition .

Notwithstanding their differences, these exercises suggest that fragmenting small systems may provide information about the caloric curve. Nevertheless, in regards of the fundamental role played by the geometry of the fragmenting system, i.e. contained systems versus free-expansions systems, a subtle question might be that of the role played by the collision in caloric curve computations. Namely, the influence of induced large correlations due to non linear dynamics on the obstacles found to compute the caloric curve.

In this study a computation is performed that reproduces experimental collisions in order to obtain the caloric curve and examine the dependence of the limit temperature on the residual size.

We use the weaponry developed in  to focus on colliding excited systems in order to obtain the caloric curve and its dependence on the residual size. The manuscript is organized in the following way: after describing the model used, the fragment recognition algorithm and the persistence in section II, Section III establishes an effective way to compute the collision dynamic stages and section IV studies the temperature in the fragmenting system participant region. Caloric curve dependence on residual size is shown in section V, and some final conclusions close the article in Section VI.

## 2  Latino Model.

LATINO Model  uses the semi-classical approximation to simulate Heavy Ions collisions via a binary interaction Pandharipande potential. This potential is formed by a linear combination of Yukawa potentials, with the coefficients fixed in order to reproduce the properties of nuclear matter. Clusters are detected with an Early Cluster Recognition Algorithm that optimizes the configurations in energy space. Ground states are produced generating a random configuration in phase space, gradually reducing the velocities of the particles confined in a parabolic potential, until the theoretical binding energy is reached.

Due to obscurities in the break up process, a complex model is needed with the auxiliary tools required. Collision evolution is modeled with a Molecular Dynamics code, though since Molecular Dynamics operates at a nucleonic level, it is necessary to transform particle information in terms of fragment information via a fragment recognition algorithm. Using caloric curve as a phase transition signature, requires the knowledge of the break up time, using a property known as partition “persistence”.

Molecular Dynamics advantages to study nuclear collisions have been established previously. In this work we study the time evolution of central collisions numerical simulation. Target consists of a 3D droplet in its “ground state”. The projectile is a randomly oriented droplet in its “ground state”, boosted on the target with different kinetic energies.

Numerical integration of the equations of motion is performed through a Verlet algorithm at time intervals that ensure energy conservation at least in a 0.05 %. Fig. 1 shows the evolution of a typical experiment.

The range for the projectile kinetic energy starts from $`E_{beam}=800`$ $`MeV`$ in the center of mass reference system with about two hundred collisions for each projectile energy. Covering a wide energy range, we obtain different dynamical evolutions: starting from events where the projectile is absorbed by the droplet surface up to events with an exponential decay mass spectra. As shown in Fig. 2, the shapes of the spectra converge to a power law in this range of intermediate energies. The temperature of the participant region is computed using Kinetic Gas Theory and system excitation is calculated with the energy deposited in the residual. When the projectile energy is increased, the multiplicity is also increased as the system enters into the phase coexistence region (Fig. 3b) .

With the fragments detected in phase-space, Mean Velocity Transfer is computed in the following way :

$$MVT_i=\underset{k}{}\left|v_{k,i}(t+dt)v_{k,i}(t)\right|$$
(1)
### 2.1  Fragment Recognition

In order to obtain information about the fragments, Molecular Dynamics microscopic data must be analyzed through a Fragment Recognition Algorithm. Although there are many, we will now describe the one used here.

In the “Minimum Spanning Tree in Energy” algorithm (MSTE), a given set of particles i,j,…,k belongs to the same cluster $`C_i`$ if: $`iC_i,`$ exists $`jC_i/e_{ij}0`$ where $`e_{ij}=V(r_{ij})+(𝐩_i𝐩_j)^2/4\mu `$, and $`\mu `$ is the reduced mass of the couple i,j. MSTE searches configuration correlations between particles, considering their relative moment. Due to its sensitivity to recognize early particles, it is extremely useful to study the pre-equilibrium energy distribution of the participant particles.

### 2.2  Partition Persistence

In order to use the caloric curve as a phase transition signature, both system excitation and temperature at phase transition are needed. Hence, break-up time is required, and is calculated in terms of the fragmentation time, $`\tau _{ff}`$, namely a time when fragments evaporate only some monomers. $`\tau _{ff}`$ can be detected comparing partitions at different times, which can be done using the “Microscopic Persistence Coefficient”, $`P`$ :

$$P[X,Y]=\frac{1}{_{cl}n_i}\underset{cl}{}n_i\frac{a_i}{b_i}$$
(2)
where $`X\{C_i\}`$ and $`Y\{C_i^{}\}`$ are different partitions, $`b_i`$ is the number of particle pairs in the cluster $`C_i`$ of partition $`X`$, $`a_i`$ is the number of particle pairs belonging to cluster $`C_i`$ that remain paired in a given cluster $`C_j^{}`$ of partition $`Y`$, $`n_i`$ is the number of particles in cluster $`C_i`$.

$`P[X,Y]`$ is equal to 1 if the microscopic composition of partition $`X`$ is equal to that of partition $`Y`$, and tends to 0 when none of the particles belonging to a given cluster $`X`$ remain together in any cluster of $`Y`$.

It is convenient to study the time evolution for the quantities:

$$\widehat{P}^+\left[X(t)\right]P[X(t),X(t\mathrm{})]_{events}$$
(3)
$$\widehat{P}^{}\left[X(t)\right]P[X(t\mathrm{}),X(t)]_{events}$$
(4)
$$\widehat{P}^{dt}\left[X(t)\right]P[X(t),X(t+dt)]_{events}$$
(5)
where $`X(t)`$ represents a partition computed at time t, $`X(t\mathrm{})`$ is an asymptotic partition, and $`_{events}`$ represents an average on the total set of collisions. In simple terms, $`\widehat{P}^+\left[X(t)\right]`$ determines whether all particles bound at time t remain bound asymptotically.

In a similar fashion $`\widehat{P}^{}\left[X(t)\right]`$ measures the reciprocal value , $`i.e.`$, the degree in which an asymptotical partition is contained in a partition occurring at time t. In conjunction, $`\widehat{P}^+`$ y $`\widehat{P}^{}`$ can be used to analyze the evolution of the microscopic composition of the most bound partition, until it reaches its asymptotic form.

Finally $`\widehat{P}^{dt}\left[X(t)\right]`$ provides information on the given partition $`activity`$, and can be thus used to define the fragmentation time, $`t_{ff}`$, once it reaches some degree of stabilization.

It is useful to work with some $`normalized`$ quantities:

$$P^+\left[X(t)\right]=\widehat{P}^+\left[X(t)\right]/P[X(t\mathrm{}),X^{^{}}(t\mathrm{})]_{events}$$
(6)
$$P^{}\left[X(t)\right]=\widehat{P}^{}\left[X(t)\right]/P[X(t),X^{^{}}(t)]_{events}$$
(7)
$$P^{dt}\left[X(t)\right]=\widehat{P}^{dt}\left[X(t)\right]/P[X(t+dt),X^{^{}}(t+dt)]_{events}$$
(8)
where $`X^{}`$ is the partition formed by the clusters belonging to partition $`X`$, and where each fragment has evaporated a particle. Normalized quantities compare real values of $`\widehat{P}^{}s`$ with an evaporative level of reference. $`P^+`$ , $`P^{}`$ and $`P^{dt}`$, refer to distinct persistence coefficients, namely forwards, backwards and differential, respectively.

## 3 Dynamic Evolution

Armed with the tools introduced before (MD, MSTE and $`P^+`$ , $`P^{}`$ and $`P_{dt}`$) we proceed to study the dynamical evolution of the collisions described in section II. Analyzing the Mean Velocity Transfer, we can characterize the collision stages and identify early emitted particles. This allows to study the excitation energy and the detection of the fragmentation time.

### A Collision Stages

Two colliding stages are observed, with an initial highly colliding stage produced when the projectile hits the droplet surface and the energy is distributed chaotically. Collisions at this initial stage form a shock wave responsible for the prompt emission of light energetic particles from the surface. As the shock wave travels into the droplet, it produces density fluctuations and internal fractures as a consequence of the momentum transferred and initiates disordered collisions leading to an excitation thermalization.

In order to reach a deeper understanding of this process, we compute the “Mean Velocity Transfer”, defined by:

$$M_j(t)=\underset{i=1}{\overset{N}{}}\left|(\stackrel{}{v}_i(t+dt)\stackrel{}{v}_i(t))\widehat{e}_j\right|_{events}$$
(9)
where j denotes incident and normal directions and $`\widehat{e}_j`$ is a unitary vector in these directions. Fig. 5 shows time evolution of $`M_x`$ for distinct values of $`E_{beam}`$ for Ni+Ag central collision.

Isotropic collisions ($`disordered`$ $`colliding`$ $`mode`$ ) are responsible of the momentum redistribution among the particles remaining in the system. This energy heats the system and constitutes an expanding collective motion. The disordered colliding mode is the only one present for uniformly excited systems where it is responsible of fragmentation and its outwards flow disperses the fragments .

MSTE algorithm can be used to study the size of the MSTE biggest fragment, total multiplicity and persistence coefficients. Figs. 3a, b and c show time evolution of these three quantities. The description suggested is that, due to an initial violent collision, some particles acquire enough energy to be released of the vicinity or interior of a cluster. This reduces the mass of the biggest MSTE cluster in this early stage and increases the total multiplicity. This tendency is maintained until the mean momentum for each particle permits the binding of particles configurationally close. After the initial reduction of the initial biggest fragment size, a coalescent behavior promotes the multiplicity decrease and the biggest cluster size grows until it reaches a maximum. This time signs the end of the initial energy distribution and can be used to define the pre-equilibrium time.

Persistence coefficients can be used to understand how the partition reaches its asymptotic microscopic composition. Fig. 3c shows the time evolution of $`P^+\left[X(t)\right]`$ and $`P^{}\left[X(t)\right]`$ computed using the MSTE partitions for a $`E_{beam}`$ .

Backwards Persistence Coefficient shows initially a decrease due to the fragment formation stage, followed by a subsequent increase due to the dynamic re-absorption lasting until $`t_{pre}`$. High values of $`P^{}`$ indicate that more particles belonging to a given cluster remain paired at time t. Consistently, $`P^+`$ coefficient shows a plateau during this re-absorption stage that is extended until time $`t_{pre}`$, followed by a monotonic increase due to an evaporative dynamics of the MSTE clusters. During this stage, MSTE algorithm detects a large biggest cluster, revealing that the system is still dense.

In summary, an initial stage characterized by the existence of a shock wave releasing particles from the surface, the reduction of the biggest fragment size and an increase of the multiplicity ($`i.e.`$ decreases $`P^{}`$). Afterwards, the shock wave crosses the droplet distributing uniformly the energy ($`M_xM_y`$) and a coalescent behavior is installed signed by the increase of $`P^{}`$ and deleting the memory of the entrance channel, reducing the multiplicity and augmenting the size of the biggest fragment until it reaches a $`80\%90\%`$ of the total mass. This is followed by an evaporative dynamics of the clusters as indicated by a monotonic increase of $`P^+`$. Fig. 4 shows colliding stages for Ni+Ag with an energy equal to 1600 MeV. A first stage ends with a peak in the kinetic energy transported by light particles, signed by a peak in the black curve. A second stage ends with the attenuation of the intermediate fragments production (pink curve). Between these peaks, a peak in Mean Velocity Transfer (orange curve) is observed. Hence, Mean Velocity Transfer (MVT) produces instabilities leading to light particles emission and subsequent intermediate fragment emission. Once the MVT is stabilized, a final stage follows characterized by the emission of light fragments and system freeze-out.

## 4 Temperature

The next step to obtain a caloric curve is the computation of the system temperature at fragmentation time, using the kinetic energy of nucleons in the participant region, relative to the center of mass. This seems to be justified as long as the excited droplets reach thermal equilibrium at break-up . Nevertheless, there exist objections with respect to the use of thermodynamic concepts in small finite transient systems .

Now we define this temperature and study its time evolution.

### A Participant Region Temperature.

In central nuclear collisions, shock waves orthogonal to the beam direction arise, separating compressed nuclear matter from cold nuclear matter. As the projectile penetrates into the target, a region with both high density and excitation is formed. The complex collectively expands in a direction normal to the beam due to the fact that pressure is null in the exterior of the system and therefore expansion in beam direction is decreased.

A discontinuity appears in the contact point between nuclei, with both nuclei traveling in opposite directions with respect to the center of mass, that evolves into two shock waves traveling in opposite directions across the projectile and the target .

A study based on the dynamic model of several fluids, showed that most of the energy associated with the stopping derived from transverse collective motion, is used for mid-rapidity fragment production . In the contact region between both nuclei, frictional forces promote energy dissipation and particle deflection, since projectile participant nucleons interact several times with target participant nucleons. . All along this discontinuity, transport effects are expected leading to a soft change in the properties instead of a discontinuous transition .

In this participant region a temperature can be obtained, considering the nucleons in this region and computing the temperature provided by the kinetic energy of these nucleons measured with respect to the center of mass. Promptly emitted, not reabsorbed nucleons are excluded, using the MSTE Early Cluster Recognition Algorithm for this sake.

### B PEPS and Energy Excitation.

As it happens in the experimental case, the observed presence of promptly emitted particles(PEP) makes possible the $`a`$ $`priori`$ computation of the fraction of beam energy transformed into excitation energy. This complication is absent in computations for both contained or infinite systems, and requires a sensible definition of the excitation energy that should depend on the PEPs.

PEPs can be defined as unbound light clusters (with mass $``$ 4) detected by the algorithm at time $`t_{pre}`$ and that remain unbound at any posterior time ($`i.e.`$ no re-absorption). Once PEPs are defined, their kinetic energy can be used to estimate the energy remaining in the system.

Figs. 6a and 6b show that the transported energy by PEPs and the number of PEPs as a function of the beam energy, respectively. With the number of PEPs and their energy quantified, the energy remaining in the system after time $`t_{pre}`$ can be also obtained.

Fig. 6c shows that energy is deposited in the target as a function of the beam energy. These figures show that the fraction of available energy leaving the system, as a consequence of this prompt emission, is considerable. Even more, the energy remaining in the system shows a saturation behavior indicating that there exists a limit for the quantity of energy that can be transferred in a collision, which does not occur in uniformly excited systems.

### C Fragmentation Time

After the $`pre`$-$`equilibrium`$ time, $`t_{pre}`$ (cf. section IIIA), the energy remaining in the system is completely distributed, and the density correlations $`𝐫𝐩`$ induced by the initial shock wave begin to constitute a collective expansion. Eventually, this expansion will transform the initial fluctuations into fragments well defined in space $`𝐫`$, recognizable for MSTE algorithm. Caloric curves should reflect the system state along the phase transition, and we consider this time as the fragmentation time, $`t_{ff}`$, associated with the stabilization of MSTE density fluctuations

## 5 Caloric Curves and Residual Size

Once the collision dynamic simulation is performed, prompt fragments are detected, fragmentation times are identified and excitation energy as well as temperature are computed, a caloric curve can be obtained. Formally, this quantity investigated experimentally  as well as computationally , is the functional relation between the temperature of the system and its excitation energy. In this study we extend this analysis to Pandharipande droplets excited by collisions.

Figs. 7-9 show the caloric curves computed for a wide energy range of excited collision systems built from participant region temperatures at time $`t_{ff}`$. As can be seen, caloric curves are similar to those obtained from uniformly excited systems . Once again, the relevant characteristic is the almost constant temperature behavior in fragmentation region. In other words, collision data provide a “rise-plateau” caloric curve. Besides, these caloric curves portray a limit temperature that diminishes as the residual size increases.

## 6 Conclusions.

Event by event “Rise plateau” Caloric curves have been obtained for different residual size bins for Ag+Ag Central Heavy Ion Collision. Temperature is computed using Kinetic Gas Theory in the Participant Region, in the moment when the persistence reaches a value close to one. When the residual size increases, a clear decrease of the limit temperature is observed. This is in agreement with recent studies of experimental data performed by Natowitz et.al. .

Therefore we have confirmed using a dynamical model, recent analysis of experimental data, where a decrease of the limit temperature is observed as the residual size is increased.

These results motivate us to perform future studies on the influence of the dynamics in the formation of a plateau characteristic of the coexistence region, for different residual sizes. Authors acknowledge financial support from NSF through PHYS-96-00038 fund and free access to the computational facilities of The U. of Texas at El Paso and UAM-A.


format: text

text:
# Additive entropy underlying the general composable entropy prescribed by thermodynamic meta-equilibrium

## Abstract

We consider the meta-equilibrium state of a composite system made up of independent subsystems satisfying the additive form of external constraints, as recently discussed by Abe \[Phys. Rev. E 63, 061105 (2001)\]. We derive the additive entropy $`S`$ underlying a composable entropy $`\stackrel{~}{S}`$ by identifying the common intensive variable. The simplest form of composable entropy satisfies Tsallis-type nonadditivity and the most general composable form is interpreted as a monotonically increasing funtion $`H`$ of this simplest form. This is consistent with the observation that the meta-equilibrium can be equivalently described by the maximum of either $`H[\stackrel{~}{S}]`$ or $`\stackrel{~}{S}`$ and the intensive variable is same in both cases.

PACS Number(s): 05.20.-y, 05.70.-a, 05.90.+m

Generalized entropic functionals which may yield consistent formulations of thermodynamics  and statistical mechanics  has attracted attention in recent years. Generalization of the standard entropy can also be a useful tool to describe effects of finiteness in physical systems . The standard thermodynamic formalism is based on Boltzmann-Gibbs-Shannon entropy, which is additive. A particular generalization of the standard formalism that has been much studied recently is popularly called Tsallis statistics . This formalism seems to provide an effective description of the meta-equilibrium states of certain complex systems. It is based on a generalized entropy that is nonadditive even when the subsystems forming the composite system are statistically independent. Tsallis-type nonadditivity may be considered as the simplest case of the most general nonadditivity rule for composable entropy, consistent with the existence of thermodynamic equilibrium . This so called equilibrium within the Tsallis approach is actually a meta-equilibrium, which is described by a generalized zeroth law.

On the other hand, it has been suggested that the thermodynamic framework based on Tsallis entropy can be mapped to the framework based on additive (extensive) entropy -. Essentially, it means that a combination of an additive entropy and an appropriate intensive variable preserves the standard thermodynamic relations as well as statistical fluctuations.

In this letter, we derive the additive entropy underlying a composable entropy $`\stackrel{~}{S}`$ by identifying the intensive variable common to subsystems in meta-equilibrium. The case of the most general form for composable entropy as treated in  is interpreted as a monotonically increasing function $`H`$ of the simplest composable entropy. This is consistent with the fact that the intensive variable should be same in the meta-equilibrium state corresponding to the maximum of either $`H[\stackrel{~}{S}]`$ or $`\stackrel{~}{S}`$.

Tsallis entropy is defined as

$$S_q^{(T)}=\frac{_{i=1}^Wp_i^q1}{1q},$$
(1)
where $`p_i`$ is the probability distribution set characterising the discrete microstates of the system and labelled $`i=1,\mathrm{},W`$. The transformation relating additive and nonadditive entropies is given by

$$S_q^{(R)}=\frac{1}{(1q)}\mathrm{ln}[1+(1q)S_q^{(T)}].$$
(2)
The additive entropy so obtained is well known as Renyi entropy .

For better understanding of this transformation, we first derive it directly by applying the generalized zeroth law for meta-equilibrium within Tsallis’ framework. In line with earlier works , we impose additive mean values as constraints. Thus we restrict to a special class of nonextensive models in which the entropy is non-additive, but the energy and other external constraints are additive. Furthermore, it is understood that the zeroth law in the standard Boltzmann-Gibbs formalism identifies an intensive variable common to the systems in mutual thermodynamic equilibrium. The generalized zeroth law serves a similar purpose within Tsallis approach.

Consider the meta-equilibrium between two systems $`A`$ and $`B`$, such that the maximum of the nonadditive entropy $`\stackrel{~}{S}_q`$ holds for the composite system $`(A+B)`$

$$\stackrel{~}{S}_q(A+B)=\stackrel{~}{S}_q(A)+\stackrel{~}{S}_q(B)+(1q)\stackrel{~}{S}_q(A)\stackrel{~}{S}_q(B),$$
(3)
under the constraints which are fixed values of $`n`$ number of additive quantities

$$E_k(A+B)=E_k(A)+E_k(B),k=1,\mathrm{},n.$$
(4)
As in standard thermodynamics, we define the meta-equilibrium entropy for the nonadditive case to be $`k_\mathrm{B}\stackrel{~}{S}_q(\{E_k|k=1,\mathrm{},n\})`$, where $`k_\mathrm{B}`$ is the Boltzmann’s constant. Apparently, $`\stackrel{~}{S}_q`$ is the explicit entropy function in terms of the set of given additive constraints $`\{E_k|k=1,\mathrm{},n\}`$ and is obtained from the optimum distribution $`\{p_i|i=1,\mathrm{},W\}`$ following from the maximisation of $`\stackrel{~}{S}_q(\{p_i\})`$ .

Now making the variations $`\delta \stackrel{~}{S}_q(A+B)`$ and $`\delta E_k(A+B)`$ for $`k=1,\mathrm{},n`$ vanish for the meta-equilibrium state, we obtain

$$[1+(1q)\stackrel{~}{S}_q(B)]\frac{\stackrel{~}{S}_q(A)}{E_k(A)}=[1+(1q)\stackrel{~}{S}_q(A)]\frac{\stackrel{~}{S}_q(B)}{E_k(B)}.$$
(5)
On rearranging (5), we can achieve separation of variables as

$$\frac{1}{[1+(1q)\stackrel{~}{S}_q(A)]}\frac{\stackrel{~}{S}_q(A)}{E_k(A)}=\frac{1}{[1+(1q)\stackrel{~}{S}_q(B)]}\frac{\stackrel{~}{S}_q(B)}{E_k(B)}=\eta _k.$$
(6)
The parameter $`\eta _k`$, by definition, is common to both the subsystems $`A`$ and $`B`$, and is identified as the intensive variable. Dropping the index $`A`$ or $`B`$, for each subsystem the following equation holds for all values of $`k`$

$$\frac{1}{[1+(1q)\stackrel{~}{S}_q]}\frac{}{E_k}\stackrel{~}{S}_q(E_1,\mathrm{},E_k,\mathrm{},E_n)=\eta _k.$$
(7)
Note that $`\frac{\stackrel{~}{S}_q}{E_k}`$ defines the Lagrange multiplier associated with the maximum entropy method and so according to the above equation, intensive variables and the corresponding Lagrange multipliers are not identical, except when the entropy is additive.

The relation with an additive entropy can be established if we assume that the intensive variable $`\eta _k`$ occuring above is the same as defined for an additive entropy $`S(E_1,\mathrm{},E_k,\mathrm{},E_n)`$ and is given by

$$\eta _k=\frac{S}{E_k}.$$
(8)
We remark here that the choice of the form of additive entropy is arbitrary here; one is free to use either BGS entropy or Renyi entropy. From (7) and (8), this implies that the nonadditive and the additive entropies are related in the following way

$$\frac{1}{(1q)}\mathrm{ln}[1+(1q)\stackrel{~}{S}_q]=S+c,$$
(9)
where the constant $`c`$ has to be independent of $`E_k`$. Thus the transformation (2) of a nonadditive entropy to an additive one is recovered for the particular choice of $`c=0`$.

Next, we explore this relation between nonadditive and additive entropies from the viewpoint of composable entropies as prescribed by the existence of meta-equilibrium. To proceed further, we recall the notion of composability of entropy . An arbitrary entropic form $`\stackrel{~}{S}`$ is defined to be composable, if the total entropy $`\stackrel{~}{S}(A,B)`$ for the composite system can be written as $`\stackrel{~}{S}(A,B)=f[\stackrel{~}{S}(A),\stackrel{~}{S}(B)]`$, where $`f[]`$ is a certain bivariate function of the $`C^2`$ class and is symmetric in its arguments, $`f[\stackrel{~}{S}(A),\stackrel{~}{S}(B)]=f[\stackrel{~}{S}(B),\stackrel{~}{S}(A)]`$.

Suppose the maximum of the general entropy determines a kind of meta-equilibrium between subsystems $`A`$ and $`B`$, under the given additive constraints (4). By equating the variations of the total entropy and total value of the constraint quantity to zero, we get

$$\frac{\stackrel{~}{S}(A,B)}{\stackrel{~}{S}(A)}\frac{\stackrel{~}{S}(A)}{E_k(A)}=\frac{\stackrel{~}{S}(A,B)}{\stackrel{~}{S}(B)}\frac{\stackrel{~}{S}(B)}{E_k(B)}.$$
(10)
To establish an intensive variable common to the two subsystems, a general case can be written in the following form

$`{\displaystyle \frac{\stackrel{~}{S}(A,B)}{\stackrel{~}{S}(A)}}`$ $`=`$ $`{\displaystyle \frac{1}{\lambda }}G[\stackrel{~}{S}(A,B)]{\displaystyle \frac{dh[\stackrel{~}{S}(A)]}{d\stackrel{~}{S}(A)}}h[\stackrel{~}{S}(B)],`$ (11)
$`{\displaystyle \frac{\stackrel{~}{S}(A,B)}{\stackrel{~}{S}(B)}}`$ $`=`$ $`{\displaystyle \frac{1}{\lambda }}G[\stackrel{~}{S}(A,B)]h[\stackrel{~}{S}(A)]{\displaystyle \frac{dh[\stackrel{~}{S}(B)]}{d\stackrel{~}{S}(B)}},`$ (12)
where $`h[]`$ is some differentiable function; $`G[]`$ is also an arbitrary function and $`\lambda `$ is a constant. Therefore, using (11) and (12) in (10) and rearranging, an intensive variable common to systems $`A`$ and $`B`$ can be defined as

$$\frac{1}{\lambda }\frac{1}{h[\stackrel{~}{S}(A)]}\frac{dh[\stackrel{~}{S}(A)]}{d\stackrel{~}{S}(A)}\frac{\stackrel{~}{S}(A)}{E_k(A)}=\frac{1}{\lambda }\frac{1}{h[\stackrel{~}{S}(B)]}\frac{dh[\stackrel{~}{S}(B)]}{d\stackrel{~}{S}(B)}\frac{\stackrel{~}{S}(B)}{E_k(B)}=\eta _k.$$
(13)
Again, $`\frac{\stackrel{~}{S}}{E_k}=\stackrel{~}{\eta }_k`$ is the Lagrange multiplier associated with system $`A`$ or $`B`$. Thus (13) defines the relation between the Lagrange multiplier and the intensive variable for a general composable entropy.

We emphasize that in the present approach the intensive variable associated with a general composable entropy is independent of the function $`G`$, which by definition, is not factorisable into contributions from systems $`A`$ and $`B`$.

For any subsystem ($`A`$ or $`B`$), we can rewrite (13) as

$$\frac{1}{\lambda h[\stackrel{~}{S}]}\frac{h[\stackrel{~}{S}]}{E_k}=\eta _k.$$
(14)
Again invoking the relation (8), we note that the function $`h`$ is related to the additive entropy $`S`$ as $`h[\stackrel{~}{S}]=\mathrm{exp}(\lambda \{S+c\})`$, where $`c`$ has been identified before. Setting $`\stackrel{~}{S}=0`$ for $`S=0`$, we obtain

$$h[0]=\mathrm{exp}(\lambda c).$$
(15)
Note that for $`\lambda 0`$, we have $`h[0]1`$, a property assumed in . Alternately, if the constant is set as $`c=0`$, then $`h[0]`$ is unity for all values of $`\lambda `$. In general, we write

$$h[\stackrel{~}{S}]=h[0]\mathrm{exp}(\lambda S).$$
(16)
This form of the function $`h`$ is solely determined by the requirement of the existence of an intensive variable associated with an additive entropy $`S`$. It is thus intrinsically independent of the form of the $`G[]`$ function.

We remark here that Tsallis type nonadditivity of entropy is obtained as a special case of the above analysis, if we put $`G`$ as a constant equal to unity, and identify for each subsystem

$$h[\stackrel{~}{S}]=1+\lambda \stackrel{~}{S},$$
(17)
which gives $`h[0]=1`$. Then equation (13) is identical to (6) for $`\lambda =(1q)`$. On combining (16) and (17), we have

$$\stackrel{~}{S}=\frac{\mathrm{exp}(\lambda S)1}{\lambda },$$
(18)
which satisfies

$$\stackrel{~}{S}(A+B)=\stackrel{~}{S}(A)+\stackrel{~}{S}(B)+\lambda \stackrel{~}{S}(A)\stackrel{~}{S}(B).$$
(19)
Next, we ascertain the relation between the additive entropy $`S`$ and the composable entropy $`\stackrel{~}{S}`$ for a general function $`G`$. Following , it is convenient to set

$$G[\stackrel{~}{S}(A,B)]=\left(\frac{dH[\stackrel{~}{S}(A,B)]}{d\stackrel{~}{S}(A,B)}\right)^1.$$
(20)
This allows to obtain the composability rule for entropy in terms of the function $`H`$. The case of Tsallis type nonadditivity correpsonds in this notation to $`H`$ as an identity function. For concreteness, let us assume that the function $`G[]`$ is positive, implying that $`H[\stackrel{~}{S}]`$ is a monotonically increasing (differentiable) function of $`\stackrel{~}{S}`$. We give an interpretation of this assumption later on. So using (20), the relations (11) and (12) are rewritten as

$`{\displaystyle \frac{dH[\stackrel{~}{S}(A,B)]}{d\stackrel{~}{S}(A,B)}}{\displaystyle \frac{\stackrel{~}{S}(A,B)}{\stackrel{~}{S}(A)}}`$ $`=`$ $`{\displaystyle \frac{1}{\lambda }}{\displaystyle \frac{dh[\stackrel{~}{S}(A)]}{d\stackrel{~}{S}(A)}}h[\stackrel{~}{S}(B)],`$ (21)
$`{\displaystyle \frac{dH[\stackrel{~}{S}(A,B)]}{d\stackrel{~}{S}(A,B)}}{\displaystyle \frac{\stackrel{~}{S}(A,B)}{\stackrel{~}{S}(B)}}`$ $`=`$ $`{\displaystyle \frac{1}{\lambda }}h[\stackrel{~}{S}(A)]{\displaystyle \frac{dh[\stackrel{~}{S}(B)]}{d\stackrel{~}{S}(B)}},`$ (22)
which further imply

$`{\displaystyle \frac{H[\stackrel{~}{S}(A,B)]}{\stackrel{~}{S}(A)}}`$ $`=`$ $`{\displaystyle \frac{1}{\lambda }}{\displaystyle \frac{dh[\stackrel{~}{S}(A)]}{d\stackrel{~}{S}(A)}}h[\stackrel{~}{S}(B)],`$ (23)
$`{\displaystyle \frac{H[\stackrel{~}{S}(A,B)]}{\stackrel{~}{S}(B)}}`$ $`=`$ $`{\displaystyle \frac{1}{\lambda }}h[\stackrel{~}{S}(A)]{\displaystyle \frac{dh[\stackrel{~}{S}(B)]}{d\stackrel{~}{S}(B)}}.`$ (24)
On integrating (23) or (24), we obtain the composability rule for entropy as

$$H[\stackrel{~}{S}(A,B)]=\frac{1}{\lambda }h[\stackrel{~}{S}(A)]h[\stackrel{~}{S}(B)]+\mathrm{constant}.$$
(25)
The constant is determined by the requirement that $`\stackrel{~}{S}(A,B)=0`$ for $`\stackrel{~}{S}(A)=\stackrel{~}{S}(B)=0`$. Finally on using (16), we obtain

$$H[\stackrel{~}{S}(A,B)]=\frac{h^2[0]}{\lambda }\left\{\mathrm{exp}[\lambda \{\mathrm{S}(\mathrm{A})+\mathrm{S}(\mathrm{B})\}]1\right\}+\mathrm{H}[0].$$
(26)
Using the fact that for an ordered system, say $`B`$, $`S(B)=0`$, which implies $`\stackrel{~}{S}(A,B)=\stackrel{~}{S}(A)`$, we obtain the transformation between an additive entropy and the general composable entropy for a subsystem ($`A`$ or $`B`$) as follows

$$S=\frac{1}{\lambda }\mathrm{ln}\left[1+\frac{\lambda }{h^2[0]}(H[\stackrel{~}{S}]H[0])\right].$$
(27)
This relation may be taken as the most general transformation connecting additive and composable nonadditive entropies compatible with the generalized zeroth law. Note that it is not implied that any composable entropy can be mapped to an additive entropy. As we argue below, the identification of an intensive variable which is independent of the function $`G`$ or $`H`$, seems to imply that the composable entropy is either $`\stackrel{~}{S}`$ which satisfies Tsallis-type nonadditivity, or it is a monotonically increasing function $`H`$ of $`\stackrel{~}{S}`$.

To illustrate this conclusion, suppose that the meta-equilibrium is described by the maximum of $`H[\stackrel{~}{S}(A,B)]`$ under the additive constraints (4). Then making the variation $`\delta H`$ vanish under $`\delta E_k(A)=\delta E_k(B)`$, the condition of equilibrium implies

$$\frac{H[\stackrel{~}{S}(A,B)]}{\stackrel{~}{S}(A)}\frac{\stackrel{~}{S}(A)}{E_k(A)}=\frac{H[\stackrel{~}{S}(A,B)]}{\stackrel{~}{S}(B)}\frac{\stackrel{~}{S}(B)}{E_k(B)}.$$
(28)
This can be rewritten as

$$\frac{dH[\stackrel{~}{S}(A,B)]}{d\stackrel{~}{S}(A,B)}\frac{\stackrel{~}{S}(A,B)}{\stackrel{~}{S}(A)}\frac{\stackrel{~}{S}(A)}{E_k(A)}=\frac{dH[\stackrel{~}{S}(A,B)]}{d\stackrel{~}{S}(A,B)}\frac{\stackrel{~}{S}(A,B)}{\stackrel{~}{S}(B)}\frac{\stackrel{~}{S}(B)}{E_k(B)},$$
(29)
Now assuming the conditions (11), (12) and on using (20) we obtain the same relation as (13). Thus the intensive variable corresponding to the maximisation of $`\stackrel{~}{S}(A,B)`$ is identical to that for the maximisation of $`H[\stackrel{~}{S}(A,B)]`$, which is in accordance with the fact that $`H`$ has been assumed to be a monotonic function of $`\stackrel{~}{S}`$ and thus the equilibrium condition obtained from each is identical under similar constraints. Again note that for $`H`$ to be identical function, the above equations go neatly to the case of the simplest composable entropy.

It may be pointed out that although the intensive variable is the same for the maximisation of either $`\stackrel{~}{S}`$ or $`H[\stackrel{~}{S}]`$, yet the relation between intensive variable $`\eta _k`$ and the corresponding Lagrange multiplier $`\stackrel{~}{\eta }_k`$ depends on the specific entropy chosen for maximisation. For the case of Tsallis entropy, the intensive variable is given by

$$\eta _k=\frac{\stackrel{~}{\eta }_k}{\{1+(1q)S_{q}^{}{}_{}{}^{(\mathrm{T})}\}}.$$
(30)
To obtain the corresponding relation for the case of most general composable entropy, we can apply the partial derivative with respect to $`E_k`$ to Eq. (27) and obtain

$$\eta _k=\frac{\stackrel{~}{\eta }_k^{(H)}}{\{h^2[0]+\lambda (H[\stackrel{~}{S}]H[0])\}}\frac{dH[\stackrel{~}{S}]}{d\stackrel{~}{S}}.$$
(31)
Summarising, it has been shown in literature that an additive entropy underlies the Tsallis entropy, which alongwith an intensive variable, preserves the standard thermodynamic structure. In this letter, we have looked at the relation between nonadditive and additive entropies from the viewpoint of the generalised zeroth law and the notion of composable entropies. We have shown the mapping between Tsallis-type nonadditive entropy (which is the simplest composable form consistent with meta-equilibrium) and an additive entropy by identifying the intensive variable common to two systems in meta-equilibrium. The abovementioned mapping becomes possible if this intensive variable is also the one associated with an additive entropy. Further, we have argued that if the meta-equilibrium is alternately described by the maximum of the most general composable entropy $`H`$, then $`H`$ can be interpreted as a monotonically increasing function of the Tsallis-type nonadditive entropy. The fact that the meta-equilibrium state is expected to be physically the same in either maximisation problem, is consistent with the observation that the intensive variable is same in both the cases.

Tha author is grateful to Professor Gerhard Soff for encouragement and support, and to Professor Sumiyoshi Abe for many useful discussions.


text:
# Competition of Percolation and Phase Separation in a Fluid of Adhesive Hard Spheres

## Abstract

Using a combination of Monte Carlo techniques, we locate the liquid–vapor critical point of adhesive hard spheres. We find that the critical point lies deep inside the gel region of the phase diagram. The (reduced) critical temperature and density are $`\tau _c=0.1133\pm 0.0005`$ and $`\rho _c=0.508\pm 0.01`$. We compare these results with the available theoretical predictions. Using a finite-size scaling analysis, we verify that the critical behavior of the adhesive hard sphere model is consistent with that of the 3D Ising universality class.

The structure of a simple liquid is well described by that of a system of hard spheres at the same effective density. To a good approximation, the effect of attractive interactions on the liquid structure can be ignored. This feature of simple liquids is implicit in the Van der Waals theory for the liquid–vapor transition, and has been made explicit in the highly successful thermodynamic perturbation theories for simple liquids perturbation . The perturbation approach becomes exact as the range of the attractive interaction tends to infinity while its integrated strength remains constant longrange . We refer to this limit as the ‘Van der Waals’ (VDW) limit vdwnote . Conversely, as the attractive forces become shorter-ranged and stronger, the perturbation approach is likely to break down. Fluids with strong, short-ranged attraction (so-called ‘energetic’ fluids energetic ) are of growing importance in the area of complex liquids. For example, short-range attractions are thought to be responsible for the transition from a ‘repulsive’ to an ‘attractive’ glass energetic , which has recently been observed experimentally in PMMA (polymethylmethacrylate) dispersions pmma .

In this Letter, we consider a model system that can be considered as the prototypical energetic fluid: a fluid of adhesive hard spheres (AHS). Introduced in 1968 model , the AHS model is a reference system for particles with short range attractions. The pair potential consists of an impenetrable core plus a surface adhesion term that favors configurations where spheres are in contact. At larger separations, there is no interaction. The AHS model can be considered as the ‘anti Van der Waals’ limit.

Baxter showed model  that the Percus–Yevick (PY) equation can be solved analytically for adhesive hard spheres. In fact, Baxter’s solution is often used to to analyze experimental results for systems as diverse as silica suspensions silica , copolymer micelles micelles , and the fluid phase of lysozyme lysozyme .

One important feature of the AHS model is that its phase diagram contains a liquid–vapor coexistence region pathological . The PY equation offers different routes to estimate the location of the liquid–vapor critical point. However, the ‘compressibility route’ model  and the ‘energy route’ energyroute  lead to estimates for the critical temperature that differ by some 20%, while the estimates for the critical density differ by almost a factor of three. For the analysis of experimental data, it is important to know the location of the critical point more accurately. The reason is that, upon cooling, adhesive hard spheres can percolate to form a gel. Information about the location of this percolation curve is available from simulation method ; lee  and from analytical estimates chiew . It is not clear, though, how percolation interferes with the fluid–fluid phase separation. If Baxter’s estimate of the critical point were right, the critical point would be on, or near, the percolation curve. However, if the estimate of Ref. energyroute  were closer to reality, then the fluid–fluid critical point would be deep inside the gel phase. Clearly, this difference has implications for the possibility to observe the fluid–fluid critical point in systems with short-ranged attraction. Also for the (attractive) glass transition, it makes a considerable difference whether or not the transition line runs close to a liquid–vapor critical point.

In spite of its obvious importance as a reference system for ‘energetic’ (as contrasted to ‘entropic’) complex liquids energetic , there exist, to our knowledge, no accurate numerical estimates of the critical point of the AHS model. This is not surprising, as computer simulations of the low-temperature AHS model are notoriously difficult, due its propensity to form large, even percolating, clusters. In this Letter, we report a Monte Carlo (MC) simulation study that allows us to locate the AHS critical point.

The AHS interaction is derived from a square well potential by taking a limit in which the well becomes infinitesimally narrow at the same time as becoming infinitely deep in such a way that the integrated Boltzmann weight of bound configurations remains finite model . The interaction is most easily defined by the expression for the Boltzmann factor as a function of pair separation, $`r`$:

$$\mathrm{exp}[U(r)/kT]=\mathrm{\Theta }(r\sigma )+\frac{\sigma }{12\tau }\delta (r\sigma ).$$
(1)
In Eq. (1), $`U(r)`$ and $`kT`$ are the formal pair potential and thermal energy, respectively, while $`\sigma `$ is the hard core diameter, which we henceforth use as the unit of length. The step function $`\mathrm{\Theta }`$ accounts for the hard sphere repulsion, and the Dirac delta function introduces the surface adhesion. The parameter $`\tau `$ determines the strength of the adhesion and can be interpreted as an effective temperature, or an ‘inverse stickiness parameter.’

To compare the properties of the AHS and VDW models, we need a measure for the strength of the attractive interactions that is meaningful in both limits. This is best achieved by comparing the reduced second virial coefficients: $`B_2^{}B_2/B_2^{\mathrm{HS}}`$, where $`B_2^{\mathrm{HS}}`$ is the second virial coefficients of hard spheres. For the AHS model, $`B_2^{}=11/(4\tau )`$, and we use this expression to define the $`\tau `$-parameter for the VDW model. In the VDW limit, the free-energy density is given by $`f_{\mathrm{VDW}}(\rho )=f_{\mathrm{HS}}(\rho )a\rho ^2`$ where $`f_{\mathrm{HS}}`$ is the free energy density of a system of hard spheres at density $`\rho `$, while the constant $`a`$ measures the strength of the attractive forces. Hence, in the VDW limit, $`\tau =kTB_2^{\mathrm{HS}}/(4a)`$ is proportional to the temperature. We henceforth refer to $`\tau `$ as the temperature parameter for both the VDW and AHS systems.

In the AHS model, particles may stick and form clusters. At sufficiently low $`\tau `$, clusters may percolate (span the system), mimicking the infinite clusters that form during gelation in the real system. Percolating clusters pose serious problems for any simulation scheme that employs volume changing moves (e.g. constant pressure MC), since all such moves are rejected as soon as percolating clusters appear.

To study the phase behavior of adhesive hard spheres, we therefore used Grand Canonical MC (GCMC) simulation. This technique can be used beyond the percolation threshold. Yet, even with GCMC, equilibration of the system at low temperatures is prohibitively slow. We therefore combine GCMC with parallel tempering to speed up equilibration, and with multiple histogram reweighting, to make optimal use of the available simulation data. All simulation techniques had to be specifically adapted for the AHS model.

The fact that the attractive term in Eq. (1) is infinitesimally narrow means that, in a simulation, the chance of generating a bound configuration by random displacement of a particle is vanishingly small. Conversely, the probability of breaking such a bond, once formed, is also negligible. Consequently, conventional Metropolis sampling cannot be applied to the AHS model.

We therefore employ a modification of the AHS-MC schemes described in Ref. method , where single particle displacements explicitly make and break up to three contacts with the test particle simultaneously. States with higher coordination numbers can be reached indirectly through suitable combinations of moves.

In our GCMC simulations, we employ particle insertion and removal steps with equal probability that together constitute 45% of MC steps. We only insert and remove particles that are not bound to others.

To speed up equilibration, we perform cluster translation moves with probability 5%. A particle is chosen at random, and the cluster to which it belongs is translated by a random amount in each Cartesian direction up to a maximum that is inversely proportional to the number of particles in the cluster.

To overcome the slow equilibration of large clusters at low $`\tau `$ or high density, we have used the parallel tempering scheme of Geyer tempering . Parallel tempering can be performed with replicas at different temperatures, different chemical potentials, or combinations of both. We have found it most convenient to gather statistics across the full density range for one temperature at a time, and therefore chose to run a hierarchy of chemical potentials, $`\mu `$, with the same value of $`\tau `$. Each replica attempts a configuration exchange with one of its neighbors in $`\mu `$ (alternating between the higher and lower neighbor) every 200 MC steps. The acceptance probability for an exchange between replicas $`i`$ and $`j`$ is $`\mathrm{min}[1,(z_i/z_j)^{N_jN_i}]`$, where $`N_i`$ is the number of particles in the current configuration of replica $`i`$, and the activity $`z`$ is related to the chemical potential by $`z=\mathrm{\Lambda }^3\mathrm{exp}(\mu /kT)`$, with $`\mathrm{\Lambda }`$ the thermal de Broglie wavelength.

The final computational tool needed to analyze the data is multiple histogram reweighting multihist . The probability of observing the system in a state with $`N`$ particles and $`B`$ binary contacts between the particles at temperature $`\tau `$ and activity $`z`$ can be exactly decomposed into the form

$$p_{NB}(\tau ,z)=\mathrm{\Omega }_{NB}\tau ^Bz^N/\mathrm{\Xi }(\tau ,z),$$
(2)
where $`\mathrm{\Omega }_{NB}`$ is an effective density of states, and $`\mathrm{\Xi }(\tau ,z)`$ is the grand canonical partition function. Knowledge of $`\mathrm{\Omega }_{NB}`$ up to a multiplicative constant is sufficient to calculate the distribution of $`N`$ or $`B`$ at any set of conditions ($`\tau ,z`$). $`\mathrm{\Omega }_{NB}`$ can be obtained over a wide range of ($`N,B`$) by combining overlapping two-dimensional histograms of $`N`$ and $`B`$ collected from simulations at different $`\tau `$ and $`z`$.

Simulations at sufficiently low $`\tau `$ show a range of activities in which the density distribution is bimodal, indicating the coexistence of a high- and a low-density fluid phase. We studied this coexistence at four different system sizes, labeled by the length $`L`$ of the cubic simulation box. For $`L=5,6`$ and $`8`$, the simulations consisted of $`10^6`$ equilibration blocks and $`10^7`$ acquisition blocks, where a block is $`L^3/\sigma ^3`$ MC trial moves, each of which may be a particle displacement, particle insertion/removal, cluster translation, or replica exchange. At each temperature, seven or eight parallel replicas were used, with the range of activities chosen to span the coexisting densities. At $`L=10`$, the numbers of equilibration and acquisition blocks had to be lowered to $`2.5\times 10^5`$ and $`2.5\times 10^6`$, respectively, to obtain results in a reasonable computational time.

The $`(N,B)`$ histograms for different $`\tau `$ and $`z`$ were combined for each system size separately. The coexisting densities were then obtained as a function of temperature by reweighting the histograms to find the activity at which the two density peaks have equal height. Figure 1 shows the resulting coexistence curves. Strong finite size effects are immediately visible. The curves obtained from different system sizes coincide at sufficiently low temperature, but deviate significantly from each other as the critical region is approached.

The size dependence of the coexistence distributions can be used to locate the critical point more accurately using the approach developed by Bruce and Wilding universality ; ljfluid . This method uses the fact that, within a universality class, the critical distribution of the order parameter is invariant up to a rescaling of the order parameter. As the AHS interactions are short-ranged, the fluid–fluid critical point of this model is expected to belong to the 3D Ising universality class. The critical distribution is known to high precision from studies of lattice systems, and accurate analytic fits are available ising .

Due to the absence of particle–hole symmetry in off-lattice models, the distribution of the density in the AHS model is not quite symmetric about its mean. Symmetry can be restored by accounting for the mixed character of the scaling fields. The appropriate order parameter is then no longer the pure particle density $`\rho =N/L^3`$, but includes a contribution from the energy density, which in the AHS model can be taken as $`u=B/L^3`$. The particle and energy density operators are replaced by the linear combinations universality

$$=\frac{ur\rho }{1sr}=\frac{\rho su}{1sr},$$
(3)
where $`s`$ and $`r`$ are system-dependent field mixing parameters that are identically zero for models with Ising symmetry. Precisely at criticality, the distribution of $``$ in a system of sufficiently large linear length $`L`$ takes on the universal form $`p_L()=p^{}(a_L[_c])`$, where $`_c=`$ evaluated at the critical temperature, and $`a_LL^{\beta /\nu }`$ with $`\beta `$ and $`\nu `$ critical exponents.

Table 1 lists the apparent critical temperature $`\tau _a`$ and mixing parameter $`s`$ at which the distribution of $``$ most closely matches $`p^{}(x)`$ for the four system sizes studied. The fact that $`s`$ is small indicates that the distribution of $`\rho `$ itself is nearly symmetric. The apparent size-dependence of $`s`$ is probably—at least partly—due to the fact that, despite the long simulations employed, statistical fluctuations in the histograms are still significant. For small $`L`$, $`\tau _a`$ is expected to show some size dependence due to finite-size corrections to scaling ljfluid . However, the differences in $`\tau _a`$ shown in Table 1 are not monotonic, and are again likely to be due to limited statistics. Rather than attempting to extrapolate to the infinite system limit, we therefore quote the critical parameters that most closely reproduce the universal critical form, with conservative error bars:

$$\tau _c^{\mathrm{sim}}=0.1133\pm 0.0005\rho _c^{\mathrm{sim}}=0.508\pm 0.01.$$

Figure 2 shows the distribution of the order parameter at the proposed $`\tau _c^{\mathrm{sim}}`$ (not the individual apparent critical temperatures $`\tau _a`$) for all four system sizes. The collapse of the data onto the Ising distribution shows that the transition is consistent with this universality class.

The size dependence of the scaling factor $`a_L`$ in Table 1 permits an estimate of the ratio of the critical exponents $`\beta `$ and $`\nu `$. A straightforward fit to the asymptotic form $`a_LL^{\beta /\nu }`$ yields $`\beta /\nu =0.50\pm 0.04`$, which is clearly compatible with the Ising value $`\beta /\nu =0.52`$.

We can now compare our numerical results with various theoretical predictions. Below, we list the estimates based on the PY compressibility model  and energy energyroute  routes, and the estimate based on the VDW (mean-field) expression meanfield :

$`\tau _c^{\mathrm{PYc}}`$ $`=0.0976`$ $`\rho _c^{\mathrm{PYc}}`$ $`=0.232`$
$`\tau _c^{\mathrm{PYe}}`$ $`=0.1185`$ $`\rho _c^{\mathrm{PYe}}`$ $`=0.609`$
$`\tau ^{\mathrm{VDW}}`$ $`=0.0943`$ $`\rho ^{\mathrm{VDW}}`$ $`=0.250.`$
The coexistence curves for the two PY routes can be obtained numerically using the analytical expressions analytic  for the chemical potential and pressure. These theoretical curves, together with the mean-field result and the simulation data at $`L=8`$ are shown in Fig. 3. The simulation results lie between the two PY predictions, but are clearly closer to those of the energy equation.

Also plotted in Fig. 3 is the percolation threshold. The theoretical result chiew  corresponds to the PY estimate of the points where the cluster size diverges. In simulations, we define the percolation threshold as the locus of points were the probability of observing a percolating cluster in a canonical simulation is 50%. This definition is relatively size independent, and our results with $`N=500`$ particles are consistent with the more elaborate analysis of Lee lee . Both the theoretical and simulated percolation lines clearly show that the critical point of the fluid–fluid transition lies well within the percolated part of the phase diagram.

The phase diagram reported in Fig. 3 is likely to be representative of that of real colloidal systems with short-ranged attraction at the same value of $`\tau `$. If so, our findings suggest that the attractive glass transition pmma ; attractive  is relatively far removed from the fluid–fluid critical point.

The work of the FOM Institute is part of the research program of FOM and is made possible by financial support from the Netherlands organization for Scientific Research (NWO). We thank Dr. R. P. Sear for the mean-field data in Fig. 3, and gratefully acknowledge insightful discussions with Prof. W. C. K. Poon and Dr. N. Kern.


text:
# Eternal expansion of closed universe

## 1 Introduction

As it is known, data on distance modulus versus redshift reveal the presence of an acceleration of the cosmological expansion . This is possible in existence of a cosmological repulsion which is described usually by $`\mathrm{\Lambda }`$-term or equivalently by the presence of an vacuum-like medium with a negative pressure. An open universe will expand forever even if $`\mathrm{\Lambda }=0`$. A closed universe at $`\mathrm{\Lambda }=0`$ will stop expanding in future and begin to recollapse without fail. If $`\mathrm{\Lambda }>0`$, the expansion becomes accelerated when the doubled vacuum density exceeds the decreasing matter density. The open universe is sure to pass this condition and its eternal expansion has not an alternative111In the work  the possibilities of recollapse of a flat universe owing to decay of cosmological constant, collision with a null singularity and formation of space-like curvature singularity during expansion have been investigated. In this work these possibilities are not considered.. Reaching of the closed universe the acceleration state and its consequent eternal expansion requires a potential barrier overcoming which is possible under certain conditions. In this work a restriction on the closed universe parameters being necessary to getting over the barrier and the universe dynamics in this case are investigated.

## 2 Condition of eternal
expansion

The closed universe dynamics is described by the Einstein field equation for curvature radius $`a`$ of 3-space in comoving with cosmological expansion nonrotating (synchronous) frame:

$$\dot{a}^2=\frac{A_m}{a}+\frac{a^2}{A_v^2}1.$$
(1)
Here and below convention $`c=1`$ is used. The evolution constants $`A_m`$ and $`A_v`$ following  are named the Friedmann integrals. The matter integral

$$A_m=\frac{8\pi }{3}G\rho _ma^3,$$
(2)
where $`\rho _m`$ is the dust-like matter density including a dark matter, $`G`$ is the gravitational constant. The vacuum integral $`A_v`$

$$A_v^2=\frac{8\pi }{3}G\rho _v=\frac{\mathrm{\Lambda }}{3},$$
(3)
where $`\rho _v`$ is the vacuum energy density, $`\mathrm{\Lambda }`$ is the cosmological constant. The equation (1) has been written without taking into account a radiation which contribution dominates at an early stage of the universe evolution and is negligible at the stage considered. Define

$$U_m=\frac{A_m}{a},U_v=\frac{a^2}{A_v^2},U=U_m+U_v$$

and write the equation (1) in form of an ”energy conservation law”

$$\dot{a}^2+U(a)=1.$$

Graphs of the functions $`U_m(a)`$, $`U_v(a)`$, $`U(a)`$ are plotted in Fig. 1.

At the point $`a=\stackrel{~}{a}`$,

$$\stackrel{~}{a}^3=\frac{1}{2}A_mA_v^2$$
(4)
the maximum $`\stackrel{~}{U}`$ of the potential barrier $`U(a)`$ is disposed,

$$\stackrel{~}{U}=U(\stackrel{~}{a})=\left(\frac{\alpha A_m}{A_v}\right)^{2/3},$$
(5)
where $`\alpha =3\sqrt{3}/2`$. At this point the balance of matter gravitation and vacuum antigravitation is achieved, $`\rho _m=2\rho _v`$. The levels $`\stackrel{~}{U}`$, $`\stackrel{~}{U}_m`$, $`\stackrel{~}{U}_v`$, $`0`$ are equidistant:

$$\stackrel{~}{U}\stackrel{~}{U}_m=\stackrel{~}{U}_m\stackrel{~}{U}_v=\stackrel{~}{U}_v=\frac{1}{3}\stackrel{~}{U}.$$

At the point $`a=\stackrel{~}{a}`$ the deceleration parameter changes the sign:

$$q=\frac{\ddot{a}a}{\dot{a}^2}=\frac{\stackrel{~}{a}^3a^3}{a^3aA_v^2+2\stackrel{~}{a}^3}.$$
(6)
At $`a<\stackrel{~}{a}`$ the matter gravitation dominates and $`q>0`$, at $`a>\stackrel{~}{a}`$ the vacuum antigravitation dominate and $`q<0`$. At the point $`a=a^{}`$,

$$a^3=2\stackrel{~}{a}^3$$

$`U_m^{}=U_v^{}`$ and $`\rho _m=\rho _v`$. If the level $`U=1>U_m^{}=U_v^{}`$ then $`A_m>A_v`$ (as represented in Fig. 1); $`A_mA_v`$ otherwise. The Friedmann integral $`A_m`$ is a maximum value of the curvature radius in the standard model $`\mathrm{\Lambda }=0`$, the integral $`A_v`$ equals an initial value $`a`$ in the de Sitter universe $`\rho _m=0`$:

$$a(t)=A_v\mathrm{cosh}(t/A_v).$$
(7)
A closed universe will expand eternally if the level $`U=1`$ is above the potential barrier summit $`U=\stackrel{~}{U}`$ (5), $`|\stackrel{~}{U}|1`$ (the equality sign corresponds to the asymptotic approaching $`a`$ to $`\stackrel{~}{a}`$ (4)):

$$\frac{\alpha A_m}{A_v}1.$$
(8)
The condition of potential barrier overcoming corresponds to the strict inequality. Thus, universe fate is determined by the Friedmann integrals ratio. The condition (8) means the limitation on the cosmological constant

$$\mathrm{\Lambda }\frac{4}{3}A_m^2.$$
$`(8^{})`$
The $`A_m`$ value (2) can be estimated very rough due to a considerable uncertainty of the scale factor $`a`$. According to  $`A_vA_m`$ and even some more, and then a fulfilment of the condition (8) turns out to be problematic.

Find out the condition (8) for the density parameters $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_v`$,

$`\mathrm{\Omega }_m`$ $`=`$ $`{\displaystyle \frac{\rho _m}{\rho _c}}=\left(1+\left({\displaystyle \frac{a}{a^{}}}\right)^3{\displaystyle \frac{a}{a^{}}}\left({\displaystyle \frac{A_v}{A_m}}\right)^{\frac{2}{3}}\right)^1,`$ (9)
$`\mathrm{\Omega }_v`$ $`=`$ $`{\displaystyle \frac{\rho _v}{\rho _c}}=\left({\displaystyle \frac{a}{a^{}}}\right)^3\mathrm{\Omega }_m=(HA_v)^2,`$ (10)
where $`\rho _c`$ is the critical density, $`\rho _c=3H^2/8\pi G`$, $`H=\dot{a}/a`$ is the Hubble constant. The density parameter $`\mathrm{\Omega }_m`$ achieves a maximum at $`a=a_1`$

$$a_1=\frac{A_v}{\sqrt{3}};$$

$`\mathrm{\Omega }_m=1`$ at $`a=0`$ and $`a=A_v`$. The parameter $`\mathrm{\Omega }_v`$ achieves a maximum (and the Hubble constant $`H`$ \- a minimum) at $`a=a_2`$

$$a_2=\frac{3}{2}A_m;$$

$`\mathrm{\Omega }_v=1`$ at $`a^1=0`$ and $`a=A_m`$. The deceleration parameter (6)

$$q=\frac{1}{2}\mathrm{\Omega }_m\mathrm{\Omega }_v;$$
(11)
$`q=0.5`$ at $`a=0`$ and $`a=a_1`$, $`q=1`$ at $`a^1=0`$ and $`a=a_2`$. The curvature radii $`a_1`$, $`a_2`$ and $`\stackrel{~}{a}`$ satisfy the ratios

$$\frac{a_2}{a_1}=\left(\frac{a_2}{\stackrel{~}{a}}\right)^{3/2}=\left(\frac{\stackrel{~}{a}}{a_1}\right)^3=\frac{\alpha A_m}{A_v}.$$

The Einstein field equation (1) connects the density parameters and the Friedmann integrals by the correlation

$$\left(\mathrm{\Omega }_m+\mathrm{\Omega }_v1\right)^3=\left(\frac{A_v}{A_m}\right)^2\mathrm{\Omega }_m^2\mathrm{\Omega }_v.$$
(12)
The condition of the eternal expansion of the closed universe (8) means

$$\left(\mathrm{\Omega }_m+\mathrm{\Omega }_v1\right)^3\alpha ^2\mathrm{\Omega }_m^2\mathrm{\Omega }_v.$$
(13)
If in accordance with the last data  we assume $`\mathrm{\Omega }0.03`$, where $`\mathrm{\Omega }=\mathrm{\Omega }_m+\mathrm{\Omega }_v`$, then the equation (12) yields $`A_v10^2A_m`$ and the condition (8), (13) is fulfilled easily.

As it is clear from (12), the ratio $`A_v/A_m`$ is connected with the degree of the universe spatial flatness, and $`A_v/A_m0`$ when $`\mathrm{\Omega }1`$. Therefore the problem of the cosmic coincidence of the Friedmann integrals $`A_m`$ and $`A_v`$ discussed in the works  seems to be not actual. Do not coincide the expansion factor $`a`$ today with constant $`A_v`$ either:

$$\left(\frac{a}{A_v}\right)^3=\frac{\mathrm{\Omega }_v}{\mathrm{\Omega }_m}\frac{A_m}{A_v}$$
(14)
and $`a/A_v\mathrm{}`$ when $`\mathrm{\Omega }1`$. There exists the problem of coincidence of $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_v`$ that should be explained by anthropic arguments .

## 3 Phase curves

Curves in Fig. 2 represents the dynamics of universes according to (12).

The point $`(1,0)`$ corresponds to a Big Bang, the point $`(0,1)`$ signifies a de Sitter universe. The evolution of an universe with the fixed Friedmann integrals ratio $`A_m/A_v`$ is described by a the curve in $`(\mathrm{\Omega }_m,\mathrm{\Omega }_v)`$ plane. One and only one curve passes through each of the point $`(\mathrm{\Omega }_m,\mathrm{\Omega }_v)`$, excepting $`(0,1)`$ and $`(1,0)`$. It means an observational determination of the values $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_v`$ fixes the universe evolution in the model circumscribed (1), (12). The phase-space area $`(\mathrm{\Omega }_m,\mathrm{\Omega }_v)`$ satisfying the condition (8), (13) is disposed between the curves 1 and $`1^{}`$ which correspond to the equality $`A_v=\alpha A_m`$. All the curves in this area corresponding the universes expansion begin at the point $`(1,0)`$; the evolutions are completed by the de Sitter universe $`(0,1)`$. On the curves 2 and 3 the correlations $`A_v=2A_m`$ and $`A_v=A_m`$ are fulfilled respectively. Below the curve 1 universes expand from the Big Bang $`(1,0)`$ till a maximum value $`a`$ defined left point of intersection of the straight line $`U=1`$ and the curve $`U(a)`$ when $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_v`$ go to infinity, afterwards the universes begin to recontract and recollapse $`(1,0)`$. Above the curve $`1^{}`$ expansion goes on the other (large $`a`$) side of the potential barrier without the Big Bang from a minimum value $`a`$ defined the right intersection point of the lines $`U=1`$ and $`U(a)`$ to the de Sitter metric $`(0,1)`$. A passage from $`(1,0)`$ to $`(0,1)`$ on a straight line corresponds a flat universe $`\mathrm{\Omega }_m+\mathrm{\Omega }_v=1`$. The dependence $`a(t)`$ (7) at $`t\mathrm{}`$ becomes exponential

$$a(t)=A_v\mathrm{exp}(t/A_v),$$

which describes the final fate of any expanding universe with the $`(\mathrm{\Omega }_m,\mathrm{\Omega }_v)`$ values above the curve 1.

Plots of the parameters $`\mathrm{\Omega }_m`$ (9), $`\mathrm{\Omega }_v`$ (10), $`q`$ (11) versus $`a/\stackrel{~}{a}`$ for the universe $`A_v=2A_m`$ (the curve 2 in Fig. 2) and for the almost flat universe $`A_v=10^2A_m`$ are given in Fig. 3(a) and 3(b) respectively.

It should be mentioned that at an early stage of an universe evolution a radiation density $`\rho _r`$ dominates, then phase-space is three-dimensional $`(\mathrm{\Omega }_m,\mathrm{\Omega }_v,\mathrm{\Omega }_r)`$. The phase curves begin at the point $`(0,0,1)`$; they are directed to the $`\mathrm{\Omega }_r`$ axis for a closed universe, and when the matter grows predominant over the radiation, the curves pass on $`(\mathrm{\Omega }_m,\mathrm{\Omega }_v)`$ plane as it is shown in Fig. 2. The vacuum energy density $`\rho _v`$ is to prevail at a still more early stage preceding radiation-dominated, and then the universe develops from de Sitter metric to de Sitter metric along a closed phase curve. This curve simplified has been plotted here for a flat universe $`\mathrm{\Omega }_m+\mathrm{\Omega }_v+\mathrm{\Omega }_r=1`$.

## 4 Summary

The measured cosmological expansion acceleration means that our Universe will expand eternally. The condition of the closed universe eternal expansion restricts the Friedmann integrals ratio by the correlation (8). Written down for the density parameters $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_v`$ this correlation (13) separate the area in the phase-space (Fig. 2). The measured parameters $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_v`$ of our Universe are found in this area. The Friedmann integrals of our Universe do not coincide. The phase curves of the eternally expanding closed universes have been plotted. The dependence of the cosmological parameters $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_v`$, $`q`$, $`H`$ on the universe curvature radius $`a`$ has been investigated.


text:
# Critical dynamics of the simple-cubic Heisenberg antiferromagnet RbMnF3: Extrapolation to q=0

## I Introduction

The dynamic critical behavior of RbMnF3, a close realization of the isotropic, simple cubic Heisenberg antiferromagnet, has been the subject of several experimental (see Coldea et alcoldea  and references therein), and theoretical studiesmazenkoTc ; mazenkoTb ; cuccoli ; hohenhal ; halhoh69 . The real time dynamics of this system is governed by coupled equations of motion for the magnetic ions. In the classification of Hohenberg and Halperinhohenhal , the critical dynamics of this system pertains to class G, for which the order parameter (the staggered magnetization) is not conserved.

Dynamic critical behavior can be characterized by a dynamic critical exponent $`z`$, and for RbMnF3 the most precise experimental estimatecoldea  is $`z=1.43\pm 0.04`$. This estimate is slightly below the predicted valuemazenkoTc ; cuccoli ; hohenhal ; halhoh69  of $`z=1.5`$ for an isotropic three-dimensional Heisenberg antiferromagnet. Our previous estimate using spin-dynamics simulationsd1  (obtained before we knew the latest experimental resultcoldea ) was $`z=1.43\pm 0.03`$ in good agreement with the experimental value. The slight disagreement between theory and both the experiment and the simulation was perplexing, but one possible explanation was that neither of these latter studies probed the asymptotic critical region sufficiently far. To check this, neutron scattering experiments would have to be performed with smaller momentum transfer $`q`$, a task that can be quite challenging. We can also resort to spin dynamics simulations of the model on larger lattice sizes which allow access to smaller $`q`$ values.

Spin dynamics simulationsDPLMKreview  have proven to be an effective tool to study dynamic behavior of magnetic systems. Direct and quantitative comparisons of magnetic excitation dispersion curves and dynamic structure factor line shapes from experimentcoldea  and spin dynamics simulations have shown good agreement, with no adjustable parameterssd1 ; sd2 . At the critical temperature $`T_c`$, both experiment and simulation find that the average dynamic structure factor $`S(𝐪,\omega )`$, for momentum $`𝐪`$ and frequency $`\omega `$, has a spin-wave peak and an additional central peak not predicted by either the renormalization-group theorymazenkoTc  or mode-coupling theorycuccoli . In a polarized neutron scattering experimentcox  below $`T_c`$, the quasi-elastic peak has been shown to be longitudinal in character. It is thus interesting to investigate whether the central peak at $`T_c`$ originates in the longitudinal, in the transverse or in both components with respect to the staggered magnetization. Another motivation for separating the longitudinal and transverse components of the dynamic structure factor is to test the theoretical predictionmazenkoTb  that below $`T_c`$ the longitudinal momentum-dependent susceptibility behaves as $`\chi ^L(q)1/q`$, whereas the $`q`$-dependence of the transverse component is $`\chi ^T(q)1/q^2`$. Although experimental measurementscoldea  are consistent with $`\chi ^L(q)1/q`$, the lack of reliable small wave vector measurements hindered a conclusive experimental test of this predicted divergence. The experimental value for the transverse component is $`\chi ^T(q)1/q^{1.91\pm 0.05}`$.

In this paper we use spin dynamics simulations to study the dynamic critical behavior of the isotropic Heisenberg antiferromagnet on the simple cubic lattice, using larger lattices than in our previous simulations, as motivated above. We investigate the longitudinal and transverse components of the dynamic structure factor, and compare the $`q`$-dependence of the susceptibility with the predicted forms both at and below $`T_c`$. For brevity, we do not present our methodology in great detail here; it has been described in earlier worksd1 ; DPLMKreview .

## II Model and Methods

We consider three-dimensional classical spins $`𝐒_𝐫`$ of unit length, defined on the sites $`𝐫`$ of $`L\times L\times L`$ simple cubic lattices. The interaction is described by a model Hamiltonian written as

$$=J\underset{<\mathrm{𝐫𝐫}^{}>}{}𝐒_𝐫𝐒_𝐫^{},$$
(1)
where the summation is over pairs of nearest-neighbor spins, and $`J>0`$ is an antiferromagnetic exchange coupling. An earlier high-resolution Monte Carlo simulationkunTc  determined $`T_c=1.442929(77)J`$, so any uncertainty in the location of the critical temperature of this model is negligible.

The dynamics of the spins is governed by coupled equations of motionDPLMKreview  and the dynamic structure factor $`S^k(𝐪,\omega )`$ is given by the Fourier transform of the space- and time-displaced correlation functions $`C^k(𝐫𝐫^{},t)=S_𝐫^k(t)S_𝐫^{}^k(0)S_𝐫^k(t)S_𝐫^{}^k(0)`$, where $`k`$ represents the spin components. The coupled equations of motion were integrated using an algorithmkrech  based on 4th-order Suzuki-Trotter decompositions of exponential operators, with a time step $`dt`$. The integrations were performed up to time $`t_{max}`$, starting from equilibrium configurations at temperature $`T`$ generated with a hybrid Monte Carlo methodsd1 . The correlations $`C^k(𝐫𝐫^{},t)`$ were computed for time displacements ranging from $`0`$ to $`t_{cutoff}`$ and the canonical ensemble was established by averaging results from $`N`$ different initial equilibrium configurations. We used periodic boundary conditions in space and thus we could only access a set of discrete values of momentum transfer given by $`q=2\pi n_q/L`$, where $`n_q=1,2,\mathrm{},L/2`$. Because of computer memory limitations we restricted our data to $`𝐪=(q,0,0)`$, $`(q,q,0)`$, and $`(q,q,q)`$, which correspond to the , , and  directions, respectively.

In an attempt to probe the true asymptotic critical region in momentum and thus provide a more accurate estimate of the dynamic critical exponent, we extended our Monte CarloMCbook  and spin-dynamics simulationsDPLMKreview ; sd1  to systems as large as $`L=72`$ at $`T_c`$. For $`L=72`$, each equilibrium configuration was generated with 2500 hybrid Monte Carlo steps, each of which consists of two Metropolis sweeps through the lattice and eight overrelaxation steps. We used $`dt=0.2/J`$, $`t_{max}=1080/J`$, $`t_{cutoff}=1000/J`$, and $`N=1000`$. We have also improved the statistics of our previous simulations for $`L=48`$ and $`60`$ by increasing the number of initial configurations to $`N=1000`$. With $`L=72`$, the smallest wave vector that we can access is $`q=2\pi /72=2\pi (0.0139)`$, whereas, in our units, the smallest $`q`$ value probed by experimentcoldea  was $`q=2\pi (0.02)`$.

The dynamic critical exponent $`z`$ can be extracted using dynamic finite-size scalingDPLMKreview , which yields $`\omega \overline{S}^k(𝐪,\omega )/\overline{\chi }^k(𝐪)=G(\omega L^z,qL,\delta _\omega L^z)`$ and $`\overline{\omega }_m^k=L^z\overline{\mathrm{\Omega }}^k(qL,\delta _\omega L^z)`$, where $`k`$ represents the polarization, $`\overline{S}^k(𝐪,\omega )`$ is the dynamic structure factor convoluted with a Gaussian resolution function with parameter $`\delta _\omega `$, $`\overline{\chi }^k(𝐪)`$ is the total integrated intensity, and $`\overline{\omega }_m^k`$ is a characteristic frequency, defined by $`_{\overline{\omega }_m^k}^{\overline{\omega }_m^k}\overline{S}^k(𝐪,\omega )𝑑\omega /2\pi =\overline{\chi }^k(𝐪)/2`$. If we set $`\delta _\omega =0`$ the exponent $`z`$ can be obtained from the slope of a graph of ln$`(\omega _m^k)`$ vs ln$`(L)`$, at fixed $`qL`$. To test the robustness of the estimate, we have also used $`\delta _\omega =0.005(72/L)^z`$, and determined $`z`$ iterativelyDPLMKreview ; sd1 . As in our previous worksd1 ; sd2 , the dynamic critical exponent is estimated using the average dynamic structure factor.

After understanding the difference in the dynamic critical exponent between experiment and theory, we can examine other unresolved issues, such as the nature of the central peak at $`T_c`$ observed in experiments, but absent in theoretical studies. To this end, we analyzed longitudinal and transverse motions of the spins with respect to the staggered magnetization. The dynamics of the isotropic Heisenberg model, given by the coupled equations of motion, conserves both the total energy and the uniform magnetization, which is the order parameter for the Heisenberg ferromagnet. However, for the antiferromagnet the order parameter (staggered magnetization) is not a constant of the motion; therefore, separating components of the spin parallel (longitudinal component) and perpendicular (transverse component) to the order parameter is challenging. Our approach to determine the individual components of the spin motion, and thus of $`S(𝐪,\omega )`$, was to rotate the frame of reference to align the $`z`$-axis parallel to the staggered magnetization before starting the time integrations, to make the $`z`$-axis coincide with the longitudinal direction. As we integrated the equations of motion, the direction of the staggered magnetization changed slightly because it is not a conserved quantity. Therefore, after each integration step we re-rotated the frame of reference to re-align the $`z`$-axis with the staggered magnetization, thereby restoring the $`z`$-axis as the longitudinal direction. In this part of our simulations we used $`L=24`$ with $`t_{max}=480/J`$, $`t_{cutoff}=400/J`$, $`dt=0.2/J`$, and $`N=12,000`$, in addition to $`L=36`$, $`48`$, and $`60`$, with $`t_{max}=880/J`$, $`t_{cutoff}=800/J`$, $`dt=0.2/J`$, and $`N=11,000`$, $`7,000`$, and $`3,000`$, respectively. We denote the longitudinal and transverse components of $`S(𝐪,\omega )`$ as $`S^L(𝐪,\omega )`$ and $`S^T(𝐪,\omega )`$, respectively.

We have also investigated $`\chi ^L(q)`$ and $`\chi ^T(q)`$, the integrated intensities of $`S^L(𝐪,\omega )`$ and $`S^T(𝐪,\omega )`$, respectively. When the decay of the dynamic structure factor is slow, as is the case of $`S^T(𝐪,\omega )`$ at $`T_c`$, computing the integrated intensity is complicated by the fact that at high frequencies \[$`\pi /(2dt)\omega \pi /dt`$\] the approximation of the integral Fourier transform as a discrete sum is highly dependent on the summation method (e.g. direct sum, trapezoidal rule, Simpson’s rule, etc.). In terms of cpu time usage, it is more efficient to integrate the equations of motion with the largest time step $`dt`$ that still yields a stable method and accurate time evolutions; however, the maximum accessible frequency is inversely proportional to $`dt`$. Simulations for $`L=24`$ at $`T_c`$ using a smaller time step ($`dt=0.1/J`$) showed that both the Simpson’s rule and the trapezoidal rule produce consistent results for $`S^L(𝐪,\omega )`$ and $`S^T(𝐪,\omega )`$ up to $`\omega 15J`$. A direct comparison between the dynamic structure factor obtained with $`dt=0.1/J`$ and that obtained with the Simpson’s rule and the trapezoidal rule using $`dt=0.2/J`$ up to $`\omega 15J`$ shows that the trapezoidal rule gives a better approximation for $`S^L(𝐪,\omega )`$ and $`S^T(𝐪,\omega )`$ at high frequencies. Estimates of $`\chi ^L(q)`$ and $`\chi ^T(q)`$ are not significantly dependent on the methods used to integrate $`S^L(𝐪,\omega )`$ and $`S^T(𝐪,\omega )`$, respectively, and we use Simpson’s rule for these integrations.

Our approach to determine $`\chi ^T(q)`$ at $`T_c`$ was to compute $`S^T(𝐪,\omega )`$ with the trapezoidal rule and then to integrate it from $`\omega =0`$ to $`8J`$. Although the intensity of $`S^T(𝐪,\omega )`$ has not decayed to zero at $`\omega =8J`$, it is a small fraction of the intensity of the spin-wave peak. Data for $`\omega =4J`$ to $`8J`$ were fitted with an exponential function, which was then used in the integration to $`\omega =\mathrm{}`$. (An exponential function was chosen for this fitting because it is the simplest one that yielded a good fitting in the range of frequency used.) To estimate the longitudinal component $`\chi ^L(q)`$ at $`T_c`$ we simply used the trapezoidal rule to determine $`S^L(𝐪,\omega )`$ and then integrated it from $`\omega =0`$ to $`\omega =10J`$, without further high-frequency corrections, because $`S^L(𝐪,\omega )`$ decays quickly in frequency.

The behavior of $`\chi ^T(q)`$ and $`\chi ^L(q)`$ below $`T_c`$ is studied with spin dynamics simulations at $`T=0.5T_c`$ using $`dt=0.2/J`$ and $`L=24`$ and $`36`$, with $`N=1000`$ for each lattice size. At this temperature both $`S^L(𝐪,\omega )`$ and $`S^T(𝐪,\omega )`$ decay quickly, hence these line shapes were obtained with the trapezoidal rule and they were then integrated from $`\omega =0`$ to $`\omega =10J`$.

For comparison, we have also computed the longitudinal and transverse components of $`S(𝐪,\omega )`$ with respect to the residual uniform magnetization. In this case, the rotation of the frame of reference to align one axis parallel to the uniform magnetization was only required once, before starting the time integrations, because the uniform magnetization vector is conserved.

## III Results

Our results for the average dynamic structure factor $`S(𝐪,\omega )`$ show a spin wave and a central peak for $`TT_c`$. The spin wave dispersion curve for $`L=72`$, $`T=T_c`$ and $`𝐪`$ vectors in the  direction, shown in Fig.1, was obtained by fitting $`S(𝐪,\omega )`$ with Lorentzian spin-wave creation and annihilation peaks centered at $`\omega =\pm \omega _s`$, and a Lorentzian peak at $`\omega =0`$.

We have obtained reasonable fittings for $`q`$ values corresponding to $`n_q=1`$ to $`7`$; for larger $`n_q`$, the spin-wave line shapes were not Lorentzian and the positions of the spin-wave peaks $`\omega _s`$ were read off directly without any fitting. We remind the reader that the first Brillouin zone edge in the  direction occurs at $`n_q=L/2`$. An estimate of the dynamic critical exponent $`z`$ is obtained by fitting the dispersion curve with a functionhohenhal  $`\omega _s=Dn_q^z`$. Since the asymptotic critical region corresponds to small values of $`q`$, we have used $`n_q=1`$ to $`9`$ in the fitting (dashed line in Fig.1), and obtained $`z=1.35\pm 0.05`$. If $`n_q=1`$ is excluded due to its large finite-size effectsd1 , the fitting yields $`z=1.45\pm 0.07`$ (solid line in Fig.1).

A better approach to determine $`z`$ is to use the dynamic finite-size scaling theory outlined above. Using this method, we obtained $`z`$ for different values of $`n_q`$, with no resolution function. Such estimates are denoted as $`z_q`$ and are shown in Fig.2.

In our previous worksd1 , we estimated $`z`$ as the average value obtained using $`n_q=1`$ and $`n_q=2`$, with a maximum lattice size of $`L=60`$. We now have larger $`L`$ and better statistics, and hence the present $`n_q=1`$ and $`n_q=2`$ correspond to $`q`$-values that are closer to the asymptotic critical region. Nevertheless, Fig.2 shows that these values of $`q`$ are not yet in the asymptotic critical region, and a better estimate of $`z`$ is given by $`z_0`$, obtained by extrapolating $`z_q`$ to the limit $`q0`$. We fitted $`z_q`$ with the function $`z_q=z_0+an_q+bn_q^2`$, where $`z_0`$, $`a`$, and $`b`$ are fitting parameters. Using $`n_q=1,2,\mathrm{},7`$ in the fitting (dashed line in Fig.2) we find $`z_0=1.48\pm 0.02`$ and excluding $`n_q=1`$, due to its large finite-size dependence, (solid line in Fig.2) we find $`z_0=1.50\pm 0.02`$. Both estimates agree with the theoretical prediction of $`z=1.5`$. Dynamic finite-size scaling theory with a small resolution function was used to determine $`z_q`$ iteratively. The results thus obtained are within a one-$`\sigma `$ error bar of the respective $`\delta _\omega =0`$ estimates.

The longitudinal and transverse components of the dynamic structure factor with respect to the staggered magnetization are shown in Figs.3a and 3b, respectively, for $`L=36`$, $`T=0.5T_c`$ and $`𝐪`$ in the  direction with $`n_q=3`$.

While the transverse component has a pronounced single spin-wave excitation at $`\omega /J1.49`$ and intensity $`30`$, as shown in the inset in Fig.3b, the structures on $`S^L(𝐪,\omega )`$ have much smaller amplitudes and comprise a central peak, and a series of two-spin-wave subtraction (a-i) and addition (k-r) peaks. Denoting the momentum and frequency of two single spin waves as $`(𝐪_1,\omega _1)`$ and $`(𝐪_2,\omega _2)`$, the excitations resulting from their addition and subtraction have frequencies $`\omega _+=\omega _1+\omega _2`$ and $`\omega _{}=|\omega _1\omega _2|`$, respectively, and momentum $`𝐪=𝐪_1+𝐪_2`$. For odd values of $`n_q`$ there are no two spin-wave peaks at $`\omega =0`$ so the central peak seen in $`S^L(𝐪,\omega )`$ is presumably due to spin diffusion.

At $`T_c`$, $`S^L(𝐪,\omega )`$ (see Fig.4a) seems to be predominantly diffusive with a central peak that is much more intense than the spin-wave peak in $`S^T(𝐪,\omega )`$ (Fig.4b). These data provide clear evidence that the central peak in the average $`S(𝐪,\omega )`$ at $`T_c`$, seen in experiment and previous simulations, but not present in renormalization-group and mode-coupling theories, appears only in the longitudinal component. Two spin-wave excitations at $`T_c`$ cannot be resolved.

Figure 5 shows log-log plots of $`\chi ^L(q)`$ and $`\chi ^T(q)`$, the integrated intensities of $`S^L(q,\omega )`$ and $`S^T(q,\omega )`$, respectively, as a function of momentum in the  direction. The momentum dependence of the integrated intensity has the form $`q^\mathrm{x}`$. At $`T=0.5T_c`$ (Fig.5a), a linear fitting in the log-log plane of $`\chi ^L(q)`$ and $`\chi ^T(q)`$ versus $`n_q`$ for $`L=36`$ using $`n_q=1`$ to $`6`$ gives $`\chi ^L1/q^{0.80\pm 0.16}`$ and $`\chi ^T1/q^{1.94\pm 0.06}`$, whereas if the $`n_q=1`$ point is dropped the linear fitting yields (solid lines) $`\chi ^L1/q^{0.88\pm 0.35}`$ and $`\chi ^T1/q^{1.92\pm 0.10}`$. These results are in agreement with both renormalization-group theory predictionmazenkoTb  and experimental resultscoldea . We have also tried to fit $`\chi ^L(q)`$ with the Ornstein-Zernike mean-field formulacollinsbook  $`\chi ^L(q)=\chi ^L(0)\kappa ^2/(q^2+\kappa ^2)`$, using $`\chi ^L(0)`$ and $`\kappa `$ as fitting parameters; however, this expression did not yield a good fitting.

At $`T_c`$ (Fig.5b), our data for $`L=60`$ indicate that $`\chi ^L1/q^{1.81\pm 0.03}`$ and $`\chi ^T1/q^{1.92\pm 0.04}`$, where $`n_q=2`$ to $`10`$ have been used in the fitting. Finite-size effects on the low-$`q`$ divergence exponents of $`\chi ^L`$ and $`\chi ^T`$ at $`T_c`$ are shown in Fig.6. A linear fitting of these exponents as a function of $`1/L`$ yields $`\chi ^L1/q^{1.88\pm 0.05}`$ and $`\chi ^T1/q^{1.95\pm 0.07}`$ for the thermodynamic limit where $`L=\mathrm{}`$. The dynamic scaling prediction for the static susceptibility at $`T_c`$ ishalhoh69 ; RitFis72  $`\chi =1/q^{2\eta }`$, where for the purpose of this comparison we can useeta\_ref  $`\eta 0.04\pm 0.01`$. We see that while $`\chi ^T`$ is consistent with the dynamic scaling prediction, our large error bars do not exclude the mean-field behavior of $`\chi 1/q^2`$. Our estimate for the divergence of $`\chi ^L`$ at small $`q`$ is slightly less rapid than predicted, but still consistent with it within a two-$`\sigma `$ error bar.

Figure 7 shows a log-log plot of the longitudinal and transverse components of the characteristic frequency as a function of $`L`$, for $`n_q=2`$ and $`\delta _\omega =0`$ at $`T_c`$. These components are denoted as $`\omega _m^L`$ and $`\omega _m^T`$, respectively, and according to dynamic finite-size scaling we have $`\omega _m^L=L^{z^L}\mathrm{\Omega }^L(qL)`$ and $`\omega _m^T=L^{z^T}\mathrm{\Omega }^T(qL)`$. For $`n_q=2`$ (see Fig.7) we obtain $`z^L=1.48\pm 0.14`$ and $`z^T=0.03\pm 0.25`$, where data for $`L=36,48`$, and $`60`$ have been included in the analysis. Our data show that the dynamic critical exponent for the transverse component is consistent with $`z=0`$, indicating that this component is not critical. In contrast, the longitudinal component is critical, with $`z1.5`$.

Separating $`S(𝐪,\omega )`$ into longitudinal and transverse components with respect to the uniform magnetization we see that the longitudinal component has a central peak, a pronounced spin-wave peak and less intense two-spin-wave peaks. In contrast, each such peak in the transverse component is split into two peaks, with frequencies $`\omega ^L\pm \mathrm{\Delta }\omega `$, where $`\omega ^L`$ is the frequency of the corresponding peak in the longitudinal component of $`S(𝐪,\omega )`$. The shift $`\mathrm{\Delta }\omega `$ in the peak frequencies corresponds to the frequency of oscillation of the staggered magnetization.

## IV Conclusions

We have used Monte Carlo and spin dynamics simulations to study the dynamic behavior of the isotropic Heisenberg antiferromagnet on the simple cubic lattice. When we use the same range of momentum $`q`$ as probed by experimentcoldea , the dynamic critical exponent obtainedsd1  is in good agreement with its experimental value, which is slightly lower than theoretical predictions. In our present work we have used a larger lattice size, and thus smaller values of $`q`$, in addition to obtaining better statistics and extrapolating finite $`q`$ results to the limit $`q=0`$. This allowed us to study systematic changes as we approach the asymptotic critical region and our improved estimate (i.e. with systematic errors largely eliminated) is $`z=1.49\pm 0.03`$, in good agreement with the renormalization group theory and dynamic scaling predictions. Presumably the values of $`q`$ used in the experimentcoldea  and in our previous simulations were not in the true asymptotic critical region, resulting in a slightly lower estimate of $`z`$. This should serve as a warning for future simulational and experimental probes of dynamic critical behavior.

Longitudinal and transverse components of $`S(𝐪,\omega )`$ with respect to the staggered magnetization are investigated separately. This required rotation of the frame of reference after each integration step because the staggered magnetization is not a conserved quantity. Below $`T_c`$, the transverse component $`S^T(𝐪,\omega )`$ has a pronounced spin-wave peak, whereas the longitudinal component $`S^L(𝐪,\omega )`$ is dominated by two-spin-wave addition and subtraction peaks, and a central peak presumably due to spin diffusion. These results are consistent with theorymazenkoTb  and experimentcox , both of which have shown that the transverse spin fluctuations are propagating, dominated by spin-waves, whereas the quasielastic peak is due to longitudinal fluctuations. At $`T_c`$, $`S^L(𝐪,\omega )`$ has a central peak that is much more intense than the spin-wave peak in $`S^T(𝐪,\omega )`$ and no central peak was seen in $`S^T(𝐪,\omega )`$. Explaining the appearance of a central peak in $`S^L(𝐪,\omega )`$ at $`T_c`$ remains a challenge for theory. We have also seen that while $`S^T(𝐪,\omega )`$ is not critical, $`S^L(𝐪,\omega )`$ is critical and it has a dynamic critical exponent $`z1.5`$.

These findings further support our earlier conclusion that a simple, nearest-neighbor, isotropic Heisenberg model describes the behavior of RbMnF3 quite well. The only limitations in the agreement appear to be at $`T_c`$, and even there all qualitative features and dynamic exponent are faithfully reproduced.

Below $`T_c`$ our results for the longitudinal and transverse components of the integrated intensities of $`S(𝐪,\omega )`$ are consistent with renormalization-group theory predictions, and not with the mean-field one. In contrast, at $`T_c`$, while the integrated intensities are consistent with the renormalization-group theory prediction, the large error bars do not allow us to exclude the mean-field behavior.

## V Acknowledgments

Fruitful discussions with A. Cuccoli are gratefully acknowledged. This research was partially supported by NSF grant DMR-0094422. Simulations were performed on the Cray T90 at SDSC and on the IBM SP at the U. of Michigan.