text:
# Surface extrapolation length and director structures in confined nematics

## I Introduction

The properties of nematic liquid crystals in confined geometries continues to be an interesting problem for both basic science and technological reasons . Nematics confined to cylindrical cavities with homeotropic radial boundary conditions can exhibit a variety of nontrivial structures depending on the competition between bulk elastic and surface energies. By minimizing the Frank elastic free energy Cladis and Kleman  and Meyer  showed that for a cylinder with a sufficiently large radius an escaped radial (ER) configuration will form. This configuration can be thought of as a planar radial (PR) structure (i.e. a disclination line) which has “escaped” along the axis of the cylinder (see Fig. 1). Later studies  showed the possibility of additional structures including planar polar (PP) and planar polar with two line defects (PPLD). In these latter configurations the director lies in a plane perpendicular to the axis of the cylinder with a strong component along a single in-plane direction (see Fig. 2). In the PP structure the local director is uniform near the axis of the cylinder and radial at the boundary. In the PPLD structure there are two half-integer disclination lines parallel to the cylinder axis.

Kralj and Zumer  carried out a very complete numerical stability analysis of the various nematic structures in a cylindrical geometry by minimizing the Frank elastic energy (including the saddle-splay elastic constant $`K_{24}`$). For most values of the elastic constants they found a phase diagram with ER and PP structures, with the ER structure stable for large radii and/or strong anchoring. In particular if the bend and splay elastic constants are equal and $`K_{24}=0`$, then the ER structure should form when $`RW/K`$ exceeds about 27 , where $`R`$, $`W`$ and $`K`$ are the cylinder radius (measured in units of intermolecular spacing), surface anchoring strength, and Frank elastic constant respectively. However, a PR structure appears if the twist constant is very large compared with splay. They also found that the PPLD structure can be stabilized if $`K_{24}`$ is nearly zero, R is approximately 100, and the half-integer defect lines are separated by a distance approximately equal to the diameter of the cylinder. An analytic study of these director structures was carried out by Burylov  with similar results.

Monte Carlo (MC) simulations of a cylindrically confined nematic with radial homeotropic boundary conditions were first performed by Chiccoli et al.  using the Lebwohl–Lasher lattice model . This model represents each mesogenic molecule or small group of molecules by a three-dimensional spin vector. These spins are free to rotate about their centers which are fixed on the sites of a lattice. Spins whose distance $`r`$ from the center of the cylinder is less than $`R`$, the radius of the cylinder, interact with their nearest-neighbors through the usual bulk Lebwohl-Lasher pair potential:

$$U_{ij}=ϵ_bP_2(\mathrm{cos}\theta _{ij})$$
(1)
where $`\theta _{ij}`$ is the angle between two spins $`i`$ and $`j`$, $`P_2`$ is a second rank Legendre polynomial and $`ϵ_b`$ is a positive constant for nearest–neighbor sites and zero otherwise. The homeotropic surface anchoring is introduced via a group of boundary layer spins located at distances $`R<r<R_c`$ which point in a fixed radial direction. The interaction of these spins with the spins inside the cylinder is described by the potential:

$$U_{ik}=ϵ_sP_2(\mathrm{cos}\theta _{ik})$$
(2)
where spin $`i`$ is located inside the cylinder and spin $`k`$ is the nearest–neighbor spin which belongs to the boundary layer. The coupling $`ϵ_s`$ is the surface anchoring strength. The total potential energy is given by the sum of the two energies above, summed over all pairs of spins. Note that in the Lebwohl-Lasher model, because the energy is invariant under a uniform rotation of all the spins, the bend, splay and twist elastic constants are equal, and the saddle splay elastic constant $`K_{24}`$ is identically zero.

Chiccoli et al. studied cylinders of radii 6 and 11 lattice spacings and heights 40, 52, 62, and 82 lattice spacings with periodic boundary conditions along the cylinder axis. The spins were placed on the sites of a simple cubic lattice and the director structures were studied at a temperature of 0.6, measured in units of $`ϵ_b`$. Various values of the ratio of the surface anchoring strength $`ϵ_s`$ to the bulk nearest-neighbor coupling $`ϵ_b`$ were considered, the largest being unity. In all cases, a planar structure was observed; however, these authors did not determine whether this structure is PP or PPLD.

Subsequently Smondyrev and Pelcovits  carried out Monte Carlo simulations on much larger Lebwohl–Lasher systems up to radii $`R=160`$, with a height of 16, again using a simple cubic lattice. For low values of the radius or anchoring strength they found a PPLD structure, whereas for $`RW/K=160,ϵ_s/ϵ_b=1`$, and a temperature of 0.9, they obtained a metastable ER configuration. But this metastable structure eventually collapsed to a planar state, apparently in contradiction with elastic theory which  predicts that the ER structure should be stable when $`RW/K`$ exceeds approximately 27. If we assume (incorrectly, as we show in Sec. II) that the ratio $`ϵ_s/ϵ_b`$ equals the corresponding Frank elastic energy ratio $`W/K`$, then Smondyrev and Pelcovits were simulating a system with $`RW/K=160`$, far in excess of the critical value of 27.

Bradač et al.  carried out molecular dynamics simulations on cylindrically confined nematics interacting with a modified induced–dipole–induced–dipole interaction. For systems quenched from the isotropic phase to temperatures deep in the nematic phase, a phase diagram was obtained that is in reasonable agreement with the predictions of elastic theory. For the largest systems studied (radii of 16 lattice spacings) the PP, PPLD and ER structures were observed with increasing values of the surface interaction. For the case of equal elastic constants (i.e. like the Lebwohl–Lasher model) the ER structure appears when the ratio of the surface to bulk interaction energies is 0.2. However, these simulations were carried out for only 10000 time steps leaving open the possibility that the ER structure ultimately collapses. Also, these simulations were carried out at apparently lower temperatures than those used in the Monte Carlo simulations.

In this paper we report the results of MC simulations of the confined Lebwohl–Lasher model which address the question of why the ER structure was not observed in previous MC simulations. We find that the ratio of the bulk–to–surface couplings in the Lebwohl–Lasher model, $`ϵ_b/ϵ_s`$, is not equal to $`K/W`$, except at very low temperatures. The latter quantity is equivalent  to the surface extrapolation length $`b`$. We find that there is significant dependence of $`b`$ on temperature in the Lebwohl–Lasher model; in particular $`b`$ grows as the NI transition temperature is approached from below (a feature observed in experimental studies of confined nematics ). Thus, identifying the parameter $`Rϵ_s/ϵ_b`$ with $`RW/KR/b`$ is in general incorrect. Except at very low temperatures $`RW/K`$ will be less than $`Rϵ_s/ϵ_b`$ (substantially so near the NI transition). Comparing the results of simulations with the predictions of elastic theory then requires consideration of this difference.

We also find, for large enough systems, that the rate at which the system is cooled from the isotropic to the nematic phase determines which director structure is ultimately observed in a simulation. Because $`b`$ is large near the NI transition, a slow cooling rate will first yield a PPLD structure in this temperature range, in agreement with elastic theory which predicts that the PPLD is stable for small values of $`R/b`$. However, because of the large free energy barrier between the ER and PPLD structures, continued slow cooling to lower temperatures (where $`b`$ is small and elastic theory predicts the appearance of the ER structure for sufficiently large $`R`$) will not yield the expected ER structure; the system becomes trapped in the metastable PPLD state. On the other hand, a relatively rapid cooling from the isotropic phase to low temperatures bypasses the PPLD state and allows the ER structure to form. To our knowledge, this is the first time a stable ER structure has been seen in a MC simulation of confined nematics.

Our paper is organized as follows. In the Sec. II we describe our measurements of the surface extrapolation length and its temperature dependence. In Sec. III we report on simulations of a cylindrically confined nematic, exploring the effect of different cooling rates on the formation of the ER and PPLD structures. We offer some concluding remarks in the final section.

## II Surface Extrapolation length

To measure the temperature dependence of the surface extrapolation length $`b`$, we consider a modification of the discussion presented in Ref. , where this length is defined by examining the director profile in a nematic confined between two plates which impose a twist distortion that is uniform across the sample except near the plates. The twist geometry is appropriate for studying planar anchoring at a surface. As we are interested in homeotropic alignment we consider instead a splay geometry. As shown in Ref.  the extrapolation length can be defined by extrapolating the constant slope portion of the director profile corresponding to the uniform twist or splay distortion; see Fig. 3. By minimizing the total elastic energy consisting of the bulk Frank energy in the one–elastic constant approximation and the surface energy, it was shown in Ref.  that the length $`b`$ so defined is related to the Frank elastic constant $`K`$ and the surface energy $`W`$ by $`b=K/W`$.

We measured $`b`$ defined in the above fashion in the Lebwohl–Lasher model on a cubic lattice of linear size 40. On the plane $`z=40`$ we introduced the homeotropic surface interaction, Eq. (2), with boundary spins aligned along the $`z`$ axis. Spins on the opposite wall, $`z=0`$, were rigidly aligned parallel to the $`x`$ axis. These spins interact with their interior neighbors via the bulk interaction, Eq. (1). Periodic boundary conditions were imposed on the remaining four sides of the lattice. Starting in the isotropic phase the system was cooled down in temperature steps of 0.05. We equilibrated at least $`200000`$ MC steps at each value of the temperature. We display results in Fig. 3 for three values of temperature ranging from deep in the nematic phase to near the NI transition. For each value of $`z`$ we plot the average value of the spin’s angle with the $`z`$ axis (averaging over all sites in the $`xy`$ plane). We note that $`b`$ appears to diverge near the NI transition temperature in agreement with experimental results . Additional data for $`b`$ as a function of $`T`$ and $`ϵ_s`$ is shown in Table I. As expected, $`b`$ increases with decreasing surface coupling.

The strong dependence of $`b`$ on temperature indicates that $`W`$, the surface energy, has a different dependence on temperature than the Frank elastic constant $`K`$. The latter was measured in the Lebwohl–Lasher model by Cleaver and Allen , who found it to be approximately proportional to $`S^2`$, where $`S`$ is the nematic order parameter. Experimental measurements of $`b`$  find that it diverges as $`(T_{NI}T)^1`$, where $`T_{NI}`$ is the NI transition temperature, which suggests that $`W`$ is proportional to $`S^4`$. Our data is not sufficiently robust to determine the dependence of $`W`$ on $`S`$.

The temperature dependence of $`b`$ helps to explain the apparent discrepancy between the results of previous simulations  and elastic theory , summarized in Sec. I. Even at a relatively low temperature $`T=0.6`$ as was used in Ref. , $`b`$ is of order 2, and thus the elastic theory parameter $`RW/K=R/b`$ that determines the relative stability of the various director structures is approximately half the value of the corresponding Lebwohl–Lasher “bare” parameter, $`Rϵ_s/ϵ_b`$. Ref.  reported on simulations carried out at $`T=0.9`$ with $`ϵ_s/ϵ_b=1.0`$, where the Lebwohl–Lasher parameter is approximately 6 times greater than $`RW/K`$. Thus the simulations carried out in this latter work at a value of $`R=160`$, in reality correspond to $`RW/K27`$, i.e., approximately the threshold value for the formation of the ER structure. On the other hand, the simulations of Ref.  were carried out deep in the nematic phase where $`RW/K`$ and $`Rϵ_s/ϵ_b`$ are nearly equal.

## III Stability of the ER and PPLD structures

With the insight gained in the previous section regarding the proper comparison of simulation results with elastic theory, we carried out additional MC simulations on the Lebwohl–Lasher model confined to a cylinder. Previous simulation studies of nematics confined to cylinders  were done on a cubic lattice. However, this choice of lattice biases the position of the two disclination lines in the PPLD structure; in particular they always appear on one of the two crystalline axes perpendicular to the long axis of the cylinder. In addition approximately half of the spins on the perimeter of the cylinder have two nearest–neighbor spins on the fixed boundary layer (see Eq. (2)), whereas the other half have only one such nearest–neighbor, resulting in a nonuniform coupling of the nematic to the cylindrical surface.

To eliminate this nonuniformity we worked instead with a honeycomb lattice. With this choice nearly all spins on the perimeter of the cylinder have only one nearest–neighbor in the boundary layer. While the honeycomb lattice does not have continuous rotational symmetry about the cylinder axis, we found that its higher symmetry about this axis (six–fold as compared to four–fold for the cubic lattice) led to the appearance of the disclination lines in the PPLD structure on axes not necessarily along the crystal lattice directions. The honeycomb lattice is bipartite allowing us to implement a standard “checkerboard” algorithm , an efficient choice for vector processors. We divided the system into two interconnected sublattices and then alternately updated spins on each.

As in Ref.  we determined the nature of the resulting director structures by computing an escaped order parameter, which measures the degree of spin tilt away from the plane perpendicular to the cylinder axis, and is a sensitive measure of the presence of the ER structure. It is given by:

$$P_2^{esc}=(1/N)P_2(\widehat{𝐮}_𝐢\widehat{𝐳})$$
(3)
where the cylinder axis is along the $`z`$ direction. The escaped order parameter is calculated as a function of the distance from the center of the cylinder, and the sum is over spins in a thin annulus with the same height as the cylinder.

We explored the stability of the PPLD and ER structures for systems of two different radii, 80 and 120, both of height 16, with $`ϵ_s/ϵ_b=0.9`$. In Table II we tabulate the values of $`R/b`$ for these two radii using the data for $`b`$ from Table I. Using the table we see that elastic theory predicts that the PPLD structure should be stable for temperatures greater than 0.6 and 0.8, for $`R=80`$ and 120, respectively, while the ER structure should be stable at temperatures below these values. Note that elastic theory by its nature is unable to predict the transition temperature separating the PPLD and ER structures from the isotropic phase. Beginning in the isotropic phase at $`T=1.4`$ we cooled the system of $`R=80`$, in temperature steps of 0.1 equilibrating $`100000`$ MC steps at each temperature. At $`T=0.9`$ a stable PPLD structure was formed. We checked the stability of this structure up to 500000 MC steps. Upon further cooling down to $`T=0.5`$ (again equilibrating 100000 MC steps at each temperature) the system remained trapped in the PPLD structure even though at this lower temperature $`b=1.1`$ and $`R/b=72.7`$, well above the threshold predicted by elastic theory for the formation of the ER structure. We found similar behavior irrespective of the cooling rate used, including instantaneous quenches directly from the isotropic phase. Conversely, if we began at $`T=0.8`$ with an ER configuration constructed by hand using the mathematical form derived in Refs. , we found this state to be stable even after 500000 MC steps. Yet for this value of $`T`$ the PPLD state should be stable ($`R/b=20.0`$). These results suggest that there is a large free energy barrier separating the PPLD and ER structures, and the system can be trapped in either state.

We found somewhat different results for the larger radius system, $`R=120`$. We cooled this system down from $`T=1.4`$ again in temperature steps of 0.1 equilibrating 300000 MC steps at each value of temperature. As expected a stable (up to 500000 MC steps) PPLD structure was formed at $`T=0.9`$ ($`R/b20<27`$). Upon further cooling the system remained trapped in the PPLD state, just as in the case of the smaller system. However, if the same system were equilibrated $`100000`$ MC steps at each temperature step as it was cooled from $`T=1.4`$ to $`T=0.8`$, then an ER structure was formed at the latter temperature and remained stable up to 500000 MC steps. In this case the more rapid rate of cooling allowed us to bypass the incipient PPLD structure at $`T=0.9`$ and reach the ER structure which should be, according to elastic theory, the globally stable one at $`T=0.8`$ (where $`b4`$, and $`R/b30`$). We suspect that our ability to form an ER structure for this system of larger radius is due to the smaller temperature range ($`0.85\pm 0.05`$) where the PPLD structure is stable compared with the case of $`R=80`$, where the PPLD structure, according to elastic theory and Table II, should be stable from $`T0.9`$ down to $`T0.7`$.

## IV Conclusions

In this paper we have addressed the question of reconciling the predictions of elastic theory with the results of numerical simulations for nematics confined to cylindrical cavities with homeotropic boundary conditions. Our main result is that care must be taken in comparing elastic theory with numerical simulations because of the strong temperature dependence of the ratio $`K/W`$, the surface extrapolation length. Except at very low temperatures, this ratio is much larger than the corresponding “bare” ratio of bulk and surface energies in a numerically simulated model. When this difference is taken into account, the results of numerical simulations are consistent with the predictions of elastic theory. In particular, the ER structure forms for values of $`RW/K`$ above the critical value of 27 predicted by elastic theory, providing that the metastable PPLD structure can be avoided. Because $`K/W`$ is very large near the NI transition, and thus $`RW/K`$ is relatively small, the system will readily form a PPLD structure if it is cooled slowly from the isotropic phase. However, with rapid cooling and for large enough systems where the PPLD phase is limited in its temperature range, this state can be bypassed and a stable ER structure formed at lower temperatures where $`RW/K`$ is sufficiently large.

## Acknowledgments

We thank Professor G. P. Crawford for helpful discussions and Dr. G. B. Loriot for computational assistance. This work was supported by the National Science Foundation under grant No. DMR98-73849. Computational work in support of this research was performed at Brown University’s Theoretical Physics Computing Facility.


format: text

text:
# New algorithm and results for the three-dimensional random field Ising Model

## I Introduction

Despite twenty five years of experimental and theoretical effort, phase transitions in systems with quenched random fields are still poorly understood. The simplest theoretical model is the random field Ising model (RFIM). The RFIM phase transition is believed to be in the same universality class as the phase transitions in diluted antiferromagnets in a uniform field and fluids in porous media. The three dimensional RFIM is known  to have an ordered phase at sufficiently low temperature and for weak random fields. As the temperature or the strength of the randomness is increased, there is a transition to a disordered phase. The nature of this transition is not well understood.

In this paper we describe a new replica exchange algorithm for simulating the RFIM and present numerical results for systems of $`L\times L\times L`$ spins with $`L`$ up to 24. We show that the qualitative features of the transition differ strongly from realization to realization of disorder. Our results can be interpreted as suggesting that the RFIM transition is first-order or that a modified version of the droplet picture holds.

The RFIM is described by the energy

$$/k_BT=\beta \underset{ij}{}S_iS_j+\mathrm{\Delta }\underset{i}{}h_iS_i+H\underset{i}{}S_i$$
(1)
where the spin variables $`S_i=\pm 1`$ reside on a lattice, $`\beta =1/k_BT`$ is the inverse temperature, the first sum is over nearest neighbor pairs on the lattice, $`H`$ is an external field and $`\mathrm{\Delta }`$ is the strength of the disorder. The random fields are independent random variables chosen from a distribution with mean zero and variance one. In this study, the random fields are Gaussian and the lattice is simple cubic.

Currently, it is not known whether the phase transition for the 3D RFIM is first-order or continuous. Experiments on magnetic systems have been plagued by problems of poor equilibration and have yielded confusing results, but it appears that there is no latent heat at the transition and this has usually been interpreted as evidence for a continuous transition. Theoretical analyses  have also favored a continuous transition, although in many cases this is an assumption rather than a conclusion, and some recent work  suggests a fluctuation driven first-order transition. The standard picture is that of a continuous transition controlled by a zero temperature fixed point. The scaling theory of this transition  has three independent critical exponents and modified hyperscaling relations. Because of the fixed point governing the transition has strong disorder, controlled renormalization group calculations have not been possible. Migdal-Kadanoff renormalization group calculations indicate a continuous transition  but also mistakenly predict that the $`q`$-state Potts transition is continuous for all $`q`$ when, for 3D it is known to be first-order for $`q3`$. Series analyses initially supported a first-order transition  but more recently point to a continuous transition, at least for weak disorder . Alternatively, it may be that the transition is continuous for weak disorder and then becomes first-order for strong disorder with a tricritical point separating the critical line from the first-order line.

Recent Monte Carlo simulations  have also been interpreted as showing a continuous transitions but with a jump in the magnetization. These simulations have been limited to system size $`16^3`$. The jump in the magnetization can be interpreted as a very small value of the magnetization exponent but might also signal a first-order transition. Simulations of systems up to $`64^3`$  were interpreted as indicating a first-order transition but were clearly not equilibrated in the transition region. A number of numerical experiments have also been carried out on the RFIM at zero temperature using polynomial time ground state algorithms. Originally these supported a continuous transition with a small magnetization exponent  but more recent studies on larger systems  show a jump in the magnetization and are thus suggestive of a first-order transition.

The numerical results presented in this paper differ from those of previous studies in two significant respects. First we use an efficient algorithm that permits us to reach equilibrium for larger systems than those of past studies. Second, we fine-tune both the temperature and the external field for each realization of disorder to be as close as possible to the transition for that realization.

To understand the motivation for this fine-tuning, let us suppose for the moment that the transition is first-order. If this is the case, then for periodic or helical boundary conditions we expect that there will be three phases in coexistence at the transition point. We call the coexisting phases $`+`$, $``$ and $`0`$. The ordered phases, $`+`$ and $``$, have long range order and finite magnetization. The disordered, $`0`$ phase has no long range order, no magnetization and is characterized by spins that are predominately aligned with their local fields. The bond energy, defined by the first term on the RHS of Eq. (1), is greater in the $`0`$ phase than the ordered phases. The expected phase diagram in the $`T`$$`H`$ plane in the vicinity of the point of three phase coexistence is shown in Fig. 1. For $`T<T_c`$ and $`H=0`$ there is phase coexistence of the $`+`$ and $``$ phases that ends at the thermal first-order transition at $`T=T_c`$ and $`H=0`$ (the black dot in the figure). Since the $`0`$ and $`+`$($``$) phases differ in both energy and magnetization, the disordered phase can be maintained in coexistence with the $`+`$($``$) phase by increasing the temperature and increasing (decreasing) the field. The two order–disorder lines corresponding to $`0+`$ and $`0`$ coexistence form the arms of the “Y” in the figure.

Since the disorder is independent and homogenous the free energy is self–averaging, hence, in the thermodynamic limit the transition occurs at a definite point ($`T_c`$, $`H=0`$) for almost all realizations of the random fields. However, for finite systems, the location and qualitative features of the transition will depend on the realization of disorder. In particular, for a system of size $`L`$, with realization of random field, $`\{h_i\}`$ and disorder strength $`\mathrm{\Delta }`$ there may be a point where three phases coexist at ($`T_c(\mathrm{\Delta };\{h_i\})`$, $`H_c(\mathrm{\Delta };\{h_i\})`$). A first guess, based on the net field due to the random fields, is that $`T_c(\mathrm{\Delta };\{h_i\})`$ and $`H_c(\mathrm{\Delta };\{h_i\})`$) should be displaced from the average value $`\stackrel{~}{T}_c(L)`$ and $`H=0`$ by an amount of order $`\mathrm{\Delta }L^{d/2}`$.

How accurately must $`H`$ and $`T`$ be fine-tuned to see three phase coexistence if it exists? Suppose first that we want the $`+`$ and $``$ phases to coexist and that the magnetizations of these phases at the transition differs by $`2m`$. If the external field deviates from the correct value by $`\delta H`$, the free energy difference between the $`+`$ and $``$ phase is $`2m\delta HL^d`$. For both phases to have a significant probability the free energy difference must not greatly exceed $`kT`$. Thus $`H`$ must be set to $`H_c(\mathrm{\Delta };\{h_i\})`$ to within an accuracy of $`kT/2mL^d`$ to have both phases represented. $`H_c(\mathrm{\Delta };\{h_i\})`$ will itself fluctuate from sample to sample as $`\mathrm{\Delta }/L^{d/2}`$. Similarly, if the entropy difference per spin between the the ordered and disordered phases is $`s`$, then $`T`$ must be fine-tuned to within $`kT/sL^d`$ to allow to these phases to coexist. Presumably sample to sample fluctuations in $`T_c(\mathrm{\Delta };\{h_i\})`$ also scale as $`1/L^{d/2}`$. Thus if a single value of $`H=0`$ and $`T=\stackrel{~}{T}_c(L)`$ is chosen for all realizations of the random field, one will almost never see more than one phase at a time.

## II Numerical methods

We use an algorithm that combines the replica exchange method, first introduced by Swendsen and Wang  and the two-replica cluster method of Redner, Machta and Chayes . Our method is also closely related to simulated and parallel tempering . The idea of this approach is to simultaneously simulate $`K`$ replicas of the system. All replicas have the same normalized random field $`\{h_i\}`$ but each replica has different values of the other parameters. Replica $`k`$ ($`k=0,\mathrm{},K1`$) has inverse temperature $`\beta _k`$, strength of randomness $`\mathrm{\Delta }_k`$ and external field $`H_k`$. The replicas form a sequence so that neighboring replicas in the sequence are nearby in the $`(\beta ,\mathrm{\Delta },H)`$ parameter space. Neighboring replicas exchange magnetization with one another according to a procedure described below. One end of the sequence of replicas is at a value of $`\beta `$, $`\mathrm{\Delta }`$ and $`H`$ that can be efficiently simulated using a known method while the other end of the sequence is at a value of the parameters that we would like to study. In our case, the replicas lie along the RFIM phase transition line starting from the pure Ising values $`\beta _0.22615`$ and $`\mathrm{\Delta }_0=H_0=0`$ as shown in Fig. 2. The replicas are equally spaced in $`\mathrm{\Delta }`$. The pure Ising replica is simulated using the Wolff single cluster algorithm . (We have also experimented with replicas lying along a line of constant $`\beta H`$ starting at the RFIM phase boundary, extending into the paramagnetic phase, and ending at a temperature high enough that the model can be efficiently simulated using the ordinary single-spin-flip Metropolis algorithm. This approach, however, is found to be less efficient than the one described above.)

Magnetization is exchanged between neighboring replicas using a generalization of the two-replica cluster method. Suppose we have two replicas at $`(\beta ,\mathrm{\Delta },H)`$ and $`(\beta ^{},\mathrm{\Delta }^{},H^{})`$ with spin configurations $`\{S_i\}`$ and $`\{S_i^{}\}`$, respectively. A site $`j`$ is considered active for this pair of replicas if $`S_jS_j^{}`$. A bond $`ij`$ between neighboring sites $`i`$ and $`j`$ is satisfied if $`S_i=S_j`$ and $`S_i^{}=S_j^{}`$. A cluster of active sites is formed starting from a randomly chosen active site. New active sites are added to the cluster by occupying satisfied bonds on the perimeter of the cluster with probability $`p(\beta ,\beta ^{})`$ where

$$p(\beta ,\beta ^{})=1e^{2(\beta +\beta ^{})}.$$
(2)
If a bond connecting a site to a cluster is occupied, the site is added to the cluster and the set of bonds on the perimeter is updated. In this way, the cluster grows until no further sites are added. The procedure is very similar to the way clusters are grown in the Wolff single cluster algorithm.

Once a cluster is identified it is flipped with a probability that depends on the change in boundary and field energy so as to satisfy detailed balance. Flipping a cluster means changing the sign of all the spins in the cluster or, equivalently, exchanging the values of the spins in the cluster between the two replicas. The probability to flip a cluster $`𝒞`$ with $`𝒞`$ sites depends on the quantity $`\mathrm{\Sigma }`$

$$\mathrm{\Sigma }=[(\mathrm{\Delta }\mathrm{\Delta }^{})h_𝒞+(HH^{})𝒞+(\beta \beta ^{})(N_{++}N_{})]S_𝒞$$
(3)
where $`N_{++}`$ and $`N_{}`$ are, respectively, the number of $`++`$ and $``$ sites that are nearest neighbors of the cluster, $`h_𝒞`$ is the net random field acting on the cluster,

$$h_𝒞=\underset{i𝒞}{}h_i,$$
(4)
and $`S_𝒞`$ is the spin value of the cluster in the unprimed replica. If $`\mathrm{\Sigma }0`$ then the cluster is flipped, otherwise it is flipped with probability $`e^\mathrm{\Sigma }`$.

It is straightforward but tedious to show that the choice of $`p(\beta ,\beta ^{})`$ (Eq. 2) together with the flipping probability defined by $`\mathrm{\Sigma }`$ is precisely what is needed to ensure detailed balance. The motivation for the choice of $`p(\beta ,\beta ^{})`$ is most easily understand by considering the limit where two replicas are at the same values of $`\beta `$, $`\mathrm{\Delta }`$ and $`H`$. In that case $`p(\beta ,\beta )=1e^{4\beta }`$ and $`\mathrm{\Sigma }=0`$. Since $`\mathrm{\Sigma }=0`$, clusters are always flipped just as is the case for the Wolff single cluster algorithm. Furthermore, we have shown in Ref.  that the active clusters percolate at the RFIM phase transition. If the transition is continuous clusters of all sizes are flipped. If the transition is first-order, there will be two distinct kinds of clusters; some clusters will be extensive and change the phase of the system while other clusters will have sizes less than or equal to the correlation length. In either case, the clusters identified by the two-replica procedure correspond to the fluctuations that actually occur in the system at the phase transition and permits large changes in the spin configuration in a single Monte Carlo sweep. When the two replicas do not have equal values of the parameters, then the clusters do not flip freely but if the replicas are close together in the parameter space, the acceptance fraction for flipping clusters will remain high.

Our method and the original replica exchange method  on which it is based is similar to parallel tempering . In all these methods, groups of spins are exchanged between neighboring replicas along a sequence. In parallel tempering the whole spin configuration is exchanged and the Boltzmann factor controlling the acceptance of the move depends on the energy difference between the replicas. In our algorithm, only some of the spins are exchanged and for a given distance in the parameter space between the replicas the acceptance fraction is larger than for parallel tempering. The consequence is that fewer, less closely spaced replicas are needed for the replica exchange method.

In order to find the phase transition temperature and external field ($`\beta _c(\mathrm{\Delta };\{h_i\})`$, $`H_c(\mathrm{\Delta };\{h_i\})`$) for a given realization and strength of disorder we use a feedback mechanism. Starting from an initial value of $`\beta `$ and $`H`$ we monitor the magnetization of the system after each Monte Carlo sweep. If the absolute value of the magnetization is less than a lower cut-off $`M_{lc}`$ the system is interpreted to be in a high temperature phase and the inverse temperature is increased by a small amount $`ϵ_\beta `$. If the absolute value of the magnetization is greater than an upper cut-off $`M_{uc}`$ the system is interpreted to be in one of the ordered phases and $`\beta `$ is decreased by $`ϵ_\beta `$. In this case, the external field is also adjusted by an amount $`ϵ_H`$ if the magnetization is positive and by $`ϵ_H`$ if the magnetization is negative. The average value of $`\beta `$ and $`H`$ is computed for the period when the feedback procedure is on and is taken to be ($`\beta _c(\mathrm{\Delta };\{h_i\})`$, $`H_c(\mathrm{\Delta };\{h_i\})`$). The feedback procedure is then turned off and these values are used for a long equilibrium simulation.

Most of our simulations were for $`24^3`$ systems. Except as otherwise noted, the simulations used $`K=16`$ replicas equally spaced in $`\mathrm{\Delta }`$ with the most disordered replica having $`\mathrm{\Delta }=.35`$. Initially the replicas lie on an elliptical curve in the temperature-disorder plane that starts at the pure Ising transition ($`\beta =0.22165`$) and ends at the zero temperature transition ($`\mathrm{\Delta }/\beta =2.35`$) as shown in Fig. 2. It was shown in Ref.  that this curve is a good estimate for the phase boundary. The feedback parameters were $`M_{lc}=2000`$, $`M_{uc}=4000`$, $`ϵ_\beta =10^4`$ and $`ϵ_H=10^5`$. The feedback mechanism was run for $`2\times 10^5`$ Monte Carlo sweeps and then, with the temperature and external field determined by the feedback procedure, data were collected for $`8\times 10^5`$ Monte Carlo sweeps. Each Monte Carlo sweep is defined as 500 cluster steps; since the average cluster size was found to be roughly 1000, a single Monte Carlo sweep corresponds to attempting to flip every spin in every replica about once. For each cluster step a random integer $`k`$ was chosen between 0 and $`K`$. If $`k=0`$, a Wolff move was performed on the pure replica, if $`1<k<K`$, a two replica cluster move was carried out between replica $`k`$ and replica $`k1`$. If $`k=K`$, a two replica cluster move was carried out between replicas $`K1`$ and $`K2`$. The procedure described above, totaling $`10^6`$ sweeps, was performed for 22 realizations of disorder, a number chosen in advance of the experiment, labeled 10 through 31 by the seed of the random number generator that produced $`\{h_i\}`$. Each realization required about six days on a 450MHz Pentium III machine. A number of additional simulations were carried out for specific realizations of disorder (not all in the range 10 through 31) using parameters that differed in some way from the above. In particular, we did a careful study of realization 1 that is described below.

## III Results

### A Energy and magnetization histograms

The main lesson of our work is that each realization of disorder has its own character. This is best seen by examining the probability distributions of magnetization and energy for several realizations of disorder. Figures 3, 4, 5, 6 and 7 show the magnetization histogram for realizations 20, 21, 25, 31 and 14, respectively. The magnetization histogram for realization 20, characterized by two broad maxima roughly symmetrical about the origin, is typical of the majority of realizations. Apart from the asymmetry in the two peaks, realization 20 and those like it are similar to the pure Ising model whose magnetization histogram is shown in Fig. 8 at the infinite system size critical temperature, $`\beta _c(0)=.22165`$. The magnetization histogram for realization 21, Fig. 4 also displays two peaks but now one peak is much sharper than the other and quite asymmetric about the origin. The magnetization histogram for realizations 25, Fig. 5 and 31, Fig. 6 are qualitatively different, displaying three distinct and well separated peaks. It should be noted that the feedback procedure was initially unsuccessful for these realizations and had to be run again with $`M_{lc}=M_{uc}=5000`$ to find all three peaks. Failure of the feedback scheme was evidenced by the presence of only one narrow peak in the magnetization histogram. In the case of realization 31, the feedback procedure with the revised parameters failed on the most disordered replica so Fig. 6 shows the second to last replica in the sequence at $`\mathrm{\Delta }=.3267`$. For both realization 25 and 31, the original run showed three peaked structures at weaker disorder giving a strong hint that the same qualitative feature would be found at stronger disorder. Finally, realization 14, Fig. 7 was unique, displaying a five peaked structure. Out of a total of 22 magnetization histograms, 19 show two peaks, 2 three peaks and 1 five peaks. It is possible but unlikely that some of the two peak realizations would display more peaks with different fine-tuning. For example, we used other alternate values of the fine-tuning parameters to search for an additional positive peak for realization 21 but none was found.

The bond energy distribution displayed two possible qualitative behaviors. Figures 9, 10, 11, 12 and 13 show the bond energy histogram for realizations 20, 21, 25, 31 and 14, respectively. For comparison, Fig. 14 shows the energy distribution for the pure Ising model at the infinite system size critical temperature. For realization 20 and those like it, the bond energy distribution has a single peak and is qualitatively like that of the pure Ising model at criticality, except that the peak is shifted to significantly lower energies. For realization 21, which has two well separated magnetization peaks, there are also two well separated peaks in the bond energy distribution. Similarly, for realizations 25 and 31 there are two peaks in the bond energy distribution. In a first-order interpretation of these realizations the $`+`$ and $``$ phases have the lower energy and the $`0`$ phase the higher. On the other hand, for realization 14, which has five magnetization peaks, there is a single broad maximum in the bond energy histogram with a significant, low energy shoulder. Altogether, 5 out of 22 realizations show two well defined peaks in the bond energy histogram while the remaining realizations show a single peak though sometimes with a significant low energy shoulder.

The joint magnetization/bond-energy distribution for realization 25 is shown in Fig. 15. For comparison, the joint distribution for the pure Ising model at the infinite system critical temperature is shown in Fig. 16. The three lobes in the distribution correspond to the three peaks in the magnetization histogram, Fig. 5. It is tempting to interpret realizations 25 and 31 in the language of first-order transitions and declare the low energy side lobes as ordered phases and the central lobe as a disordered phase.

The importance of fine-tuning the magnetization and temperature for each realization of disorder is illustrated in Figs. 17, 18, 19 and 20. The magnetization histogram for realization 1 at the fine-tuned parameter values, $`\beta _c=0.268385`$ and $`H_c=0.00127049`$ is shown in Fig. 17. Note, that this value of $`H_c`$ is not the same as the value that would be obtained by summing the local fields. For realization one, that value is almost $`50\%`$ larger, $`\mathrm{\Delta }_ih_i/N=0.00178774`$. Figure 18 shows realization 1 with same external field but with the temperature raised by $`5\%`$ leaving only the disordered phases. Figure 19 shows realization 1 with the same external field but the temperature decreased by $`5\%`$. Although, the external field is too negative here, it is nonetheless clear that the two ordered phases are now dominant over the disordered phase, which has nearly vanished. Finally, if the external field is set to 0 and the feedback temperature is used, only the negative peak survives, as shown in Fig. 20.

We considered several realizations at stronger disorder. Figure 21 shows the magnetization histogram for realization 1 at disorder $`\mathrm{\Delta }=0.433`$ at the value of external field and temperature determined from the fine-tuning procedure. Comparing Fig. 21 to Fig. 17 we see that the three-peaked structure is conserved and the peaks becomes sharper as the strength of the disorder is increased. In this experiment, we used 16 replicas with maximum disorder strength, $`\mathrm{\Delta }=0.5`$ however, the fine-tuning procedure failed for the two most disordered replicas and only yielded the minus phase. In general, the feedback procedure was less stable for stronger disorder. Figure 22 shows realization 14 at disorder $`\mathrm{\Delta }=0.5`$ at the values of the temperature and external field determined from the fine-tuning procedure. A comparison to Fig. 7 shows that increasing the disorder has sharpened the five peaked structure of the magnetization histogram. It would be interesting to carry out a systematic study of stronger disorder to see if there is a trend toward more “first-order” like realization however this will require developing a better fine-tuning mechanism.

Overall, out of 21 realizations (seeds 10-30) at disorder strength $`\mathrm{\Delta }=.35`$ we found that the average values of the critical parameters is $`\beta _c=0.2663`$ and $`H_c=0.0006`$ (compared with initial values before the feedback procedure of $`\beta =0.2670`$ and $`H=0`$) with a standard deviations of $`\sigma _\beta =0.0037`$ and $`\sigma _H=0.0031`$. The standard deviation of $`H_c`$ is consistent with the anticipated value $`\sigma _H\mathrm{\Delta }L^{d/2}=0.0030`$.

### B Dynamics of the algorithm

We have studied the dynamics of the replica exchange algorithm and compared it to the Metropolis algorithm. There are two ways of comparing the algorithms. First, time can be measured by Monte Carlo sweeps. This approach favors the replica exchange algorithm since a single Monte Carlo sweep involves flipping spins in many (here 16) replicas and since growing clusters is computationally more intensive than flipping single spins. The second approach is to compare actual running time for the two algorithms on the same computer. This approach ignores the possibility that the one algorithm is better optimized and it is not an appropriate way to measure fundamental quantities such as the dynamic exponent however it does give a practical comparison and is useful for deciding which algorithm to choose for a given system size. For $`24^3`$ systems with 16 replicas with maximum disorder $`\mathrm{\Delta }=.35`$ on a 450MHz Pentium III machine we find that the replica exchange algorithm runs at about $`0.56`$ sec/sweep and that our implementation of the Metropolis algorithm runs at about $`0.0088`$ sec/sweep. As is the case for the equilibrium properties of the system, we find great differences in the dynamics depending on the realization of disorder. These differences are much more pronounced for the Metropolis algorithm. We measured the integrated autocorrelation time for the magnetization for realization 20 for both algorithms with the result $`\tau ^{\mathrm{Metropolis}}=6000`$ and $`\tau ^{\mathrm{replica}}=1100`$ when measured in Monte Carlo sweeps. The advantage of the replica exchange algorithm is lost however when these times are converted to running times; $`\tau ^{\mathrm{Metropolis}}=53`$s and $`\tau ^{\mathrm{replica}}=616`$s. The joint energy/magnetization histogram for realization 20 is smooth and without gaps suggesting that there are no global energy barriers separating regions of phase space. At the other extreme is realization 31 for which the joint energy/magnetization distribution is split into three pieces separated by gaps where the probability is very small. Using the parameters for realization 31 determined by fine-tuning for $`\mathrm{\Delta }=.3267`$, we did a long Metropolis run of 100 million Monte Carlo sweeps starting from random initial spins. The simulation stayed mainly in the “$`0`$” state with two excursions to the “$`+`$” state. The “$``$” state was never visited. On the other hand, two Metropolis runs of 20 million sweeps starting from all spins up and all spins down went to the “$`0`$” and “$``$” states, respectively, and stayed there for the entire simulation. Figure 23 shows the magnetization time series from the replica exchange algorithm for realization 31 for 800,000 sweeps. Although it is difficult to estimate the autocorrelation time from this series, it is clear that the simulation samples each of the states many times. The conclusion is that for most realizations of disorder, at size $`24^3`$ there is no great practical advantage to using the replica exchange algorithm but that for cases where well separated regions are present in the joint distribution, only the replica exchange algorithm is able to reach equilibrium within reasonable simulation times.

## IV Discussion

One of the chief motivations for this study was to determine the nature of the RFIM phase transition. Unfortunately, our results fail to settle this question. What would we have expected from the different scenarios for the phase transition? If the transition is first-order we should find three phases in coexistence at the transition point. The high temperature phase, would have little magnetization and a large bond energy while the two low temperature phases would have large absolute values of magnetization, one positive and the other negative, and small bond energies. This situation is exactly what is seen in roughly 10$`\%`$ of realizations. For systems much smaller than $`24^3`$, three phase coexistence is not seen simply because the magnetization distribution is too broad to display three distinct peaks. One hypothesis is that, as system size becomes large, the fraction of systems displaying three phases increases and that, in the thermodynamic limit, the transition at $`\mathrm{\Delta }=.35`$ is first-order. However, the existence of realization 14 with its five peaks is difficult to reconcile with this viewpoint.

If the transition is continuous and follows the droplet model scenario  we would expect two peaks in the magnetization histogram corresponding to two states of a single critical phase. The width of each peak should scale as $`L^{\gamma /\nu }`$ and the separation between the peaks should scale as $`L^{\beta /\nu }`$. The majority of our disorder realizations display two peaks and an alternative to the first-order transition hypothesis is that, as system size increases, the fraction of realizations with two peaks in the magnetization histogram approaches $`100\%`$. Although the original droplet model envisioned two states with roughly the same energy and differing by flipping a critical cluster of spins, it is possible to imagine more complicated droplet pictures where there are sometimes more than two states in the critical phase and therefore more than two peaks in the magnetization histogram. However, our five realizations that display more than one peak in the bond energy histogram are difficult to reconcile with the droplet model since one expects the states, however many there are, to be nearly degenerate in bond energy.

Ultimately, the essential difference between the first-order scenario and the droplet model continuous-transition scenario is whether there are distinct phases at the transition point or whether the peaks in the magnetization histogram are different states that are part of a single critical phase. The comes down to the question of whether the peaks in magnetization histogram move toward the origin as $`L\mathrm{}`$ or remain at finite values. Unfortunately, the consensus is that the magnetization exponent, if it is not actually zero, is very tiny so that the decrease in magnetization with system size would be far too weak to observe in simulations.

## V Conclusions

We have studied the phase transition of the random field Ising model by Monte Carlo simulation of systems of size $`24^3`$ using a new cluster algorithm. We find major qualitative differences between the behavior of different realizations of disorder. Some realizations display a broad two-peaked magnetization histogram consistent with a continuous transition, while a small fraction display a three-peaked structure consistent with a first-order transition. Our main conclusion is that more work needs to be done to determine the nature of the transition. It would be very useful to study larger system sizes and stronger disorder to determine whether there are trends in the qualitative features of the phase transition. Does the fraction of “first-order” systems increase as disorder or system size increases? It will require improvements of the algorithm including the fine-tuning technique or substantially more computer power to resolve these questions.

## VI Acknowledgements

We thank Pozen Wong and Robert Swendsen for useful discussions. This work was supported by NSF grants DMR-9978233 (JM) and DMS-9971016 (LC) and NSA grant MDA904-98-1-0518 (LC). JM thanks the Santa Fe Institute for hospitality during a visit in which some of this work was carried out.


text:
# Monte Carlo Simulations of Star Clusters

## 1. Introduction

Very detailed, recent observations of globular clusters suggest very close interplay between stellar evolution, binary evolution and dynamical interactions. This interplay is far from being understood. Monte Carlo codes, which use a statistical method of solving the Fokker–Planck equation provide all the necessary flexibility to disentangle the mutual interaction between physical processes important during globular cluster evolution. The codes were developed by Spitzer (1975, and references therein) and Hénon (1975, and references therein) in the early seventies, and substantially improved by Marchant & Shapiro (1980, and references therein) and Stodółkiewicz (1986, and references therein) and recently reintroduced by Giersz (1998, 2000), Heggie et al. 1999, Joshi et al. (1999a), Joshi et al. (1999b, hereafter JNR) and Rasio (2000, this volume). The Monte Carlo scheme takes full advantage of the undisputed physical knowledge of the secular evolution of (spherical) star clusters as inferred from continuum model simulations. Additionally, it describes in a proper way the graininess of the gravitational field and the stochasticity of the real $`N`$–body systems and provides as detailed as in direct $`N`$–body simulations information about movement of any objects in the system. This does not include any additional physical approximations or assumptions which are common in Fokker–Planck and gas models (for e.g. conductivity). Because of this, the Monte Carlo scheme can be regarded as a method which lies between direct $`N`$–body and Fokker–Planck models and combines most of their advantages. Moreover, Monte Carlo codes are simple, very fast, easily parallelized and easily scalable to the physical units. There is no need for special hardware or supercomputers to efficiently simulate evolution of realistic star clusters . However, as any numerical method, the Monte Carlo method suffers from some disadvantages. It can only deal with spherically symmetrical systems, and only small–angle two–body interactions. The galactic tidal field can only be approximated by the tidal cut-off, and unfortunately, cross–sections for some physical processes are needed (e.g. three–body binary formation). Additionally, the physical processes, which evolve on time–scales comparable to the crossing time–scale, can not be properly investigated, and the method has difficulty with the proper definition of local parameters (e.g. density, velocity dispersion). Despite all of these disadvantages the Monte Carlo method can be easily and efficiently used to simulate evolution of realistic globular clusters. The comparison between numerical simulations and observations will help to infer the initial parameters of proto–clusters and help to disentangle the interplay between physical processes involved in cluster evolution. Moreover, the Monte Carlo method can be used to simulate dynamical formation of massive black holes in dense spherical stellar systems (e.g. galactic nuclei).

## 2. Results

The Monte Carlo code is described in detail in Giersz (1998), which deals with simulations of isolated single–mass systems. Here and in Giersz (2000) the Monte Carlo code is extended to include the following additional physical processes:

* multi–mass systems described by the power–law initial mass function:
$`N(m)dm=Cm^\alpha dm,m_{min}mm_{max}`$, where $`C`$ and $`\alpha `$ are constants.
* stellar evolution introduced according to prescription given by Chernoff & Weinberg (1990, hereafter CW) or Taut et al. (1997).
* three–body binaries described by the suitably modified Spitzer’s formula (Spitzer 1987, Giersz 2000).
* binary–binary interactions introduced according to Mikkola (1983, 1984) and Stodółkiewicz (1986).
* tidal field simulated by tidal cut-off with energy and/or apocenter criterion.

The results of Monte Carlo simulations of star cluster evolution will be presented in the next two subsections.

### 2.1. Family $``$ 1

The initial conditions were chosen in a similar way as in a collaborative experiment (Heggie et al. 1999). The positions and velocities of all stars were drawn from a King model. All standard models have the same total mass $`M=60000M_{}`$ and the same tidal radius $`R_c=30`$ pc. Masses are drawn from the power–law mass function described above. The minimum mass was chosen as $`0.1M_{}`$ and maximum mass as $`15M_{}`$. Three different values of the power–law index were adopted: $`\alpha =1.5`$, $`2.35`$ and $`3.5`$. The set of initial King models was characterized by $`W_0=3`$, $`5`$ and $`7`$. Additional models of CW’s Family 1 were computed to facilitate comparison with results of other simulations (minimum mass equal to $`0.4M_{}`$, $`\alpha =1.5`$, $`2.5`$, $`3.5`$ and total mass $`M=90685M_{}`$, $`99100M_{}`$, $`103040M_{}`$, respectively).

In Table 1 the comparison between available results of $`N`$-body, Fokker–Planck and Monte Carlo simulations is presented. The standard models show a remarkably good agreement with $`N`$-body results (Heggie 2000). See columns labeled by G and H-0.1 in Table 1. Only models with a flat mass function show some disagreement. These models are difficult for both methods. Violent stellar evolution and induced strong tidal striping lead to troubles with time–scaling for the $`N`$-body model and proper determination of the tidal radius for the Monte Carlo model. Generally, the same is true for Monte Carlo models of Family 1. Results of these models show good agreement with results of CW, Aarseth & Heggie (1998) and Takahashi & Portegies Zwart (1999). JNR’s results, particularly for strongly concentrated systems, disagree with all other models. This can be connected with the fact, that JNR’s Monte Carlo scheme is not particularly suitable for high central density and strong density contrast. Too large deflection angles adopted by JNR and consequently too large time-steps can lead to too fast evolution in these models.

All standard models, for which mass loss due to violent stellar evolution of the most massive stars does not induce quick cluster disruption, evolve in a very similar way. The rate of mass loss, evolution of the central potential, evolution of the average mass does not depend much on the initial central concentration of the system. They depend strongly on the index of the mass function. Models of Family 1, on the contrary show, as well, dependence on the initial concentration. Very high initial mass loss across the tidal boundary, connected with evolution of the most massive stars (for models of Family 1, there are more massive stars than for standard models), forces the system to substantial changes of its structure and in consequence to different evolution of the total mass, anisotropy, etc.. Models which are quickly disrupted show only small signs of mass segregation. Models with larger central concentration survive the phase of rapid mass loss and then undergo core collapse and subsequent post–collapse expansion in a manner similar to isolated models. The expansion phase is eventually reversed when tidal limitation becomes important. As in isolated models, mass segregation substantially slows down by the end of the core collapse. After a core bounce there is a substantial increase in the mean mass in the middle and outer parts of the system, caused by the preferential escape of stars of low mass and tidal effects. Standard models, which are not quickly disrupted, show modest initial build up of anisotropy in the outer parts of the system. As the tidal stripping exposes inner parts of the system, anisotropy gradually decreases and eventually becomes slightly negative. The central part of the system stays nearly isotropic. Models of Family 1, from the very beginning, develop in the outer parts of the system modest negative anisotropy. It stays negative until the time of cluster disruption, when it becomes slightly positive (during cluster disruption most stars are on radial orbits).

### 2.2. 1,000,000 body run

For the first time, the Monte Carlo simulation, on the ”star–by–star” bases, of evolution of 1,000,000 body star cluster is presented. The initial conditions were as follows: total mass equal to $`319,305M_{}`$, tidal radius equal to $`33.57pc`$, power–low index of mass function equal to $`2.35`$, minimum and maximum mass equal to $`0.1M_{}`$ and $`15.0M_{}`$, respectively and King model parameter $`W_0=5`$.

The 1,000,000 body run shows basically the same features as, discussed above, models of Family $``$ 1. As an example of the overall cluster evolution the time dependence of Lagrangian radii, core radius and tidal radius are presented in Figure 1. The three different phases of evolution can be clearly distinguished. First, short phase of violent mass loss due to stellar evolution leads to overall cluster expansion. Even the innermost Lagrangian radius expands, the contraction connected with mass segregation is not strong enough to dominate the expansion. Second phase is characterized by the slow core collapse. Tidal effects are small and cluster behaves in a similar way as an ordinary isolated system. Then in the third phase, post–collapse evolution is superposed with growing tidal striping effects. The cluster nearly homogeneously contracts. The central parts of the system show clear signs of the gravothermal oscillations.

In Figure 2 the density profiles for the different epochs are presented. It is clear, that in the central parts of the system the density profile shows steeper slope than $`2.2`$ (line labeled by 6), the standard value for single–mass systems. This is in agreement with results of CW. The power–low index is a function of the ratio of mass of the most massive stars to the average mass. The larger the ratio the smaller the power–law index. The core is mainly populated by the most massive stars (massive white dwarfs, neutron stars and black holes), whose masses are larger than the average mass in the vicinity of the core. So as a consequence, the power–law index is smaller than $`2.2`$. The density profiles for an advanced collapse phase (lines labeled by 4, 5 and 6) show a bump in the middle part of the system (close to $`0.08`$ in x–axis). The bump originates because of the growing influence of low mass stars on the determination of the local density. The position and the size of the bump is in a good agreement with Fokker–Planck results (Takahashi & Lee 2000).

In order to perform simulations of real globular clusters several additional physical effects have to be included into the code. The tidal shock heating of the cluster due to passages through the Galactic disk, interaction with the bulge, shock–induced relaxation, primordial binaries, physical collisions between single stars and binaries are one of them. Inclusion of all these processes do not pose a fundamental theoretical or technical challenge. It will allow to perform detailed comparison between simulations and observed properties of globular clusters and will help to understand the globular cluster formation conditions and explain how peculiar objects observed in clusters can be formed. These kinds of simulations will also help to introduce, in a proper way, into future $`N`$–body simulations all necessary processes to simulate on the star–by–star basis evolution of real globular clusters from their birth to death.

#### Acknowledgments.

I would like to thank Douglas C. Heggie for making the $`N`$–body results of standard model simulations available.

## References

Aarseth S.L. & Heggie D.C., 1998, MNRAS, 297, 794

Chernoff D.F. & Weinberg M.D., 1990, (CW), ApJ, 351, 121

Giersz M., 1998, MNRAS, 298, 1239

Giersz M., 2000, MNRAS, submitted

Heggie D.C., 2000, private information

Heggie D.C., Giersz M., Spurzem R., Takakashi K., 1999, in Highlights of Astronomy, vol 11, ed. J. Andersen, 591

Hénon M., 1975, in Hayli A., ed., Dynamics of Stellar Systems, Reidel: Dordrecht, 133

Joshi K.,J., Nave C. & Rasio F.A., 1999a, ApJ submitted (astro-ph/9912155)

Joshi K.,J., Rasio F.A. & Portegies Zwart S., 1999b, (JNR), ApJ submitted (astro-ph/9909115)

Marchant A.B. & Shapiro S.L., 1980, ApJ, 239, 685

Mikkola S., 1983, MNRAS, 205, 733

Mikkola S., 1984, MNRAS, 208, 75

Rasio F.A., 2000, in Spurzem R., et al., eds, Dynamics of Star Clusters and the Milky Way, (ASP Conference Series)

Spitzer L., Jr., 1975, in Hayli A., ed, Dynamics of Stellar Systems, Reidel: Dordrecht, p.3

Spitzer L., Jr., 1987, Dynamical Evolution of Globular Clusters. Princeton Univ. Press, Princeton, 77

Stodółkiewicz J.S., 1986, Acta Astr., 36, 19

Takahashi K. & Lee H.M., 1999, MNRAS submitted (astro-ph/9909006)

Takahashi K. & Portegies Zwart, S. F., 1999, ApJ submitted (astro-ph/9903366)

Tout C.A, Aarseth S.J., Pols O.R. & Eggleton P.P, 1997, MNRAS, 291, 732


text:
# Confidence levels of evolutionary synthesis models

## 1 Introduction

Since the development of evolutionary synthesis models by Tinsley (1980) several groups have developed models trying to reproduce observable quantities in systems in which the stellar content is not resolved. These models provide a powerful tool to study galactic and extragalactic Hii regions. Several types of models can be found in the literature. For instance, Leitherer et al. (1996) give an extensive and comprehensive review on them.

The variety of external inputs used in the models (evolutionary tracks, stellar atmosphere libraries, etc) and the treatment of these input parameters, (for example, the interpolation methods), lead to slightly different model outputs from different authors, although “grosso modo” all of them produce similar results. One way to check the applicability of these models is to compare the results using different inputs. In this direction we highlight the work of Bruzual (2000) in which the author presents an extensive study on how the external inputs to the models influence the output results.

These problems can be exacerbated when modelling systems with small number of stars given that, in general, synthesis models use continuous functions that do not reproduce exactly the discontinuous nature of star formation, especially in systems with low number of stars. In this letter we address this issue by studying the intrinsic variations of the evolutionary synthesis models outputs due to the deviations of the analytical Initial Mass Function (IMF) when modelling Hii regions.

Here, we study the importance of such fluctuations in some typical observables of extragalactic Hii regions. We characterise our results showing confidence levels of the distribution of a given parameter as a function of age and system mass. We focus on the study of the H$`\beta `$ equivalent width in emission (EW(H$`\beta `$), that is extensively used as an age indicator of star forming regions. The results for the number of ionizing photons, $`Q(H^0)`$, and the ratio of the blue Wolf-Rayet (WR) bump luminosity to the H$`\beta `$ luminosity are also presented. We differ the study of other quantities using different metallicities, IMF slopes, star formation regimes and evolutionary tracks to a future paper.

## 2 Synthesis code and the Initial Mass Function

The study of the IMF has been broadly covered in the astronomical literature. See Scalo (1986) for a comprehensive review. We define the IMF in the following way:

$$\mathrm{\Phi }(m)=\frac{dN}{dm}=Am^\alpha $$
(1)
where $`\alpha `$ is the IMF slope, and $`A`$ is a normalisation factor. This function gives us the number of stars in a given mass range. The widely used Salpeter IMF corresponds to $`\alpha =2.35`$ with this definition. The total mass of the system will be

$$M=_{m_{\mathrm{low}}}^{m_{\mathrm{up}}}m\mathrm{\Phi }(m)𝑑m=_{m_{\mathrm{low}}}^{m_{\mathrm{up}}}m𝑑N$$
(2)
In an evolutionary synthesis code, the amount of stars produced is usually generated binning the IMF (or using Montecarlo simulations). The evolution of each bin (or star) is followed along the corresponding evolutionary track and the final result is computed from the overall evolution of all the stars (or masses) considered. In such a procedure, variations in the IMF affect directly the output results. The IMF and the synthesised parameters are usually normalised to the total mass transformed into stars in the synthetic cluster. The comparison with real data is always carried out comparing an observed quantity with the corresponding one obtained with the code. It is assumed that the synthesised quantity is a good approximation to the theoretical analytical value and usually no further considerations are taken into account when comparing data and model outputs.

We have used an updated version of the evolutionary synthesis code of Cerviño & Mas-Hesse (1994), that uses Montecarlo realisations of the IMF to produce star forming regions of a given mass. We have also used the code to generate analytical simulations of clusters of the same mass, where analytical means a mass binning approximation of the IMF. We have tested that the analytical and the Montecarlo outputs show consistent results. The dispersion between them is lower than 1.5% when a large number of stars is used ($`5\times 10^5`$ stars, i.e. $`3\times 10^6`$ M transformed into stars in the mass range 2 – 120 M). In this version of the code we use the solar metallicity evolutionary tracks of Meynet et al. (1994). The analytical results have also been compared with the predictions of the code of Leitherer et al. (1999), that uses the same evolutionary tracks. The results obtained with both codes are in general agreement.

We have performed 600 Montecarlo realisations of a cluster in which 103 M are transformed into stars in the mass range 2 – 120 M; 400 realisations of a cluster of 104 M and 200 of a cluster with 105 M. The synthesis code is the same in all Montecarlo realisations, the only difference being the stochasticity of the Montecarlo process which reflects the fluctuations of the actual mass distribution with respect to the analytical continuous IMF.

## 3 Distribution of observed parameters: Confidence levels

In this section we discuss the resulting distributions of the observable parameters obtained from the different Montecarlo runs for the masses mentioned above. We note that these distribution are in general non-Gaussian. The confidence levels presented hereafter are computed integrating the resulting distributions to contain the 90% (68%) central values of the given parameter. We refer to these values of the parameters as the 90% (68%) confidence levels of the resulting parameter distribution obtained from the Montecarlo realisations.

Fig. 1. shows the 90% and 68% (equivalent to 1$`\sigma `$) confidence levels for EW(H$`\beta `$) resulting from the simulations in which $`10^4`$ M are transformed into stars. This is the mass range in which Hii regions with L(H$`\alpha `$) $``$ $`10^{39}`$ erg s-1 lie. This value of luminosity of H$`\alpha `$ is usually used as the boundary between giant and normal Hii regions.

It can be seen from the figure that the confidence levels change with time. The largest spread in EW(H$`\beta `$) is obtained between 3 and 9 Myrs (see also fig. 3) when the fluctuations in the ionizing flux and the continuum can lead to variations of the order of 20% at the 1$`\sigma `$ level. The dispersion is reduced in the age range between 9 and 13 Myrs and increases again to an almost constant level for older ages.

The mass dependence of such fluctuations can be seen in Fig. 2, in which the 90% and 68% confidence levels are displayed for a fixed age of 4.5 Myr since the onset of the burst as a function of the mass transformed into stars. The figure shows clearly that the dispersion is not symmetric around the analytical value. The lower confidence level shows an almost linear correlation with the logarithm of the mass transformed into stars. The upper level does not show this linear behavior. An important implication regards the different interpretations of observational data when compared to models assuming that they provide an exact value. As can be seen from the figures, the 90% or 68% confidence levels of the parameter distributions are uncomfortably large, especially for low masses. Due to this behaviour of the EW(H$`\beta `$) evolution, the age determination between 3 and 5 Myrs are not well defined at the 90% confidence level. This situation may be improved using more diagnostic observables to estimate the age since the onset of the burst.

Fig. 3 presents the evolution with age of the 90% confidence level of the EW(H$`\beta `$) distribution for different amount of mass transformed into stars. This figure demonstrates the importance of taking into account the confidence levels of an observable parameter for a given mass range. For instance, for a 105 M cluster, a value of the EW(H$`\beta `$) equal to 50Å, is compatible with any age between 3.5 and 5 Myrs within the 90% confidence levels.

In Fig. 4 the evolution of the ionizing flux, $`Q(H^0)`$, is displayed against the mass transformed into stars. For a 103 M system the fluctuations allow values that differ by an order of magnitude for a given age. The 90% confidence level values for the distribution of this observable show a relatively broad distribution up to 7 Myrs. Comparing this figure with the previous one, it is possible to figure out that the the scatter in the continuum flux around 4861 Å, $`L_c(\mathrm{H}\beta )`$, is the driving parameter producing the time evolution of the EW(H$`\beta `$) distribution. We have checked this point with the actual values and confirmed it to be the case. As a consistency check, we have also corroborated that the dispersion in $`L_c(H\beta )`$ combined with the dispersion in $`Q(H^0)`$ approximately reproduce the dispersion in EW(H$`\beta `$).

Fig. 5 shows the evolution of the ratio of the Wolf-Rayet bump and H$`\beta `$ luminosities, $`L(\mathrm{WRbump})/L(H\beta )`$, presenting the 90% confidence levels of its distribution. We would like to point out that the analytical value, close to 1, is 2.3 times higher that the one obtained in Cerviño & Mas-Hesse (1994). We have checked this value and confirm that it is the result of using different evolutionary tracks with enhanced mass loss rates. It is important to note that at the 90% confidence level, for a 103 M cluster, values consistent with no detection of WR stars are allowed. The distribution at 1$`\sigma `$ level is also compatible with such non detection.

This study opens an additional question. If stochasticity of the stellar mass distribution of low mass clusters imply such dispersions in the distribution of observable parameters, what would its effect be in a system composed of several smaller groups? In order to address this question, we have simulated multiply-composed 105 M clusters by adding randomly 6 sets of 100 clusters with 103 M and 6 sets of 10 clusters 104 M to obtain the observables for a variety of ages. We have plotted the results for the EW(H$`\beta `$) and the $`L(\mathrm{WRbump})/L(H\beta )`$ together with the analytical values and the 90% confidence level of the distributions for a 105 M cluster in Fig. 6.

The figure shows that all the simulated points fall within the 90% confidence level for the 105 M simulation. However, this is not the case for the 68% confidence level. In conclusion, the distribution of the observed parameters due to the stochasticity of the star formation process shrinks towards the analytical values of evolutionary synthesis models as larger masses are considered more or less independently of whether these systems are monolithic unique systems or composed by several coeval subgroups. We would like to stress that our previous simulations assumed the same age for the subgroups that form a system. We have not addressed the consequences of different ages for the different subgroups.

As a final note we would like to point out that the total masses used here have been computed using stellar masses in the range from 2 to 120 M, which are the ones relevant for the observable parameters discussed. The corresponding total masses should be corrected by multiplying by a factor 1.36 for stellar masses in the range from 1 to 120 M; multiplying by 1.5 for the mass range 0.8 to 120 M and multiplying by 3.75 for the mass range 0.08 to 120 M. In all the cases a power law IMF with the Salpeter exponent has been assumed.

## 4 Conclusions

The analysis presented in this letter allow us to draw the following conclusions:

* The intrinsic deviations between analytical and stochastic IMF must be taken into account when comparing synthesis models with observational data. The inherent uncertainties depend on the amount of gas transformed into stars.
* These deviations are independent of the synthesis code and represent a lower limit to the total uncertainties. The total deviations will depend on other parameters like the star formation history, the input ingredients, the numerical approximations, etc.
* The deviations depend on the evolutionary tracks used.
* The widths of the parameter distributions compared to the analytical values are proportional to the mass transformed into stars in the stellar cluster or groups of clusters. They also depend on the IMF slope and the evolutionary status of the cluster.
* The specific distribution of the deviation from the analytical value varies from observable to observable. As an extreme example, the 68% confidence level of the $`L(\mathrm{WRbump})/L(H\beta )`$ ratio is compatible with no detection of WR stars.

###### Acknowledgements.

We want to acknowledge useful discussion with M. Gómez-Flechoso, J. Knödlseder, J.F. Gómez, J.M. Mas-Hesse and D. Schaerer. VL also acknowledges the Observatoire Midi-Pyrénées for providing facilities to conduct this study.


text:
# Slow logarithmic relaxation in models with hierarchically constrained dynamics

## Abstract

A general kind of models with hierarchically constrained dynamics is shown to exhibit logarithmic anomalous relaxation, similarly to a variety of complex strongly interacting materials. The logarithmic behavior describes most of the decay of the response function.

In the last years, an anomalous slow relaxation that can be accurately described by a logarithmic law, has been found in the time evolution of a variety of complex strongly interacting materials. This includes spin-glasses , granular materials , structural glasses , and protein models . In all the above cases, the decay of some perturbation $`\varphi (t)`$ can be fitted at intermediate times by an expression of the form

$$\varphi (t)=C_1C_2\mathrm{ln}t,$$
(1)
where $`C_1`$ and $`C_2`$ are constant, and may depend on the macroscopic variables (temperature, compactivity, …) characterizing the state of the system. The time range in which the logarithmic law accurately describes the relaxation covers most of the change of $`\varphi (t)`$ from its initial value to the final one.

Palmer, Stein, Abrahams, and Anderson (PSAA) in a seminal paper  described some dynamical models for relaxation in complex systems with strong interactions. The models were based in a hierarchically constrained dynamics so that the relaxation of the system involves a sequential series of correlated processes. Although PSAA concentrate in those physical systems showing “stretched exponential” relaxation, it is clear that the general picture emerging from hierarchically constrained dynamics can be translated into very different physical situations. As a consequence, the above models provide a very useful tool to investigate many kinds of “anomalous relaxation” in complex systems.

In particular, the picture of hierarchical constrained dynamics seem to be physically appropriate to describe many of the complex systems showing logarithmic relaxation. Let us consider as a particular example the densification of powders and structural glasses at high hydrostatic pressure studied in Refs. . External pressure causes some selected mesoscopic regions of the system to reach a more compact state very quickly. The nearest neighbor regions are influenced by this rearrangement, and can also reorganize themselves. Of course, these new structural changes will affect the regions outside the nearest neighbors, allowing them to relax, and so on. The above discussion suggests that the relaxation in this kind of systems is slow because it consists of a large number of correlated reorganization steps with an increasing characteristic time. In other words, it has a hierarchical structure, with the faster degrees of freedom (associated to the “inner” regions) controlling the relaxation of the slower degrees of freedom (associated to the “outer” regions). A similar picture has been already suggested in Refs. . The aim of this paper is to introduce a quite general class of hierarchical models, which seem appropriate to describe this kind of systems, in order to derive the logarithmic relaxation of Eq. (1).

Like PSAA we will formulate our model in the most general way, hoping it characterizes the behavior of a specific class of systems. We assume that the state of the system can be given by means of $`N`$ modes. They do not correspond to a particle description of the system, but to some mesoscopic level of description. In this sense, $`N`$ does not scale with the size of the system, but it is assumed to have a definite value, independent of the size, for a large enough (macroscopic) system. Of course, the specific value depends on the system and the physical situation we are dealing with.

The $`N`$ modes can be grouped into a discrete series of levels, $`n=0,1,2,\mathrm{}`$. The number of modes belonging to level $`n`$ will be denoted by $`N_n`$, so that $`N=_nN_n`$. Moreover, modes in level $`n`$ will be indicated by the variables $`S_i^{(n)}`$, $`i=1,2,\mathrm{},N_n`$. For the densification experiments of Refs. , we can regard $`S_i^{(n)}`$ as a pseudospin variable, with its two possible values corresponding to different volumes of the mesoscopic regions. The dynamics of the system is defined in a hierarchical way as follows. A mode in level $`n+1`$ can only modify its state if the configuration of modes in level $`n`$ belongs to a well defined subset of all the possible configurations of the modes in that level. Clearly, this slows down the modes in a given level, as compared with those of the previous level. More specifically, the average transition rate $`W_{n+1}`$ for a mode in level $`n+1`$ is assumed to have the form

$$W_{n+1}=W_nP_n,$$
(2)
where $`P_n`$ is the probability that the modes in level $`n`$ are in any of the configurations allowing to relax level $`n+1`$. Our choice for $`P_n`$ will be

$$P_n=\mathrm{exp}(N_n\mathrm{\Delta }\mu ),$$
(3)
with $`\mathrm{\Delta }\mu `$ being a dimensionless free energy defining the activation barrier per mode. Eq. (3) appears quite naturally if $`P_n`$ is evaluated from the equilibrium distribution of the system under consideration, provided that the subset of configurations of level $`n`$ allowing to relax level $`n+1`$ involve all the modes (or an extensive part of them) in level $`n`$. Of course, the free energy $`\mathrm{\Delta }\mu `$ must be appropriately interpreted in each particular situation. For instance, in the densification experiments at high pressure  it would be an activation enthalpy divided by the temperature of the sample. On the other hand, in a granular system, it could be understood as an “effective volume” divided by the compactivity . Thus, the corresponding average relaxation time of the mode is

$$\tau _{n+1}=\tau _n\mathrm{exp}(N_n\mathrm{\Delta }\mu ).$$
(4)
For the sake of simplicity, the barrier $`\mathrm{\Delta }\mu `$ is assumed to be independent of the level $`n`$ of the mode and also of the number of modes in the level $`N_n`$. Iteration of Eq. (4) yields

$$\tau _n=\tau _0\mathrm{exp}\left(\mathrm{\Delta }\mu \underset{m=0}{\overset{n1}{}}N_m\right).$$
(5)
The range of variation of $`n`$ and the values of the populations $`N_n`$ characterize the complexity of the time evolution of the system. In principle, there is no need for a maximum value $`n_{\text{max}}`$ of $`n`$, i.e. $`n_{\text{max}}`$ may be infinite. On the other hand, the systems we are interested in show slow relaxation, which requires a large relaxation time. In this way, it is necessary that

$$\tau _{\text{max}}=\underset{nn_{\text{max}}}{lim}\tau _n\tau _0.$$
(6)
This condition for slow relaxation was already pointed out by PSAA .

We are interested in the relaxation of the quantity

$$Q(t)=\frac{1}{N}\underset{n}{}\underset{i=1}{\overset{N_n}{}}S_i^{(n)}(t),$$
(7)
where the angular brackets denote configuration average. Note that in the context of the densification experiments  $`Q(t)`$ would be proportional to the actual volume of the sample at time $`t`$, with the physical meaning given above to the variables $`S_i^{(n)}`$. The corresponding relaxation function is defined in terms of the asymptotic long time value $`Q(\mathrm{})`$ as

$$\varphi (t)=\frac{Q(t)Q(\mathrm{})}{Q(0)Q(\mathrm{})}=\underset{n}{}\omega _n\mathrm{exp}\left(\frac{t}{\tau _n}\right),$$
(8)
with $`\omega _n=N_n/N`$ being the fraction of modes in level $`n`$. Combination of Eqs. (5) and (8) gives

$$\varphi (t)=\underset{n=0}{\overset{n=n_{\text{max}}}{}}\omega _n\mathrm{exp}\left[\frac{t}{\tau _0}\mathrm{exp}\left(\zeta \underset{m=0}{\overset{n1}{}}\omega _m\right)\right],$$
(9)
where

$$\zeta =N\mathrm{\Delta }\mu .$$
(10)
Next we will assume that $`\omega _n`$ changes very smoothly with $`n`$. This allows to introduce a continuum limit for the distribution of relaxation times, which is expected to be closer to real systems than the discrete level model considered up to now. We define a function $`\stackrel{~}{\omega }(x)`$ by $`\stackrel{~}{\omega }(nϵ)=\omega _n`$, with $`ϵ1`$, so that the variable $`x_n=nϵ`$ is almost continuum and the sums over $`n`$ can be replaced by integrals. Let us introduce the new variable ($`x_{\text{max}}=n_{\text{max}}ϵ`$)

$$u=\frac{_0^x𝑑x^{}\stackrel{~}{\omega }(x^{})}{_0^{x_{\text{max}}}𝑑x^{}\stackrel{~}{\omega }(x^{})},$$
(11)
$`0u1`$, which measures the fraction of the total number of modes $`N`$ which are in levels up to a value $`n=x/ϵ`$. It is then a simple task to transform Eq. (9) into

$$\varphi (t)=_0^1𝑑u\mathrm{exp}\left(\frac{t}{\tau _0}e^{\zeta u}\right),$$
(12)
showing that the decay of $`\varphi `$ can be characterized by only two parameters, namely $`\tau _0`$ and $`\zeta `$. The former fixes the time scale, while the latter determines the actual shape of the relaxation. The elementary relaxation times given by Eq. (5) are $`\tau (u)=\tau _0\mathrm{exp}(\zeta u)`$, and, in particular, $`\tau _{\text{max}}=\tau _0\mathrm{exp}\zeta `$. Thus, the condition for slow relaxation, Eq. (6), implies that $`\mathrm{exp}\zeta 1`$.

The analysis of Eq. (12) is greatly simplified by realizing that it is equivalent to

$$\varphi (t)=\frac{1}{\zeta }\left[E_1\left(\frac{t}{\tau _0}e^\zeta \right)E_1\left(\frac{t}{\tau _0}\right)\right],$$
(13)
where $`E_1(z)`$ is the integral exponential function . There are two natural time scales in Eq. (13). One is that defined by $`t/\tau _0`$, while the other is a slow time scale defined by $`s=\frac{t}{\tau _0}e^\zeta \frac{t}{\tau _0}`$. In the following, we will use the asymptotic expansions of the exponential integral function from Ref. . For very short times such that $`t/\tau _01`$, Eq. (13) is equal to

$$\varphi (t)=1\frac{1e^\zeta }{\zeta }\frac{t}{t_0}+𝒪\left(\frac{t}{\tau _0}\right)^2,$$
(14)
i.e., $`\varphi (t)`$ decays very little from its initial unity value in this region. In the opposite limit of very large times for which $`s1`$, it is

$$\varphi (t)\frac{1}{\zeta }\frac{e^s}{s},$$
(15)
the response function has already decayed to very small values. One main conclusion of the above discussion is that the relevant part of the relaxation, in which $`\varphi (t)`$ changes from values very close to unity to very close to zero, takes place between the two asymptotic time regimes we have analyzed. In this way, we are led to to consider the distinguished limit $`t/\tau _01`$, $`s1`$. In this intermediate time window it is found that

$$\varphi (t)1\frac{1}{\zeta }\left(\gamma +\mathrm{ln}\frac{t}{\tau _0}\right).$$
(16)
This is precisely the logarithmic decay of Eq. (1). Moreover, the previous discussion of the short and large time limits, clearly indicates that most of relaxation of the system is expected to be accurately described by Eq. (16). As an example, we present in Fig. 1 the relaxation function $`\varphi (t)`$ for $`\mathrm{exp}\zeta =10^3`$ and $`\mathrm{exp}\zeta =10^5`$. The integral expression given by Eq. (12) (circles) is well fitted by the logarithmic asymptotic form of Eq. (16) (solid lines). The agreement extends over several decades of intermediates times, measured in units of the basic time $`\tau _0`$. The fit improves as the value of $`\zeta `$ increases, consistently with the condition for time scales separation.

It is also interesting to compute an average relaxation time $`\overline{\tau }`$, given by

$$\overline{\tau }=_0^{\mathrm{}}𝑑t\varphi (t)=\tau _0\frac{e^\zeta 1}{\zeta }\tau _0\frac{e^\zeta }{\zeta },$$
(17)
showing an exponential increase with $`\zeta `$. Taking into account the definition of $`\zeta `$, Eq. (10), the above result can be understood as corresponding to an activation energy barrier whose height is proportional to the complexity of the system, measured by the total number of modes $`N`$ actively involved in the relaxation. As already mentioned, for a given specific situation, the dimensionless parameter $`\mathrm{\Delta }\mu `$ might have a well defined physical interpretation, being for instance the ratio between the Gibbs free energy barrier per mode and the thermal energy. In this way, it is also interesting to note that the coefficient of $`\mathrm{ln}t`$ in Eq. (16) is $`\zeta ^1`$, which is then proportional to the temperature $`T`$. This would lead to the $`T\mathrm{ln}t`$ behavior observed in many physical systems, for instance in Refs. .

In summary, we have studied a simple class of models which exhibit an anomalous logarithmic relaxation. The dynamics is hierarchically constrained, and they are adequate to study slow relaxation because there is a clear separation of time scales. The slow logarithmic relaxation is similar to what is experimentally observed in a variety of complex systems, fitting the data in most of the relevant part of the relaxation, i.e., the time window where the relaxation function is neither too small nor too close to its initial value. Perhaps the main assumption of this work is the “extensiveness” with $`N_n`$ of the effective barrier for modes in level $`n+1`$, as expressed by Eq. (3). This is quite a general condition and, therefore, the class of hierarchical models presented here may provide an explanation of the appearance of linear logarithmic relaxation in very different physical situations. In this context, the role of hierarchically constrained dynamics may also be relevant in order to understand other kinds of anomalous relaxation, as the inverse logarithmic behavior found in compaction experiments of granular systems . Also, it seems interesting to investigate the behavior of hierarchical models in the temperature cycling experiments usually carried out in glassy systems .

###### Acknowledgements.

This research was partially supported by Grant No. PB98-1124 from the Dirección General de Investigación Científica y Técnica (Spain).


text:
# Untitled Document

###

The thermodynamics is extented to spacetimes with spin-torsion density.Impplications to Einstein-Cartan-de Sitter inflationary phases are discussed.A relation between the spin-torsion density,entropy and temperature is presented.A lower limit for the radius of the Universe may be obtained from the spin-torsion density and the Planck lenght.

## 1 Introduction

###

Riemannian extensions to the second law of thermodynamics to General Relativity have been considered by H.P. de Oliveira and D.Pavon  in recent years.Since the spin dynamics seems to be fundamental on the mechanism of statistical mechanics and themodynamics and since spin dynamics is reasonably well explained in the context of Einstein-Cartan theory of gravity,the recent interest of thermodynamics on the investigation of cosmic background microwave radiation (CBMR) distribution in the Universe motivated us to investigate here the thermodynamics in the context of some inflationary models.We first solve the Gibbs equation in the realm of Einstein-Cartan gravity for a closed cosmological model and solve it for the radiation and inflationary eras.The next step is to solve the Einstein-Cartan cosmological equations for the de Sitter inflationary case.In this example we obtain a type of Ising model like a relation between the spin-torsion density and the absolute temperature and entropy density,producing a coherent model for thermodynamics in cosmology where the temperature and entropy increases with the spin-torsion density.Early attempts to include thermodynamics into the scope of Einstein-Cartan gravity were given by R.Amorim and L.L.Smalley  by making use of a Lagrangean treatment.Nevertheless no attempts had been made by those authors to apply their theories to Cosmological problems.The COBE data  of the CMBR to compute the spin-torsion fluctuation since it is shown that the temperature fluctuation measured by COBE is proportional to the spin-torsion fluctuation.

###

## 2 Thermodynamics in Closed Einstein-Cartan Cosmological Models

Let us start by defining the entropy density $`s`$ as

$$s=\frac{\gamma \rho }{T}$$
(2.1)
From the conservation law in EC gravity

$$\dot{\rho }+3H(\rho +p)=2\pi (\dot{\sigma ^2}+3H\sigma ^2)$$
(2.2)
where we have addopted $`G=1`$ units and the equation of state is given by

$$p=(\gamma 1)\rho $$
(2.3)
where $`\gamma `$ is the adiabatic index.Substitution of equation (2.3) into (2.2) yields

$$\dot{\rho }+3H\gamma \rho =2\pi (\dot{\sigma ^2}+3H\sigma ^2)$$
(2.4)
From the conservation law (2.3) one obtains

$$\frac{d(\rho ^{}R^3)}{dt}=3p^{}R^3$$
(2.5)
where

$$\rho ^{}=\rho 2\pi \sigma ^2$$
(2.6)
and

$$p^{}=p2\pi \sigma ^2$$
(2.7)
are respectively the effective energy density and pressure including the spin-torsion density

$$\sigma ^2=<S_{ij}S^{ij}>=\frac{n^2h^2}{R^6}$$
(2.8)
where $`R(t)`$ is the cosmic scale factor and $`S_{ij}`$ is the spin-torsion density tensor.Here $`i,j=0,1,2,3`$ and $`n`$ is the nucleon number and $`h`$ is the Planck constant.From the Gibbs relation

$$T\frac{dS}{dt}=\frac{dE}{dt}+p\frac{d(2\pi ^2R^3)}{dt}$$
(2.9)
Here $`2\pi ^2R^3`$ is the volume of the spatially open $`k=+1`$ cosmological model  and

$$S=sR^3$$
(2.10)
where S is the total entropy and s is the entropy density.By making $`p^{}=0`$ into expression(2.5) one obtains

$$T\frac{d(sR^3)}{dt}=0$$
(2.11)
Substitution of (2.1) into (2.11) yields

$$\frac{d\frac{(\gamma \rho R^3)}{T}}{dt}=0$$
(2.12)
Expansion of expression (2.12) yields

$$\frac{\dot{\gamma }\rho R^3}{T}\frac{\dot{T}}{T}\gamma \rho R^3+\frac{\gamma \dot{\rho }R^3}{T}+3\frac{R^2\dot{R}\gamma \rho }{T}=0$$
(2.13)
Substitution of the conservation law into expression (2.13) yields

$$\frac{\dot{\gamma }}{\gamma }\frac{\dot{T}}{T}3H(\gamma 1)+2\pi [\dot{\sigma ^2}+3H\sigma ^2]=0$$
(2.14)
Taking the factor $`\gamma =\gamma _0=constant`$,after some algebra we obtain

$$lnT=3(1\gamma _0)R(t)+\frac{1}{\pi }\frac{dS_{0}^{}{}_{}{}^{2}}{R^2}$$
(2.15)
where here $`S_0^2=2\pi ^2\sigma ^2R^3`$ is the total spin-torsion quantity.Expression (2.15) yields

$$T(t)=e^{3(\gamma _01)R(t)+\frac{1}{\pi }{\scriptscriptstyle {\scriptscriptstyle \frac{dS_{0}^{}{}_{}{}^{2}}{R^2}}}}$$
(2.16)
Substitution of the relation

$$dS_0^2=\frac{6\pi n^2h^2dR}{R^4}$$
(2.17)
which upon integration and substitution on the last term on the RHS of equation(2.16) yields

$$T(t)=e^{3(\gamma _01)R(t)+\frac{30\pi n^2h^2}{R^7}}$$
(2.18)
which shows that the spatially closed Friedmann-Robertson-Walker model with spin-torsion density would cool down as the spin-torsion density increases, even for values of this density as big as the one corresponding to th nucleon number of the Universe $`n=10^{80}`$ nucleons.Another important physical result of the above result is that in the torsion free case as shown by Starkovich and Cooperstock  during the inflation era $`\gamma <\frac{2}{3}`$ ,$`\dot{T}>0`$ while during the radiation era $`\gamma =\frac{4}{3}`$ ,$`\dot{T}<0`$ and $`T`$ reaches a maximum at the beginning of the standard model era.Let us now examine what happens in the more general case of spacetimes with spin-torsion density.In this case we have

$$\frac{\dot{T}}{T}=2\pi [\dot{\sigma ^2}+3H\sigma ^2]+3H(1\gamma )$$
(2.19)
On the radiation era this equation reduces to

$$\frac{\dot{T}}{T}=2\pi [\dot{\sigma ^2}+3H\sigma ^2]H$$
(2.20)
Therefore the condition $`\dot{T}<0`$ together with equation (2.20) implies

$$[\dot{\sigma ^2}+3H\sigma ^2]<\frac{3H}{2\pi }$$
(2.21)
Substituting the value of the spin-torsion density in terms of the cosmic scale distance $`R(t)`$ yields $`R_{min}^{}{}_{}{}^{3}>nh=nc^3l_P^2`$ where $`l_P`$ is the Planck lenght.Therefore the spin-torsion density places a lower limit in the radius of the Universe avoiding the cosmological Singularity!As early as 1978 Kuchowicz  has obtained similar results withour recurring to a thermodynamical model.Using Friedmann-like models and spin=spin interactions he was able to obtain for example for an Euclidean model a minimum radius $`R_{min}=\frac{S_{0}^{}{}_{}{}^{2}}{4D}`$ where D is an integration constant.

## 3 Thermodynamics on de Sitter Inflationary Phase in Einstein-Cartan Gravity

###

In this section I shall consider an example of a de Sitter Cosmological model with spin-torsion density and thermodynamics.From the EC equations for this case one obtains

$$H^2=\frac{8\pi }{3}(\rho 2\pi \sigma ^2)$$
(3.1)
The other equation is given by

$$H^2=\frac{4\pi }{3}((\rho +3p)4\pi \sigma ^2)$$
(3.2)
Since $`p=(\gamma 1)\rho `$ from equation (3.2) one obtains

$$H^2=\frac{4\pi }{3}((3\gamma 2)\rho 4\pi \sigma ^2)$$
(3.3)
By equating equations (3.1) and (3.3) one obtains

$$\sigma ^2=\frac{\gamma \rho }{2\pi }$$
(3.4)
Substitution of expression (2.1) into (3.4) yields

$$\sigma ^2=\frac{sT}{2\pi }$$
(3.5)
This expression seems to be analogous to the Ising model for the behaviour of spins on a ferromagnet.Besides it shows that the increase on the spin-torsion density represents an increase on the temperature and entropy density.Expressing the spin-torsion density in terms of the cosmic scale $`R(t)`$ allow us to obtain an expression between the temperature total entropy and the cosmic scale factor given by

$$T=\frac{2\pi n^2h^2}{SR^3}$$
(3.6)
From the standard cosmology we know that

$$T=\frac{1}{R}$$
(3.7)
which yields

$$\frac{\delta T}{T}=\frac{\delta R}{R}$$
(3.8)
From the expression for $`\sigma ^2`$ one obtains

$$\delta \sigma ^2=6\frac{(nh)^2\delta R}{R^7}$$
(3.9)
which finally implies from the COBE data that

$$\frac{\delta T}{T}=\frac{1}{6}\frac{\delta \sigma ^2}{\sigma ^2}=10^5$$
(3.10)
which allow us to determine the spin-torsion density fluctuation from the COBE data .To better compare our thermodynamics to their ferromagnetic analogous Amorim-Smalley Lagrangean could be used to obtain a Hamiltonian in terms of the spin-torsion density and to compare the results with the Ising model.This can be done elsewhere. Acknowledgements: I am very much indebt to Prof.L.L.Smalley and H.P.de Oliveira for their interest in our work.Financial support from CNPq (Brazilian Government Agency) and UERJ are gratefully acknowledged.