format: text

text:
# Universal scaling functions for bond percolation on planar random and square lattices with multiple percolating clusters

## I Introduction

Percolation is related to many interesting scientific phenomena Stauffer94 . In recent years percolation problems with multiple percolating clusters have attracted much attention \[2-19\]. Most of the simulational studies of such problems have been restricted to percolation on lattices hw97note . However, many physical systems with multiple percolating clusters such as Carbino disks used in the study of quantum Hall effects Rusin96 , or oil fields confronted with drilling problems, do not have underlined regular lattice structures. Thus, it is of interest to know the relationship between the quantities for percolation on regular lattices and the quantities for percolation not on regular lattices, such as random lattices. In the present paper, we use Monte Carlo simulations to study bond percolation on $`L_1\times L_2`$ planar random lattices, duals of random lattices, and square lattices with free and periodic boundary conditions in vertical and horizontal directions, respectively, and with various aspect ratio $`L_1/L_2`$. We calculate the probability for the appearance of $`n`$ percolating clusters, $`W_n`$, the percolating probabilities, $`P`$, the average fraction of lattice bonds (sites) in percolating clusters, $`<c^b>_n`$ ($`<c^s>_n`$), and the probability distribution function for fraction $`c`$ of lattice bonds (sites), in percolating clusters of subgraphs with $`n`$ percolating clusters, $`f_n(c^b)`$ ($`f_n(c^s)`$). Using a small number of nonuniversal metric factors, we find that $`W_n`$, $`P`$, $`<c^b>_n`$ ($`<c^s>_n`$), and $`f_n(c^b)`$ ($`f_n(c^s)`$) for random lattices, duals of random lattices, and square lattices, have the same universal finite-size scaling functions. We also find that nonuniversal metric factors are independent of boundary conditions and aspect ratios. Furthermore, this study is related to recent developments in the universality and scaling of critical phenomena.

Universality and scaling are two important concepts in modern theory of critical phenomena stanley71 ; Fisher71 ; Privman84 , and percolation models are an ideal system for studying critical phenomena Stauffer94 . Thus, universality and scaling have been actively studied in recent decades, especially for percolation models Okabe00 . In 1992, Langlands et. al. Lang92  proposed that for bond and site percolation models on square (sq), planar triangular (pt), and honeycomb (hc) lattices, the critical existence probability (also called crossing probability or spanning probability) is a universal quantity, when aspect ratios of sq, pt, and hc lattices have relative ratios, $`1:\sqrt{3}:`$ $`\sqrt{3}/2`$. From 1995$``$1996, Hu, Lin and Chen (HLC) Hu95 ; hl96  calculated existence probability, $`E_p`$, percolation probability, $`P`$, and probability for the appearance of $`n`$ percolating clusters, $`W_n,`$ of bond and site percolation models, on sq, hc, and pt lattices with aspect ratios, $`1:\sqrt{3}:\sqrt{3}/2`$, and showed that all their scaled data fall on the same universal scaling functions, by selecting a very small numbers of nonuniversal metric factors, and maintaining similar nonuniversal metric factors under two boundary conditions, free and periodic boundary conditions. By using renormalization group theory, Hovi and Aharony in 1996 Hoviprl96  also pointed out that scaling functions, for the spanning probability are universal at the fixed point for every system with the same dimensionality, spanning rule, aspect ratio and boundary conditions. Okabe and Kikuchi in 1996 Okabeijmp96 , extended the work of HLC to a two-dimensional Ising model on planar regular lattices. In 1997, Hu and Wang hw97  found that the lattice percolation and continuum percolation of hard and soft disks, have the same universal scaling functions for $`W_n`$. Using the connection between an Ising model and a bond-correlated percolation model hu84 , Tomita, Okabe, and Hu in 1999 Tomita99 , calculated the probability for the appearance of $`n`$ percolating clusters, $`W_n,`$ the percolating probabilities, $`P`$, the average fraction of lattice sites in percolating clusters, $`<c>_n,`$ and the probability distribution function for the fraction $`c`$ of lattice sites in percolating clusters of subgraphs with $`n`$ percolating clusters, $`f_n(c)`$ , for bond-correlated percolation model on sq, hc, and pt lattices, with aspect ratios of $`1:\sqrt{3}:\sqrt{3}/2`$. Using a small number of nonuniversal metric factors, they found that $`W_n`$ , $`P`$, $`<c>_n,`$ and $`f_n(c)`$ for sq, hc, and pt lattices have the same universal finite-size scaling functions.

However, the studies mentioned above are mostly focused on regular lattices, with fixed coordination numbers hw97note . In 1999 Hsu and Huang (HH) Hsu1999  determined the percolation thresholds and critical exponents, and demonstrated explicitly that the ideas of universal critical exponents and universal scaling function with nonuniversal metric factors can be extended to bond percolation on $`L\times L`$ periodic planar random lattices, duals of random lattices, and square lattices, for both existence and percolating probabilities, and mean cluster size. This paper will study bond percolation on $`L_1\times L_2`$ planar random lattices, duals of random lattices, and square lattices, in more detail, and consider the case that the lattices have free and periodic boundary conditions, in vertical and horizontal directions, respectively, as in hl96 . Percolating probability is defined by the ratio of the number of bonds in the percolating clusters to the total number of bonds in Hsu1999 . Here, we consider two different definitions of the percolating probability, in terms of bonds and in terms of sites; the latter was also used in hl96  and Tomita99 . We calculate the probability $`W_n`$ for the appearance of $`n`$ percolating clusters, the percolating probability $`P`$, the average fraction of lattice bonds (sites) in percolating clusters, $`<c^b>_n`$ $`(<c^s>_n)`$, and the probability distribution function for the fraction $`c`$ of lattice bonds (sites) in percolating clusters of subgraphs with $`n`$ percolating clusters, $`f_n(c^b)`$ $`(f_n(c^s))`$, for various values of aspect ratios $`L_1/L_2`$, and finally check the universal finite-size scaling behaviors for these quantities. In Hsu1999 , HH used two nonuniversal metric factors, $`D_2`$ and $`D_3`$, to fix universal finite-size scaling functions, for percolating probability in terms of bonds. In the present paper, we calculated two nonuniversal metric factors of percolating probability, in terms of sites and obtained previous known values of nonuniversal metric factors determined by HH, to check whether we could have universal scaling functions for $`W_n`$, $`P`$, $`<c^b>_n`$ $`(<c^s>_n)`$, and $`f_n(c^b)`$ $`(f_n(c^s))`$, of bond percolation on random lattices, duals of random lattices, and square lattices.

Dirichlet and Voronoi Itzykson92  first used the concept of random lattices in condensed matter theory and Christ, Friedberg, and Lee (CFL) Christ82a  used another type of random lattices to formulate quantum field theory. Here, we adopt the CFL algorithm and give a brief review of the construction of planar random lattices and their duals. First, we randomly generate $`N`$ sites in the $`L_1\times L_2`$ rectangular domain with periodic boundary conditions. Next, we choose three nearby sites arbitrarily, and draw a circle to go through the three sites. If there are no lattice sites inside the circle, the three sites are connected by links to form a triangle. A planar random lattice is constructed by repeating the process until all sites are connected by links. The whole rectangular domain is divided into $`2N`$ non-overlapping triangles, whose vertices are sites of the random lattice, and circle centers with triangles are the sites of dual lattice. Thus, there is a one to one correspondence between triangles and dual lattice sites. Because a link of the random lattice is shared by two triangles, the two corresponding dual lattice sites are connected by one dual link. There is a one to one correspondence between links and dual links. The whole rectangular domain is partitioned into $`N`$ non-overlapping planar convex polyhedra, which are formed by dual links and the vertices of $`N`$ polyhedra are sites for dual lattice. There is also a one to one correspondence between the lattice sites and polyhedra on dual lattice. Examples of a planar random lattice with its dual, under periodic boundary conditions, in both vertical and horizontal directions, are shown in Fig. 1(a).

This paper is organized as follows: In Sec. II, we present the simulational results for $`W_n`$, $`P`$, $`<c^b>_n`$ ($`<c^s>_n`$), and $`f_n(c^b)`$ ($`f_n(c^s)`$) for bond percolation, on $`L_1\times L_2`$ random lattices, duals of random lattices, and square lattices, under free and periodic boundary conditions in vertical and horizontal directions with $`L_1/L_2=4`$. The boundary bonds which cross the rectangular domain in vertical direction on random lattices, due to periodic boundary conditions, are eliminated because of free boundary conditions in vertical direction considered in this paper. We adopt the method of HH Hsu1999 , to find percolating clusters. Only the first kind of percolating cluster paths, without boundary bonds in the vertical direction (the clusters extend from top to bottom), should be identified, and an example of this is shown in Fig. 1(b). In Sec. III, we use finite-size scaling theory to check the scaling behaviors of various quantities, and to show that such quantities have universal finite-size scaling functions for regular lattices and random lattices. A summary is provided in Sec. IV.

## II $`W_n(L_1,L_2,p)`$, $`f_n(c)`$, and $`<c>_n`$

We see the bond percolation on a lattice $`G`$, with linear dimensions $`L_1`$ and $`L_2`$ in horizontal and vertical directions, respectively; the probability for the appearance of $`n`$ top-to-bottom percolating clusters, $`W_n(L_1,L_2,p),`$ is defined by hl96

$$W_n(L_1,L_2,p)=\underset{G_n^{}G}{}p^{b(G_n^{})}(1p)^{Eb(G_n^{})}.$$
(1)
Here, the percolating cluster is defined as a cluster extending from top to bottom in $`G`$, $`G_n^{}`$ denotes a percolating subgraph with $`n`$ percolating clusters, $`b(G_n^{})`$ is the number of occupied bonds in $`G_n^{}`$, and $`E`$ is the total number of links in $`G`$. The existence probability $`E_p`$ can be expressed obviously as

$$E_p=\underset{n=1}{\overset{\mathrm{}}{}}W_n,$$
(2)
with $`W_0=1E_p`$.

To obtain more detailed information about the contents of the percolating cluster, following Tomita et. al. Tomita99 , we decompose $`W_n`$ as

$$W_n=_0^1f_n(c)𝑑c,$$
(3)
where $`n=1,\mathrm{}\mathrm{}`$, $`c`$ denotes the fraction of lattice bonds (sites) in percolating clusters, and $`f_n(c)`$ is the probability distribution function of $`c`$ in subgraphs with $`n`$ percolating clusters. The probability distribution function of $`c`$ in all subgraphs is the overall summations of $`f_n(c),`$ i.e.,

$$f(c)=\underset{n=1}{\overset{\mathrm{}}{}}f_n(c).$$
(4)
In terms of $`f_n(c)`$, the average fraction of lattice bonds (sites) in subgraphs with $`n`$ percolating clusters can be expressed as

$$<c>_n=_0^1cf_n(c)𝑑c,$$
(5)
where $`n=1,\mathrm{},\mathrm{},`$ and percolating probability $`P`$ can be written as

$$<c>=\underset{n=1}{\overset{\mathrm{}}{}}<c>_n=_0^{\mathrm{}}cf(c)𝑑c=P.$$
(6)
To generate subgraphs, we use the random bond occupation process with equal occupation probability for each link. The simulations are performed on $`128\times 32`$, $`256\times 64`$, and $`512\times 128`$ planar random (pran) lattices, and their duals (dpran), with free and periodic boundary conditions in the vertical and horizontal directions, respectively. To compare the results with regular lattices, we also perform simulations on square (sq) lattices of the same sizes. On each lattice, we take $`60`$ occupation probabilities around the critical percolation threshold for every $`0.002`$ increment, and use the random bond occupation process to generate $`10^510^6`$ configurations for each occupied probability, $`p`$. We calculate $`W_n(L_1,L_2,p)`$ , $`<c^b>_n`$, and $`<c^s>_n`$, where $`c^b`$ denotes the fraction of bonds in percolating clusters, and $`c^s`$ denotes the fraction of sites in percolating clusters; and the results are shown in Figs. 2 and 3. The calculated results of the percolating probabilities in terms of bonds, $`P^b,`$ and in terms of sites, $`P^s`$, are also shown in Fig. 3. We calculate $`f_n(c^b)`$ and $`f_n(c^s)`$ at $`p=p_c`$ and take $`p_c=0.3333`$ for planar random lattices, and $`p_c=0.6667`$ for dual lattice Hsu1999 . The results are shown in Fig. 4. The differences between bond and site contents in percolating clusters are shown in Figs. 3 and 4; here, for the clarity of presentation only, the results for $`512\times 128`$ lattices are plotted in the figures.

## III Universal Finite-size Scaling Functions

The finite-size scaling theory was first formulated by Fisher in 1971 Fisher71 . According to the theory, for a physical quantity $`X,`$ which scales as $`X(t)t^\rho `$ in a thermodynamic system near critical point $`t=0`$, then the same quantity in a finite system with linear dimension $`L`$, $`X_L(t)`$, should obey the general law,

$$X_L(t)L^{\rho /\nu }F(tL^{1/\nu }).$$
(7)
Here, $`F(x)`$ with $`x=tL^{1/\nu }`$ labelled as scaling function, with $`\nu `$ as correlation length exponent. In 1984, Privman and Fisher Privman84  considered universal finite-size scaling functions and nonuniversal metric factors, and proposed that the singular part of the free energy of a critical system can be written as

$$f_L^s(t)L^dY(DtL^{1/\nu }),$$
(8)
where $`d`$ is the spatial dimensionality of lattice, $`Y`$ is a universal scaling function, and $`D`$ is a nonuniversal metric factor.

At the critical point, $`p=p_c`$, there also exists a finite-size scaling form for the distribution function of $`X_L(t)`$ Tomita99 :

$$Q(X_L(t=0))L^{\rho /\nu }Y(X_L(t=0)L^{\rho /\nu }).$$
(9)
In Hu95 ; Hsu1999 , three nonuniversal metric factors $`D_1`$, $`D_2`$, and $`D_3`$ were used for regular lattices and random lattices, to describe the universal scaling functions of existence probability $`E_p`$, and the percolating probability $`P`$, i.e.

$$E_p(p,L)=F(x),$$
(10)
with $`x=D_1(pp_{c)}L^{1/\nu }`$ and

$$D_3P(p,L)=L^{\beta /\nu }S_p(z),$$
(11)
with $`z=D_2(pp_c)L^{1/\nu }.`$

Following Hu et. al. hl96 , we use the evaluated percolation thresholds $`p_c`$ Hsu1999 , and the exact values of critical exponent $`\nu =4/3`$, to plot $`W_n`$ as a function of $`y=(pp_c)L^{1/\nu }`$, for planar random lattices and their duals in Figs. 5(a) and 5(b), respectively. We can see from these results, that the scaled data for $`W_n`$ can be described by a single scaling function $`F_n(r,y)`$ with $`r=L_1/L_2`$, and $`F_n(r,y)`$ for $`n2`$ as a symmetric function of $`y`$. In Fig. 6 we plot $`W_n(L_1,L_2,p)`$ as a function of $`x`$ for bond percolation on $`512\times 128`$ random lattice, dual of random lattice, and square lattice, where $`x=D_1(pp_c)L^{1/\nu }`$, with $`D_1`$ taken from Table I and $`L=(L_1\times L_2)^{1/2}`$. Fig. 6 shows that the calculated results for each $`n`$ can be well described by a single universal scaling function, $`U_n(x)`$.

In Hu and Lin’s paper hl96 , the scaling functions $`F_n(r,y)`$ were calculated for bond percolation on square lattice for various values of aspect ratios, $`r`$. We will examine whether the same nonuniversal metric factors, $`D_1`$, can be extended to different aspect ratios. We calculate $`W_n`$ for $`L_1\times L_2`$ random lattices, dual of random lattices, and square lattices with $`r=L_1/L_2=1,`$ $`2,`$ $`\mathrm{},`$ $`6`$ and determine the universal scaling functions, $`U_n(r,x),`$ where $`x=D_1(pp_c)L^{1/\nu }`$ and where $`D_1`$ is taken from Table I. The results for $`n=1`$ and $`2`$ are shown in Figs. 7(a) and 7(b), respectively. We can see that the scaled data for each $`r`$ can be described by a single universal scaling function very well. The results of $`U_n(r,x)`$ as a function of $`r`$ for $`n=0,1,\mathrm{},`$ $`4`$ at the critical point $`p=p_c`$, are presented in Fig. 8(a) which show that the three lattices provide similar results. We also calculate the average number of percolating clusters $`C(r,x)`$, defined by

$$C(r,x)=\underset{n=1}{\overset{\mathrm{}}{}}U_n(r,x)n.$$
(12)
$`C(r,0)`$ for random lattice, dual of random lattice, and square lattice as a function of $`r`$ are shown in Fig. 8(b). Fig. 8(b) shows that $`C(r,0)`$ increases linearly with an increasingly large, $`r`$ and that different lattices have the same slope of approximately 0.43.

Tomita et. al. Tomita99  had obtained universal finite-size scaling functions of $`<c>_n`$, $`<c>,`$ $`f_n(c),`$ and $`f(c)`$ for a bond-correlated percolation model, corresponding to an Ising model on planar regular lattices. It is of interest to extend such study to bond random percolation on random lattices. From Eqs.(5), (6), and (11), the universal scaling function of $`<c>`$ and $`<c>_n`$ can be expressed as

$$D_3<c(p,L)>=L^{\beta /\nu }G(z),$$
(13)
and

$$D_3<c(p,L)>_n=L^{\beta /\nu }G_n(z),$$
(14)
with $`z=D_2(pp_c)L^{1/\nu }.`$ At $`p=p_c,`$ the universal scaling function of $`f(c)`$ and $`f_n(c)`$ are expressed as

$$D_3^1f(c)=L^{\beta /\nu }H(z^{}),$$
(15)
and

$$D_3^1f_n(c)=L^{\beta /\nu }H_n(z^{}),$$
(16)
with $`z^{}=D_3cL^{\beta /\nu }.`$ To check the finite-size scaling and universality of these quantities, we use simulation results for $`256\times 64`$ and $`512\times 128`$ square lattices, planar random lattices and their duals. In Hsu1999 , the percolating probability $`P`$ is defined in terms of the bond number in percolating clusters, and the nonuniversal metric factors $`D_2=D_2^b`$ and $`D_3=D_3^b`$ are used. To evaluate factors $`D_2^s`$ and $`D_3^s`$, we adopt the same procedure as in Hsu1999 , plotting $`P^b/L^{\beta /\nu }`$ and $`P^s/L^{\beta /\nu }`$ as functions of $`y=(pp_c)L^{1/\nu }`$, as shown in Fig. 9. All of the nonuniversal metric factors for different types of lattices used in this paper are listed in Table I.

We plot $`D_3P/L^{\beta /\nu }`$ and $`D_3<c>_n/L^{\beta /\nu }`$ as a function of $`z=D_2(pp_c)L^{1/\nu }`$, in Figs. 10(a) and 10(b) for bond content and site content, respectively, with $`D_2`$ and $`D_3`$ taken from Table I. At $`p=p_c,`$ the scaled data $`D_3^1f(c)/L^{\beta /\nu }`$, and $`D_3^1f_n(c)/L^{\beta /\nu }`$ as functions of $`z^{}=D_3cL^{\beta /\nu }`$, are presented in Figs. 11(a) and 11(b), respectively, for the bond content and site content. Figs. 10 and 11 show that the bond percolation processes on square lattices, random lattices and their duals have universal finite-size scaling functions.

## IV Summary and Discussion

Having used nonuniversal metric factors from Table I in this paper, we have found that the universal finite-size scaling functions for $`W_n`$ (Figs. 6 and 7), $`<c^b>_n`$ and $`P^b`$ (Fig. 10(a)), $`<c^s>_n`$ and $`P^s`$ (Fig. 10(b)), $`f_n(c^b)`$ (Fig. 11(a)), and $`f_n(c^s)`$ (Fig. 11(b)). Fig. 7 includes results for different aspect ratios $`r`$, i.e. $`6r1`$. The values of nonuniversal metric factors, $`D_1`$, $`D_2^b`$, and $`D_3^b`$ of Table I are consistent with the corresponding values of Ref. Hsu1999 , where the boundary conditions are different from the boundary conditions of the present paper. These results suggest that in random lattices, the nonuniversal metric factors are also independent of the boundary conditions and aspect ratios, as in the case of regular lattices Hu95 . Please also note that $`D_2^s`$ and $`D_3^s`$ in Table I are consistent within numerical errors.

Many interesting problems are related to the properties of multiple percolating clusters. It is of interest to extend the study of the present paper to higher spatial dimensions. In particular, the further study of multiple percolating clusters in three dimensions could be related to an oil drilling problem.

###### Acknowledgements.

This work was supported in part by the National Science Council of the Republic of China (Taiwan) under Contract No. NSC 89-2112-M-001-084.

Table I: The values of metric factors $`D_1`$, $`D_2^b`$, $`D_3^b`$, $`D_2^s`$ and $`D_3^s`$, for square lattices, random lattices and their duals, with free and periodic boundary conditions in vertical and horizontal directions, respectively.    Lattices   Square   Planar random   Dual of planar random    $`D_1`$   $`1`$   $`1.166\pm 0.020`$   $`1.177\pm 0.016`$  $`D_2^b`$   $`1`$   $`1.164\pm 0.014`$   $`1.176\pm 0.015`$  $`D_3^b`$   $`1`$   $`1.512\pm 0.008`$   $`0.778\pm 0.002`$  $`D_2^s`$   $`1`$   $`1.186\pm 0.012`$   $`1.180\pm 0.014`$  $`D_3^s`$   $`1`$   $`1.062\pm 0.001`$   $`1.005\pm 0.002`$

FIGURES

FIG. 1. Examples of (a) $`L_1\times L_2=8\times 4`$ planar random lattice (solid lines), with its dual (dashed lines) on $`L_1\times L_2=8\times 4`$ rectangular area, with periodic boundary conditions, and (b) a first kind of percolating cluster path, without boundary bonds (bold solid lines) on random lattice.

FIG. 2. $`W_n(L_1,L_2,p)`$ on square lattices, planar random lattices and their duals of size $`128\times 32`$, $`256\times 64`$, and $`512\times 128`$, which are represented by dotted, dashed, and solid lines, respectively.

FIG. 3. (a) $`P^b`$ and $`<c^b>_n,`$ and (b) $`P^s`$ and $`<c^s>_n`$ on square lattice, planar random lattice and its dual of same size $`512\times 128`$.

FIG. 4. At $`p=p_c`$, (a) $`f(c^b)`$ and $`f_n(c^b),`$ and (b) $`f(c^s)`$ and $`f_n(c^s)`$ on square lattice, planar random lattice and its dual of same size $`512\times 128`$.

FIG. 5. The scaled results of $`F_n(r,y)=`$ $`W_n(L_{1,}L_2,p)`$ as a function of $`y=(pp_c)L^{1/\nu }`$ for (a) planar random lattices and (b) their duals of size $`128\times 32`$, $`256\times 64`$, and $`512\times 128.`$ The monotonic decreasing function is for $`F_0(r,y).`$ The S shaped curve is for $`F_1(r,y)`$. The bell shaped curves from top to bottom are for $`F_n(r,y)`$, with $`n`$ as $`2`$, $`3`$, and $`4`$, respectively.

FIG. 6. The scaled results of $`U_n(x)=`$ $`W_n(L_{1,}L_2,p)`$, as a function of $`x=D_1(pp_c)L^{1/\nu }`$ for square lattice (solid curves), planar random lattice (dashed curves), and the dual of planar random lattice (dotted curves) of same size $`512\times 128.`$

FIG. 7. $`U_n(r,x)`$ for square lattice (solid curves), planar random lattice (dashed curves), and the dual of planar random lattice (dotted curves). (a) $`n=1,`$ the intersection of the curves on $`x=0`$ axis, from up to down are for $`r=L_1/L_2=1,2,3,\mathrm{},6.`$ (b) The width of curves, from small to large are for $`r=L_1/L_2=1,2,3,\mathrm{},6.`$

FIG. 8. (a) $`U_n(r,0)`$ as a function of $`r=L_1/L_2`$, for a number of percolating clusters (npc) run from $`0`$ to $`4`$, and (b) $`C(r,0)`$ as a function of $`r=L_1/L_2`$, with slope of the fitting line $`0.43`$. Square lattice $`(\mathrm{}),`$ planar random lattice $`(\mathrm{})`$ and the dual of planar random lattice $`(\times )`$, all have horizontal periodic boundary conditions.

FIG. 9. The scaled results of $`P(p,L)/L^{\beta /\nu }`$, in terms of bonds $`(P=P^b)`$ and sites $`(P=P^s)`$, for different types of lattice size $`256\times 64`$ and $`512\times 128`$, as a function of $`y=(pp_c)L^{1/\nu }`$.

FIG. 10. The scaled results of $`G_n(z,L)=`$ $`D_3<c(p,L)>_n/L^{\beta /\nu }`$, as a function of $`z=D_2(pp_c)L^{1/\nu }`$, for square lattices, planar random lattices and their duals of same size $`256\times 64`$, and $`512\times 128`$. (a) $`c=c^b`$, $`D_2=D_2^b`$, $`D_3=D_3^b`$, and (b) $`c=c^s`$, $`D_2=D_2^s`$, $`D_3=D_3^s`$.

FIG. 11. The scaled results of $`H_n(z^{},L)=D_3^1f_n(c)/L^{\beta /\nu }`$, as a function of $`z^{}=D_3(pp_c)L^{\beta /\nu }`$ for square lattices, planar random lattices, their duals of same size $`256\times 64`$, and $`512\times 128`$, with fitting results (solid curves); (a) $`c=c^b`$, $`D_3=D_3^b`$, and (b) $`c=c^s`$, $`D_3=D_3^s`$.


text:
# Non equilibrium effects in fragmentation

## I Introduction

The possibility of getting information about the thermodynamics of nuclear matter from the analysis of intermediate energies heavy ion collisions has triggered a lot of interest in this field. Starting from the pioneering work of the Purdue group, when for the first time a power law was used to fit the mass spectra resulting from highly excited nuclear systems and thus suggesting that a critical phenomena was taking place, it has open the study of a rather new branch of thermodynamics i.e the study of phase transitions (liquid-vapor) in finite systems.

One of the most challenging features that nuclear multi-fragmentation phenomena presents, is that signals of phase transitions of different orders can be extracted from experimental data. On one hand, the analogy between the nuclear force and a van der Waals interaction suggests that the nuclear equation of state should reproduce the main features that characterize a liquid-gas phase transition. Different indications of this kind of transition have been reported . On the other, signatures of a second order phase transition or critical behavior have been obtained assuming Fisher-like scaling relations for fragment distributions . Moreover, in recent works  is suggested that the observed critical behavior is compatible with a first order phase transition, and it is due exclusively to finite size effects.

An important point to be kept in mind is that the above mentioned approaches to the fragmentation problem are based on descriptions where not only the concept of thermodynamical equilibrium, but also the macroscopic constraints imposed to the system play a relevant role. For example, several statistical descriptions of the nuclear multi-fragmentation process, e.g. the statistical multifragmentation model (SMM)  and the microcanonical Metropolis Monte Carlo model (MMMC) , employ the concept of a freeze-out volume inside of which the existence of a thermodynamical equilibrated ensemble of fragments is assumed. These statistical models have been widespreadly used by the nuclear community in a successful manner to described some experimental observations. However, from the experience gained in numerical simulations, such a concept as freeze-out volume, even as an approximation, does not seem to be completely correct. It is then important to study which kind of differences arise as a consequence of not assumming a finite volume scenario, as it is not a priori evident that the evolution of a fragmenting system confined in a finite volume would produced the same macroscopic observable when compared with a non confined one, in which an expansive motion is present as an extra collective degree of freedom.

In previous works  we have studied the fragmentation of a simple classical system where the dynamics is governed by a Hamiltonian with a two body interaction Lennard-Jones term. A microscopic description, employing molecular dynamics techniques, was used in order to adequately handle the possible presence of a non equilibrium behavior. It was shown that a fragmentation time can be defined, after which a certain degree of local equilibrium is achieved in the system. This fact allowed us to calculate a caloric curve for our expanding-fragmenting system, which is characterized by the absence of a vapor branch.

The aim of the present communication is to study how the restriction of a finite volume, and then the imposition of equilibration, affects some of the results obtained in the unconstrained case. We will show that one of the main effects is seen in the behavior of the caloric curve (CC). For the constrained system it clearly shows a vapor-branch, which is absent in the free expanding case. This behavior is a direct consequence of the presence of a confining volume, which destroys the velocity correlations that, in the case of the free expanding system, build up an expansive flux that acts as a heat sink. Nevertheless for both cases (at rather low densities for the constrained system) a local maximum and a loop in the CC can be seen. This feature can be associated with a negative branch in the corresponding thermal response function (TRF), that might signal a phase transition of first order. A gradual smoothing of the mentioned indicators is observed for the confined system at higher densities. Differences in fragment mass distributions are also reported and related to the microscopic origin of such behavior.

This paper is organized as follows. In Section II we will describe the model used in our simulations. A brief review of the results already obtained for the unconstrained expanding system is included. Section III is devoted to the study of the caloric curves. In Section IV we study possible choices for the freeze out volume for the constrained case. In Section V we calculate the thermal response functions and the kinetic energy fluctuations in order to signal the presence of possible phase transitions. A microscopic correlations study is performed in Section VI analyzing the results of two different clusterization algorithms. Finally, in Section VII, conclusions are drawn.

## II Numerical simulations

The system under study is composed by excited drops made up of particles interacting via a 6-12 Lennard Jones potential, which reads:

$$V(r)=\{\begin{array}{cc}4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\left(\frac{\sigma }{r_c}\right)^{12}+\left(\frac{\sigma }{r_c}\right)^6\right]rr_c\hfill & \\ 0r\text{¿}r_c\hfill & \end{array}$$
(1)
We took the cut-off radius as $`r_c=3\sigma `$. Energies and distances are measured in units of the potential well ($`ϵ`$) and the distance at which the potential changes sign ($`\sigma `$), respectively while the unit of time used is: $`t_0=\sqrt{\sigma ^2m/48ϵ}`$. We integrated the set of classical equations of motion using the well known Verlet algorithm , taking $`t_{int}=0.002t_0`$ as the integration time step. Initial conditions were constructed using the already presented method of cutting spherical drops composed of 147 particles out of equilibrated and periodic 512 particles per cell L.J. system .

A broad energy range was considered such that the asymptotic mass spectra of the fragmented drops for the unconstrained system changes from a ”U shaped” pattern to an exponentially decaying one. Somewhere in between this two extremes a power law like spectra can be obtained.

### A Analysis of the unconstrained expanding case

For the sake of completeness we summarize in this section the main results obtained in previous work for the non-constrained fragmenting system (see  for details).

Due to the fact that we are dealing with a non-stationary process the determination of the time at which fragments are formed becomes one of the key ingredients in the analysis. To accomplish this task one can use a simple and intuitive cluster definition, that is based on correlations in configuration space: a particle $`i`$ belongs to a cluster $`C`$ if there is another particle $`j`$ that belongs to $`C`$ and $`|𝐫_𝐢𝐫_𝐣|r_{cl}`$, where $`r_{cl}`$ is a parameter called clusterization radius (In this work we used $`r_{cl}=r_{cut}=3\sigma `$). The recognition algorithm introduced by this definition is known as minimum spanning tree (MST) fragment recognition method.

In  it was shown that a different recognition method , the Early Cluster Formation Model and its companion practical realization the so called Early Cluster Recognition Algorithm (ECRA), outperforms the MST clusterization algorithm in the sense that it finds the asymptotic fragmentation pattern in phase space at early stages in the evolution, when fragments are still not observables in configuration space. At that time the system still looks like a rather compact piece of excited matter. Instead of being defined through a proximity criteria, the ECRA fragments are associated with the set of clusters $`\{C_i\}`$ for which the sum of the fragment internal energies attains its minimum value:

$`\{C_i\}`$ $`=\text{min}\left\{C_i\right\}[E_{\{C_i\}}=_iE_{int}^{C_i}]`$ (2)
$`E_{int}^{C_i}`$ $`=`$ $`{\displaystyle \underset{i}{}}[{\displaystyle \underset{jC_i}{}}K_j^{c.m.}+{\displaystyle \underset{j,kC_ijk}{}}V_{j,k}]`$ (3)
where the first sum in ( 3) is over the clusters of the partition, $`K_j^{c.m.}`$ is the kinetic energy of particle $`j`$ measured in the center of mass frame of the cluster which contains particle $`j`$, and $`V_{ij}`$ stands for the inter-particle potential. We dub the partition found by ECRA as the most bound density fluctuation in phase space (MBDF) and we define as time of fragment formation ($`\tau _{ff}`$) the time at which the MBDF attain microscopic stability (see appendix for details). In this way $`\tau _{ff}`$ is related to the time at which the systems switches from a regime dominated by fragmentation to one in which the dominant decay mode is evaporation of light aggregates (mostly single particles) by the excited fragments. In this way the ECFM-ECRA outperforms the MST not only in terms of the time at which fragments are detected but also in it’s capability of unveiling the nature of fragmentation process. It shows that fragments are formed in phase space as a consequence of correlations in both q-p space.

Once the time of fragment formation is determined it is possible to calculate several properties of the system at fragmentation time. For this purpose the expanding system is decomposed in concentric shells and the mean radial velocity is calculated. It can be seen that the expansion is almost linear with the distance to the center of mass of the system, and that a local temperature can be defined as the fluctuations of the velocity around the local expansive collective motion. Moreover, the isotropical character of those fluctuations supports the idea that local equilibrium is achieved .

It is worth noticing that, on one hand at $`\tau _{ff}`$ most of the system is still interacting, and on the other the local temperature of the inner shells attain a rather constant value that can be consistently considered as the temperature of the system at fragmentation time. This can be seen in figure 1 where we show, for different energies, the dependence of the local temperature, $`T_{loc}(r)`$, and the local density, $`\rho (r)`$, as a function of the distance to the center of mass at fragmentation time.

### B Analysis of the Constrained System

In order to study the consequences of imposing a finite volume constraint to our system we used a spherical confining ‘wall’. The considered external potential behaves like $`V_{wall}(rr_{wall})^{12}`$ with a cut off distance $`r_{cut}=1\sigma `$ , where it smoothly became zero along with its first derivative. A rather broad range of values for $`r_{wall}`$ were used.

Inside this potential, a highly excited drop was initialized in the way already described above and the corresponding equations of motion were integrated. Once the transient behavior was over we performed a microcanonical sampling of particle configurations every $`5t_0`$ up to a final time of $`140000t_0`$.

In this case the standard prescription for temperature calculation in the microcanonical ensemble was used, i.e. the kinetic energy ($`K`$) was related to the temperature of our $`N`$-particles system using:

$`T={\displaystyle \frac{2}{3(N1)}}K`$

## III Caloric curves

One of the thermodynamics measurements that remains useful in the small system limit (hundred of particles) is the caloric curve, i.e. the functional relationship of the system temperature with its excitation energy.

The resulting caloric curve (CC) for the unconstrained system is displayed in Fig. 2 a). In this figure circles denote the temperature of the system, already defined in Section II A, as a function of the energy. On the same figure a curve denoted by squares is also present, for this one the total kinetic energy, i.e. including both the fluctuations around the collective motion and the collective motion itself, is used instead. Of course this is not a temperature , just a fraction of the total kinetic energy, but it is then obvious that collective motion begins to be noticeable around a value of $`E=0`$ , and it becomes dominant at around $`E=2.ϵ`$. Two main features are to be noticed, in first place the CC develops a maximum, and second the CC has no ”vapor branch” but develops a rather constant behavior. These features are in contrast with the standard expectations inherited from the thermodynamics of infinite systems.

In Fig.3) we show the behavior of $`T(E,\rho )`$ for the constrained system (see also Fig.4a), where for the sake of clarity, we include just three different values of $`\rho `$).

For all density values a vapor branch is clearly observed. Moreover, comparing Fig.2a) with Fig.3) and Fig.4a), it is clear that the vapor branch is due to the presence of the confining volume which inhibits the formation of radial collective motion that behaves like a heat sink in the non constrained case. Inside the volume the emitted aggregates can do nothing but interact among themselves, and with the constraining wall, until thermal equilibrium is attained.

The effect of the boundaries on the CC is not important as long as the total energy in the system is small enough, and the expansive collective motion (in the unconstrained case) can be neglected. But this is not the case for highly excited finite systems, for which the presence of the constraining wall prohibits the formation of the expansive radial motion. For that range of energies the local equilibrium features that can be found in the fragmentation of the unconstrained system are replaced by global equilibrium ones when the dynamics is confined to finite volumes. This means that one has to be at least extremely careful if one intends to analyze expanding systems using confining volumes or any other model that rely on a global equilibrium hypothesis.

Another interesting feature to be notice in Fig.3) and Fig.4a) is the presence of a loop at the beginning of the gas rise in the CC for diluted enough systems. This back-bending behavior is gradually smeared out as the system density is increased. In Section V we will relate this to the behavior of the thermal response function of the system and in Section VI we will find that this can also be understood with the aid of the mass spectrum resulting from the MST and ECRA analysis.

## IV Freeze out volume approximation

As was already shown, due to the intrinsic non-equilibrium character of the fragmentation process in the non-confined case, an assumption of local equilibrium (‘mounted’ over an expansive radial flux) instead of a global one has been found to be more appropriate. Accordingly, distributions of density and temperature, and not unique values, are necessary to described the system at fragmentation time in a rigorous way(see Fig.1). It is then not expected to find an exact mapping in the ($`T,\rho `$) space between the constrained and non-constrained dynamics.

Nevertheless if one insists on establishing such a comparison, two different criteria can be adopted in order to choose the appropriate volume for the constrained case (which we will call freeze-out volume, $`V_{fo}`$). On one hand it can be seen from Fig.1) that at $`\tau _{ff}`$ the local density value for the inner cells (where most of the mass is present at the time of fragment formation) remains rather constant ($`\rho 0.08\sigma ^3`$) for the broad energy range presented in the figure. Accordingly, one can choose $`V_{fo}`$ in order to attain such density value in the constrained case. This election, that could be considered as a ‘density-guided’ choice, corresponds to $`r_{wall}8\sigma `$ and a freeze out density value of $`\rho _{fo}^+\rho _0/10`$ .

For the other hand a different approach can also be adopted looking at Fig.2a) and Fig.4a). It can be verified that for energy values where the onset of the fragmentation process occurs ($`2ϵE0ϵ`$) the thermal energy in the expanding system is clearly lower than for the $`r_{wall}=8\sigma `$ case in the constrained case. In fact, similar temperatures can be found in the constrained case only for much diluted situations, i.e. for $`r_{wall}=12\sigma `$. In this sense, from a ‘temperature point of view’, a different freeze out density $`\rho _{fo}^{}0.02\sigma ^3\rho _0/40`$ can be established.

As we mentioned above this ambiguity in the definition of a freeze out volume (or density) is a consequence of trying to reproduce the behavior of an expanding fragmenting system using a global equilibrium scenario.

## V Thermal response function

In recent works  a lot of attention has been paid to the role played by the behavior of the specific heat (or more generally speaking: the thermal response function, $`TRF`$) as a signal of the occurrence of a phase transition in finite systems. Moreover, it has been shown (see ) that the presence of a negative branch in the $`TRF`$ can be related to a first order phase transition taking place in an isolated finite system. This kind of analysis can straightforwardly be performed over the systems under the current study, using the already calculated CC and taking into account that the respective TRF’s can be calculated as:

$$TRF=(\frac{T}{E})^1$$
(4)
In Fig.2b) we show the TRF for the unconstrained system. Two poles can be seen. The first one ($`E0.5ϵ`$) signals the entrance of the system into the multifragmentation regime, while the second one is related to the leveling off in the corresponding CC, and can be related to the increasing limit imposed by the strong flux to the ’thermalization’ of the total available energy.

Fig.4b) shows the corresponding TRF for the constrained case. We include the curves for the two limiting cases discussed in Section IV as possible freeze out choices, i.e. $`\rho _{fo}^+`$ and $`\rho _{fo}^{}`$, and an intermediate $`\rho ^0\rho _0/20`$ value ($`\rho _{fo}^{}<\rho ^0<\rho _{fo}^+`$). We notice that two poles and a negative branch can be observed for the $`\rho _{fo}^{}`$ TRF curve. In this case we can relate the first pole with the onset of the transition, and the second one with the entrance in the gas phase, the energy distance between them being related to the associated latent heat. For $`\rho =\rho ^0`$ it can be seen that the curve exhibits a qualitatively similar behavior, the distance between the two poles is just reduced. Nevertheless, for higher density cases, such as $`\rho _{fo}^+`$, the two poles merge into a single ‘singularity’ limited by finite-size effects. In that case, the TRF remains positive for all energies, showing a peak as a signature of the transition, as a consequence of the fact that the corresponding caloric curve does not display a loop but instead a simple change of slope (see next section for an analysis of this change of behavior).

Before leaving this section we would like to discuss, using an approach introduced in , the behavior of the fluctuations in the kinetic energy as an indicator of the occurrence of a phase transition in our finite system. In  it is shown that for an isolated and equilibrated system in which the total energy can be decomposed as: $`E=E_1+E_2`$ the heat capacity can be calculated as:

$$C\frac{C_1^2}{C_1\sigma _1^2/T^2}$$
(5)
where $`C_1`$ is the heat capacity associated with the subsystem-$`1`$, $`\sigma _1`$ is the fluctuation of the partial energy $`E_1`$, and $`T`$ the temperature of the system. In this context, the presence of poles and negative values in TRF’s can be associated to abnormally large fluctuations of the partial energy stored in subsystem-$`1`$, i.e. $`\sigma _1^2C_1T^2`$, during the phase transition.

In Fig.5) we show, for the constrained case, the relative fluctuation, $`A(E)`$, for the kinetic energy defined as:

$$A(E)=N\frac{\sigma _K^2|_E}{<K>|_E^2}$$
(6)
where $`N`$ is the number of particles in the system, and $`<>`$ stands for an average over the microcanonical sampling of configurations. We show the direct calculations using equation 6 for both, $`\rho _{fo}^+`$ and $`\rho _{fo}^{}`$ , and we also include the estimation of $`A(E)`$ derived for the respective CC’s in Fig.5). A rather good agreement can be seen between both ways of calculating this magnitude. In the same graph a reference level marks the canonical value for $`A(E)`$. As expected, for the appropriate range of energies, ‘unusually’ large relative fluctuations for the $`\rho _{fo}^{}`$ case can be observed, while for the $`\rho _{fo}^+`$ case a local maximum, that does not exceed the canonical value can be seen. We postpone an analysis of this behavior until the next section, where we will study the properties of the system clusters.

We now try to extend the above performed analysis for the unconstrained system. In order to do that, we considered the central region (or core) defined by $`r<=6\sigma `$, and we calculated the total energy, $`E_{core}`$, and the number of particles $`N_{core}`$ that remains inside the core at fragmentation time. Binning the whole set of events according to these variables we can sample unconstrained fragmenting events in a pseudo-microcanonical way in what concerns core-quantities at $`\tau _{ff}`$. In Fig.6) we show the behavior of the relative fluctuation of the core kinetic energy, $`A(E_{core})`$, calculated using eq. 6 with $`K=K_{core}`$ and $`N=N_{core}=60\pm 5`$. In the same figure we include the canonical reference level. Even with a severe reduction of statistics imposed by the binning, a region of large fluctuations in $`K_{core}`$, for $`N_{core}`$ and $`E_{core}`$ fixed, can be easily recognized. This finding suggest that the ’large kinetic energy fluctuation signal’, counterpart of the behavior observed for the corresponding CC and TRF, can still be found for the expanding case.

## VI Cluster Distribution

In previous sections we have analyzed our system studying thermodynamical features like the behavior of the CC and the TRF. In this part of the work we will focus our attention to a different aspect of the process, more directly related to microscopic correlations, that can be studied using the clusterization algorithms already introduced in Section II A.

Let us begin with the simple MST analysis. As this cluster recognition method is exclusively based on spatial correlations one expects for high density situations a big MST-cluster to be present. On the other hand, for lower densities a deviation from the U-shaped behavior in the mass spectra can be expected. As was already stressed, MST underestimates the number of clusters because it does not take into account the relative velocities of the cluster constituent particles. For unconstrained systems this effect was apparent when we considered the time of fragment formation, at this time the system was already fragmented in phase space, but in configuration space, a big cluster was still present. Due to the fact that the system was free to expand the asymptotic stage was such that MST converged to ECRA. But in the constrained case such an expansion is prohibited and then MST will not converge to ECRA.

Nevertheless, an interesting point to be notice is that the MST algorithm can still provide useful information about the limit imposed by the constraining finite volume to the formation of well defined fragments in configurational space because it reflects the size of the interacting subsystems present at a given time.

To further explore this we studied the behavior of the MST clusterization when applied to constrained systems at high energies. We observed that once we enter the ‘vapor-branch’, the spectra remain unchanged. In Fig.7) we show the obtained MST cluster distribution for four different densities at energies that situate the systems in the vapor-branch in the respective CC. It can be seen that the MST-spectra go from a U-shape up to an exponential decaying behavior as the system density decreases. Notice that for densities around $`\rho \rho _{fo}^{}`$ a power-law like shape can be recognized. That means that from the spatial point of view the system admits, at least in principle, fragment configurations formed by well separated clusters of almost all possible sizes.

In Fig.8) we show the mean value of the size of the maximum MST-fragment as a function of the energy for the two densities taken as a reference. It can be seen that for the dense case, $`\rho =\rho _{fo}^+`$, the maximum fragment comprises almost all of the mass of the system irrespective of the value of the energy deposited in it. On the other hand for $`\rho =\rho _{fo}^{}`$ the size of the biggest fragment decreases with increasing energy. This results support the idea that for the high density constrained case, no surfaces can be built in order to get small MST-fragments at high densities, whereas at low densities MST-fragments are well defined structures in configuration space and the surfaces appear.

We will now discuss the results obtained within the ECFM model. In Fig.9) we show the ECRA and MST mass distribution for the unconstrained system ,both calculated at $`\tau _{ff}`$. We also show the same quantities calculated for the constrained case at densities $`\rho =\rho _{fo}^+`$ and $`\rho =\rho _{fo}^{}`$.

First we can notice that for the dense constrained case MST always gives essentially the same U shaped curve, while the ECRA spectra show the usual transition from U shaped to exponentially decaying behavior. On the other hand, for the low density constrained case and the unconstrained case both the MST and ECRA results show such a transition.

It is interesting to note that when comparing the most energetic case the ECRA mass distributions corresponding to the constrained cases present a steeper slope than the corresponding to unconstrained situation (figures 9 (d), (h), and (l)). This happens as a direct consequence of the existence of boundaries. For the constrained cases the violent inter particle collisions, that are present at all times due to the reflections of particles with the constraining walls, destroy almost every correlation in phase space. On the other hand, for the unconstrained case the system develops a collective expansive motion which limits the ‘thermal’ component of the kinetic energy, allowing for the build up of well correlated density fluctuations in phase space, i.e.ECRA-clusters (note also that while the ECRA clusters are microscopically stable for the unconstrained case, the ones for the constrained system are not). It is in this regime where the differences between a local equilibrium picture and a global one become important.

As long as the collective expansive energy in the free case can be neglected ($`E0.5ϵ`$) the ECRA spectra present a rather similar shape for the three shown cases. Moreover, a clear signature of a first order phase transition can be recognized in the shape of the ECRA spectra noticing the simultaneous presence of a ‘liquid-phase’ (mass fragments $`m30`$) and a ‘gas-phase’ (mass fragments $`m<10`$) for $`0.1ϵE0.4ϵ`$

This findings clarify the shapes of the corresponding CC. For the dense constrained system, $`\rho _{fo}^+`$, we have a big configurational MST-cluster all the time, inside of which an (unstable) ECRA partition can be defined that exhibits features expected for a first order phase transition. The fact that the whole system is strongly interacting (there are no spatial surfaces dividing ECRA-phases) is responsible of the reduction of the back-bending behavior in the corresponding CC and the lack of strong kinetic energy fluctuations. For the more dilute case, $`\rho _{fo}^{}`$, the system is capable of ‘reproduce’ phase space surfaces (at least at some degree) in configurational space and the mentioned signals start to be noticeable.

An extra support to the idea that fragmentation do take place in phase space not only for the free expanding system (see ) but also for the constrained case can be gained if we remove the external ‘wall’ potential and let the system expand. In Fig.10) we compare the asymptotic spectra resulting from such a process of removal of the constraints, with the ECRA one obtained with boundaries for $`\rho _{fo}^{}`$ and $`\rho _{fo}^+`$ cases at $`E=0.4ϵ`$ and $`E=0.7ϵ`$ . A quite remarkable agreement can be found. The fact that this agreement is better for lower densities reflects that different degree of spatial-fragmentation is achieved at different densities. In this way, for $`\rho _{fo}^+`$, the differences in both spectra (slight suppression of intermediate mass fragments favoring bigger clusters) can be understood as consequence of an aggregation-like process occurring during the early expansive stage due to the strong interacting nature of the system.

In Section V we showed for the confined system that the diverging behavior of the TRF is related to abnormally large fluctuations of the kinetic energy. In order to understand the origin of those fluctuations we show in Fig.11) the system kinetic energy $`K`$ and the mass of the biggest ECRA-cluster as a function of time for $`\rho =\rho _{fo}^{}`$ at $`E=0.4ϵ`$. Also the mass of the biggest MST-cluster is shown. It is clearly seen that the first two quantities are well correlated while the MST biggest cluster is not. This indicates that from a microscopic point of view the systems is sometimes mainly liquid (big ECRA biggest cluster) and others it is mainly vapor (small ECRA biggest cluster). This phase alternation is the kind of coexistence expected for finite system (see ).

## VII Conclusions

In this work we have presented a complete analysis of Caloric Curves, Thermal Response Functions and MST and ECRA fragment distributions for constrained and unconstrained excited liquid drops. We have shown that the effect of constraints is quite important, and that its main result is to allow the system to reach thermal equilibrium. In particular the Caloric Curves for constrained cases show a vapor branch which is clearly an ‘equilibrium effect’ that is not present for unconstrained systems.

We have also verified, studying the behavior of ECRA-cluster mass distributions, that fragmentation does occur in phase space in all the studied cases. Moreover, phase coexistence can be found in all ECRA spectra at appropriate energies (we have shown that the ECRA fragments in constrained systems correspond closely to the MST asymptotic clusters when the constraints are removed. This reinforces the validity of the ECFM-ECRA approach to study fragmentation phenomena occurring in finite volumes).

It is interesting to notice that for the constrained case, as the system density is increased, the expected signals of a first order phase transition in the CC, TRF and kinetic energy fluctuations are gradually smoothed, and eventually disappear. This can be associated to the fact that at low densities ”internal surfaces” can be developed in the constrained system allowing the transition to be traced in configurational space (i.e. MST clusters are formed reflecting the fact that well separated aggregates appear in coordinate space). Such a feature is not possible in the dense case, all of the time the system is composed by a big configurational cluster that comprises more than 95% of the total mass.

Moreover for the constrained case, a relation between the fluctuations in size of the maximum ECRA-fragment and the kinetic energy fluctuations was established. As a consequence, the presence of abnormally large fluctuations in the system kinetic energy could directly be linked with a phase coexistence phenomenon taking place in a finite system.

As a final remark we want to emphasize that, even some similarities can be found between the behavior of low density constrained systems and unconstrained ones at low energies, such a comparison fails at and above the energy that corresponds to the onset of the fragmentation process, i.e. when the collective radial modes begin to drive the evolution of the system.

As we mentioned, the role of this collective motion is to behave as a heat sink, precluding the system from developing a vapor branch, and freezing the most bound density fluctuations in phase space allowing them to become the asymptotic clusters. This is why a ‘local’ equilibrium picture, and not a ‘global’ one, is necessary to describe the fragmentation process correctly, i.e. taking into account the effects of the radial flux. This teaches us that when dealing with fragmentation phenomena of the kind appearing in nuclear multifragmentation, caloric curves which display vapor branches should be taken with caution, and consequently the results obtained from models that display such a feature should be at least critically reexamined.

### Acknowledgment

This work was done under partial financial support from the University of Buenos Aires via Grant No.TW98, and CONICET via Grant No. PIP 4436/96. A.Chernomoretz acknowledges CONICET for financial support and the warm hospitality of the Laboratoire de Physique Nucléaire, Département de Physique, Université Laval (Quebec, Canada). M.Ison aknowledges finantial support from UBA via a student scholarship.

### A Appendix

Once the clusters have been calculated using the ECFM we have to determine the time at which the fragmentation is over and the systems enters an evaporation regime (i.e the fragmentation pattern is formed and the fragments undergo a simple evaporation process). We will then define the time of fragment formation to that time at which the MBDF attain microscopic stability. This property of the clusters has been calculated using what we call the Short Time Persistence, in which we calculate the stability of MBDF against evaporation and coalescence. We define the time of stabilization in the following way: Given a configuration resulting from the ECRA analysis at a given time $`t`$ we analyze the microscopic stability of each fragment $`C_i^t`$ of size $`N_i^t`$ by searching on all the fragments $`C_j^{t+dt}`$ present at time $`t+dt`$ for the biggest subset $`N_{\mathrm{max}}]_i^{t+dt}`$ of particles that belonged to $`C_i^t.`$ We then we assign to this fragment a value $`STP_d=`$ $`\frac{N_{\mathrm{max}}]_i^{t+dt}}{N_i^t}.`$ In this way we are taking into account what we call the ”evaporation process”. We also have to take into account that there can be some realization for which cases the $`N_{\mathrm{max}}]_i^{t+dt}`$ does not constitute a subset of the original cluster $`C_i^t`$ but is embedded in a bigger fragment of mass $`N_i^{t+dt}`$. We include this effect by defining $`STP_i=\frac{N_{\mathrm{max}}]_i^{t+dt}}{N_i^{t+dt}}`$ (with $`i`$ standing for $`inverse)`$. Finally the Short time persistence reads :

$$STP(t,dt)=\left[\frac{STP_d(t,dt)+STP_i(t,dt)}{2}\right]_j_m_e$$
(7)
where $`\mathrm{}_m`$ is the mass weighted average over all the fragments with size $`N>3.`$ And $`\mathrm{}_e`$ is the average over an ensemble of fragmentation events at a given energy E.

A reference value for $`STP(t,dt)`$ can be obtained considering that the fragments undergo only a simple evaporative process. In this way we say that when $`STP`$ reaches this reference value the system goes from fragmentation to evaporation. We call this time $`\tau _{ff}`$.