text:
# An Analytic Model of Black Hole Evolution and Gamma Ray Bursts

## 1 INTRODUCTION

As is well known, the Blandford-Znajek(BZ) process was proposed originally as a possible energy mechanism of quasars and active galactic nuclei(AGNs;Blandford & Znajek 1977; Macdonald & Thorne 1982, hereafter MT82; Rees 1984). Recently, Lee, Wijers, & Brown (2000) proposed that the BZ process can be used as a central engine for powering gamma-ray bursts (GRBs), where the rotating energy of a stellar black hole (BH) with a magnetic field of $`10^{15}G`$ is extracted along the magnetic field lines supported by a transient magnetized accretion disk. Very recently, Lee and Kim (2000, 2002) proposed a model of the evolution of a rotating BH at the center of GRBs by considering the effects of the BZ process with the transient disk (hereafter the LK model), and some constraints to the parameters of the BH accretion disk are given.

In this paper an analytic model of BH evolution is proposed based on the LK model and is improved in two respects: (1) A parameter $`\lambda `$ is introduced for the position of the inner edge of the disk/torus, which is located between the innermost bound orbit and the last stable orbit. (2) Another parameter $`k`$ is used to indicate the ratio of the angular velocity of the magnetic field lines to that of the BH horizon. It is shown that the ratio $`k`$ being a little greater than 0.5 affects the evolutionary characteristics of the BH significantly, and the BH spin increases rather than decreases in its evolutionary process provided that the initial value of the BH spin is located in an appropriate value range determined by the ratio $`k`$. Our calculations show that the system of a BH accretion disk with $`k=0.6`$ might provide a much higher output energy in a shorter timescale for GRBs than the same system with $`k=0.5`$.

This paper is organized as follows: In §2 we derive the basic equations of BH evolution and the corresponding characteristic functions by introducing two parameters $`\lambda `$ and $`k`$. In §§3 and 4 we discuss the effects of the parameters $`k`$ and $`\lambda `$ on BH evolution and GRBs by using the parameter space and the corresponding characteristic functions, respectively. Finally, in §5 we discuss some problems concerning our model. Geometrized units ($`G=c=1`$) are used in this paper.

## 2 BASIC EQUATIONS OF BH EVOLUTION AND CHARACTERISTIC FUNCTIONS

Now we are going to describe our model for the evolution of a rotating BH surrounded by a magnetized accretion disk based on the LK model. First, a parameter $`\lambda `$ is introduced to indicate the position of the inner edge of the disk in the following equation:

$$\chi _{in}=\chi _{mb}+\lambda (\chi _{ms}\chi _{mb}),0\lambda 1,$$
(1)
where $`\chi _{in}`$ is the dimensionless radial parameter corresponding to the inner-edge radius $`r_{in}`$ of the disk, which is located between $`r_{mb}`$ and $`r_{ms}`$ (Abramowicz, Jaroszynski, & Sikora 1978; Abramowicz & Lasota 1980). The quantities $`\chi _{in},\chi _{mb}`$, and $`\chi _{ms}`$ are dimensionless radial parameters:

$$\chi _{in}\sqrt{\frac{r_{in}}{M}},\chi _{mb}\sqrt{\frac{r_{mb}}{M}},\chi _{ms}\sqrt{\frac{r_{ms}}{M}}.$$
(2)
It is found from equation (1) that the inner edge of a thick disk varies continuously from the innermost bound orbit to the last stable orbit as $`\lambda `$ changes from zero to unity, and the disk becomes a thin disk as $`\lambda `$ attains unity.

Usually the evolutionary state of a Kerr BH can be described by two parameters: BH mass $`M`$ and spin $`a_{}`$, and the latter is related to the BH mass $`M`$ and angular momentum $`J`$ by

$$a_{}\frac{J}{M^2},0a_{}<1.$$
(3)
Both $`\chi _{mb}`$ and $`\chi _{ms}`$ depend only on $`a_{}`$ (Novikov & Thorne 1973; Abramowicz et al. 1978):

$$\chi _{mb}=1+\sqrt{1a_{}},$$
(4)
$$\chi _{ms}^46\chi _{ms}^2+8a_{}\chi _{ms}3a_{}^2=0.$$
(5)
The specific energy $`E_{in}`$ and the specific angular momentum $`L_{in}`$ corresponding to the inner-edge radius $`r_{in}`$ of the disk/torus are expressed by (Novikov & Thorne 1973)

$$E_{in}=\frac{(12\chi _{in}^2+a_{}\chi _{in}^3)}{(13\chi _{in}^2+2a_{}\chi _{in}^3)^{1/2}},$$
(6)
$$L_{in}=\frac{M\chi _{in}(12a_{}\chi _{in}^3+a_{}^2\chi _{in}^4)}{(13\chi _{in}^2+2a_{}\chi _{in}^3)^{1/2}}.$$
(7)
In the LK model, the effects of the BZ process on the evolution of a rotating BH are taken into account in the environment of a transient magnetized accretion disk. Due to the lack of knowledge of the remote astrophysical load in the BZ process, it is doubtful that the matching condition $`k=0.5`$ can be satisfied exactly. Therefore, we introduce a parameter $`k`$ to describe the rotation of the magnetic field lines relative to the spinning BH,

$$k\frac{\mathrm{\Omega }_f}{\mathrm{\Omega }_h},0<k<1,$$
(8)
where $`\mathrm{\Omega }_f`$ is the angular velocity of the magnetic field lines and $`\mathrm{\Omega }_h`$ is that of the horizon and is related to the horizon radius $`r__h`$ by

$$\mathrm{\Omega }_h=\frac{a_{}}{2r__h},r__h=M(1+q).$$
(9)
By using a modified equivalent circuit for the BZ process based on MT82, we derive the expression for the rate of extracting energy from the rotating BH in the BZ process (hereafter the BZ power; Wang, Xiao, & Lei 2002) as follows:

$$P_{BZ}=4k(1k)P_{BZ}^{optimal},$$
(10)
$$P_{BZ}^{optimal}=B_h^2M^2[Q^1arctanQ\frac{1+q}{2}],$$
(11)
$$Q\sqrt{\frac{1q}{1+q}},q=\sqrt{1a_{}^2},$$
(12)
where $`P_{BZ}^{optimal}`$ is the optimal BZ power corresponding to the impedance matching with $`k=0.5`$ and $`B_h`$ is the magnetic field on the horizon. Both $`Q`$ and $`q`$ are the parameters depending only on the BH spin. Taking $`k=0.5`$ in equation (10), we have the same expression for the optimal BZ power as given by Lee et al. (2000), who pointed out that the BZ power had been underestimated by a factor of 10 in previous works.

Based on the conservation of energy and angular momentum, the basic evolutionary equations of a rotating BH with disk accretion by the BZ process can be written as (Park & Vishniac 1988; Moderski & Sikora 1996; Moderski, Sikora, & Lasota 1997; Wang, Lu, & Yang 1998; Lee & Kim 2000)

$$\frac{dM}{dt}=E_{in}\dot{M}_dP_{BZ},$$
(13)
$$\frac{dJ}{dt}=L_{in}\dot{M}_d\frac{P_{BZ}}{k\mathrm{\Omega }_h}.$$
(14)
Incorporating equations (13) and (14), we obtain the evolutionary equation of the BH spin:

$$\frac{da_{}}{dt}=M^2(L_{in}\dot{M}_d\frac{P_{BZ}}{k\mathrm{\Omega }_h})2M^1a_{}(E_{in}\dot{M}_dP_{BZ}).$$
(15)
The accretion rate $`\dot{M}_d`$ of the transient magnetized disk is related to $`B_h`$ by (Lee & Kim 2000)

$$\dot{M}_d=\frac{Mr_{in}^3}{\varpi ^2r_h^2}M^2a_{}^2B_h^2=\frac{B_h^2M^2Q^2\chi _{in}^2}{1+(a_{}/\chi _{in}^2)^2(1+2/\chi _{in}^2)},$$
(16)
where $`\varpi `$ is the Kerr metric parameter at the inner edge. According to the LK model, $`B_h`$ depends on the mass loss from the disk (Lee & Kim 2000):

$$B_h^2=B_h^2(0)D(t),$$
(17)
$$D(t)=1\frac{_0^t\dot{M}_d𝑑t}{M_d(0)},$$
(18)
where $`M_d(0)`$ and $`B_h(0)`$ are the initial disk mass and magnetic field, respectively.

Substituting the above expression for $`E_{in}`$, $`L_{in}`$, $`\chi _{in}`$, $`\mathrm{\Omega }_h`$, $`P_{BZ}`$ and $`\dot{M}_d`$ into equations (13) and (15), we obtain the following evolutionary equations:

$$\frac{dM}{dt}=f(a_{},k,\lambda )\dot{M}_d,$$
(19)
$$\frac{da_{}}{dt}=M^1g(a_{},k,\lambda )\dot{M}_d.$$
(20)
From equations (19) and (20), we find that $`dM/dt`$ and $`da_{}/dt`$ have the same sign as $`f(a_{},k,\lambda )`$ and $`g(a_{},k,\lambda )`$, and hereafter the two functions are referred to as the characteristic function of BH mass (CFBHM) and that of BH spin (CFBHS), respectively.

## 3 EFFECTS OF THE PARAMETER $`k`$ ON BH EVOLUTION AND GRBs

The characteristics of BH evolution can be well described by using the CFBHM and CFBHS, and the curves of these two functions varying with $`a_{}`$ for a thin disk ($`\lambda =1`$) with different values of $`k`$ are shown in Figure 1. It is shown in Figure 1a that $`f(a_{},k,1)`$ is negative when $`a_{}>0.9981,0.9984,`$ and $`0.9999`$ for $`k=0.5`$, 0.6, and 0.8, respectively. From Figure 1b we find that $`g(a_{},k,1)`$ is always negative for $`k=0.5`$, while it is positive for $`k=0.6`$ with $`0.3114<a_{}<0.8067`$, and for $`k=0.8`$ with $`0.1205<a_{}<0.9598`$. These results imply that the BH mass might decrease as the BH spin approaches unity very closely, and the BH could be spun up in the duration of powering the GRB, provided that the parameter $`k`$ is a little greater than 0.5.

The effect of the variation of $`k`$ on the evolution of the BH can be discussed more visually in the parameter space consisting of $`a_{}`$ and $`k`$ with the curves represented by $`f(a_{},k,1)=0`$ and $`g(a_{},k,1)=0`$, as shown in Figure 2.

It is found that the parameter space is divided into three regions by two thick solid curves, which are represented by $`g(a_{},k,1)=0`$ and $`f(a_{},k,1)=0`$. The former is the common boundary of regions I and II (curve abc), as shown in Figure 2a, and the latter is the common boundary of regions I and III, as shown in Figure 2b. Regions I and III are each further divided into two sub-regions, IA and IB, and IIIA and IIIB, by the horizontal thin line bd, as shown in Figures 2a and b, respectively. The segment bd is tangent to the curve $`g(a_{},k,1)=0`$ at the bottom point b(0.5610, 0.5118). In the parameter space, each filled circle with an arrow is referred to as a representative point (RP), which represents one evolutionary state of BH. We can use RP displacement in the parameter space to describe the evolution of the BH. According to the positions of RPs, we have five different BH evolutionary states corresponding to five subregions, as shown in Table 1.

From Table 1 we find that the five evolutionary states of a rotating BH can be easily determined by RP position in the parameter space, and the details are given as follows:

(1) The RPs in regions IA and IB represent the evolutionary states of the BH with increasing mass and decreasing spin, and the BH with the RP in region IA will never stop rotating, while the BH with the RP in region IB might evolve to a Schwarzschild BH with zero spin in long enough time. The evolution state given in the LK model is just described by the RP in region IB.

(2) The RPs in region II represent the evolutionary states of the BH towards the equilibrium spin with increasing mass and spin increasing toward equilibrium, which makes the BH a more powerful central engine for GRBs.

(3) The RPs in regions IIIA and IIIB represent the evolutionary states of the BH with decreasing mass and decreasing spin, which implies that the efficiency of transferring the accreted mass into GRB energy might be greater than unity in these two evolutionary states.

(4) As shown in Figure 2 and Table 1, there are two possible evolutionary terminals for the BH, if the evolution time is long enough. The terminal with equilibrium spin $`a_{}^{eq}`$ corresponds to the RPs in regions IA, II, and IIIA, while that with zero spin corresponds to the RPs in regions IB and IIIB.

From the above discussion, we find that evolutionary characteristics of a rotating BH in the BZ process are very sensitive to the parameter $`k`$, especially as $`k`$ varies across the vertical ordinate of the bottom point $`k^{bot}0.5118`$ in Figure 2a.

The value range of $`a_{}`$ in region II is expressed by $`a_{}^{turn}<a_{}<a_{}^{eq}`$, where $`a_{}^{turn}`$ and $`a_{}^{eq}`$ are, respectively, the left and right intersections of the curve abc with a horizontal line above the segment bd. Inspecting Figure 2a, we find that the more the horizontal line is above the segment bd or the greater $`k`$ is than $`k^{bot}`$, the greater the value range of $`a_{}`$ is in region II. The left part of curve abc is the set of $`a_{}^{turn}`$, while the right part of curve abc is the set of $`a_{}^{eq}`$. Obviously, $`a_{}^{turn}`$ is unstable, while $`a_{}^{eq}`$ is stable for the evolutionary states.

In order to discuss the effects of the parameter $`k`$ on the evolution of a rotating BH and the related characteristics of GRBs, we compare the results of the BH accretion disk system by taking $`k=0.6`$ and $`k=0.5`$, respectively. Our main results are shown in Figures 3-8, in which the solid lines are used for the results with $`k=0.6`$, and the dashed lines for those with $`k=0.5`$. In our calculations, the initial values of the concerned quantities are taken as those given in the LK model: $`M_h(0)=7M_{}`$, $`M_d(0)=3M_{}`$, and $`B_d(0)=10^{15}G`$.

(1) In Figure 3 the BH spin $`a_{}`$ increases with the evolution time if its initial value is taken to be $`0.3114<a_{}(0)<0.8067`$ for $`k=0.6`$, while it always decreases with the evolution time no matter what initial value it takes for $`k=0.5`$. As shown in Figure 1b, 0.3114 and 0.8067 are exactly the values of $`a_{}^{turn}`$ and $`a_{}^{eq}`$ for $`k=0.6`$, respectively.

(2) In Figure 4 $`a_{}^c`$ can never be greater than $`a_{}^{turn}=0.3114`$ no matter how much the initial disk mass is for $`k=0.6`$, while it increases monotonically with the initial disk mass for $`k=0.5`$. Since $`a_{}^c`$ is defined as the critical BH spin to ensure that the BH will stop spinning after all the initial disk mass accretes onto the BH in the LK model, our result implies that a BH with its initial spin greater than $`a_{}^{turn}`$ can never be spun down to a Schwarzschild BH.

(3) In Figure 5 the output energy $`E`$ from a system with $`k=0.6`$ (solid lines) is always greater than that from a system with $`k=0.5`$ (dotted lines). In the former case, $`E`$ will monotonically increase with the initial disk mass $`M_d(0)`$ without an upper limit, provided that the initial BH spin is greater than $`a_{}^{turn}`$, while in the latter case, it always keeps constant for $`M_d(0)`$ greater than one critical value. Another difference between these two systems lies in the derivative of $`E`$ with respect to $`M_d(0)`$. As shown in Table 2 the first derivative $`dE/dM_d(0)`$ always increases with $`M_d(0)`$ (in units of solar mass) provided that the system has the parameter $`k=0.6`$ with $`0.3114<a_{}(0)<0.8067`$, while it always decreases with $`M_d(0)`$ either for $`k=0.6`$ with $`0<a_{}(0)<0.3114`$ or for $`k=0.5`$ with any value of $`a_{}(0)`$. These results can be easily explained by the RPs in the parameter space (Figure 2a) and the relation of the BZ power to the BH spin.

(4) In Figure 6 the output energy from the system with $`k=0.6`$ is always greater than that with $`k=0.5`$ for any initial values of the BH spin. Defining the relative difference of output energy between the above two systems as $`\delta (E_{out}^{0.6}E_{out}^{0.5})/E_{out}^{0.5}`$, we have the curve of $`\delta `$ versus $`a_{}(0)`$ as shown in Figure 7, and the maximum is $`\delta _{max}2.44`$ at $`a_{}(0)0.294`$. Thus, the output energy for GRBs can be remarkably augmented by a little increase above $`k^{bot}`$.

(5) In Figure 8 the evolutionary timescales $`t_{evl}`$ of GRBs versus $`a_{}(0)`$ are obtained by incorporating equations (16) $``$ (18) and the initial values of the systems with $`k=0.5`$ and $`k=0.6`$, where the cutoff of $`t_{evl}`$ is taken as $`T_{90}`$ defined in the LK model. A peak value of $`t_{evl}`$ is found for each system, which is located at $`a_{}(0)0.298`$ for $`k=0.6`$ and at $`a_{}(0)0.402`$ for $`k=0.5`$, respectively. It is shown that the $`t_{evl}`$ of the system with $`k=0.6`$ is always less than that with $`k=0.5`$ for $`a_{}(0)>0.366`$. These results might be helpful to explain GRBs with the most energetic power.

## 4 EFFECTS OF THE PARAMETER $`\lambda `$ ON BH EVOLUTION AND GRBs

Taking the parameter $`\lambda `$ into account, we can discuss the evolutionary characteristics of a BH surrounded by a thick disk/torus, and the parameter spaces with boundary curves $`g(a_{},k,\lambda )=0`$ and $`f(a_{},k,\lambda )=0`$ corresponding to different values of $`\lambda `$ are shown in Figure 9. It is found that the above five evolutionary states of a BH surrounded by a thin disk $`(\lambda =1)`$ remain valid in the case of a thick disk/torus with its inner edge positioned by the parameter $`\lambda `$ $`(0\lambda <1)`$. It is interesting to note that both the bottom point of the curve $`g(a_{},k,\lambda )=0`$ and the leftmost point of the curve $`f(a_{},k,\lambda )=0`$ vary nonmonotonically with the parameter $`\lambda `$, and we have the curves of $`k^{bot}`$ and $`a_{}^{left}`$ versus $`\lambda `$ as shown in Figures 10 and 11, respectively. The quantity $`k^{bot}`$ attains its maximum of 0.598 at $`\lambda 0.137`$, and $`a_{}^{left}`$ attains its minimum of 0.995 at $`\lambda 0.261`$. The output energy $`E_{out}`$ from the systems and the evolutionary timescales $`t_{evl}`$ corresponding to the different values of $`\lambda `$ with $`k=0.6`$ are compared in Figures12 and 13, respectively. Finally, the relative differences of output energy defined by $`\delta (E_{out}^{0.5}E_{out}^{0.5})/E_{out}^{0.5}`$ corresponding to the different values of $`\lambda `$ are shown in Figure 14. The main effects of the parameter $`\lambda `$ on BH evolution and GRBs are summarized as follows.

(1) The BH will be spun up in the case of a thick disk $`(0\lambda <1)`$ as well as in the case of a thin disk $`(\lambda =1)`$, provided that we have $`k>k^{bot}`$ with $`a_{}^{turn}<a_{}(0)<a_{}^{eq}`$, where $`a_{}^{turn}`$ and $`a_{}^{eq}`$ have the same meaning as given in §3.

(2) The efficiency of transferring the accreted mass into GRB energy might be greater than unity provided that the RP of BH evolutionary state is located in the corresponding region III as shown in Figure 9b.

(3) From Figures 12 and 13 we find that the output energy from a system with a thick disk/torus might be greater, and the concerned evolutionary timescales might be longer, than the corresponding values of a system with a thin disk if $`a_{}(0)`$ is greater than some value. Therefore, our model with different values of $`\lambda `$ might adapt to the GRBs with different output energy and different timescales.

(4) From Figure 14 we find that the output energy from a system with $`k=0.6`$ is always greater than that from a system with $`k=0.5`$, no matter what value $`a_{}(0)`$ has. This is true for a system with a thick disk/torus as well as for a system with a thin disk.

In summary, compared to the effects of $`k`$ the characteristics of BH evolution and GRBs are not remarkably affected by variation of the parameter $`\lambda `$.

## 5 DISCUSSION

In this paper the evolutionary characteristics of a rotating BH surrounded by a transient magnetized disk are discussed by considering the coexistence of disk accretion with the BZ process, and the two parameters $`\lambda `$ and $`k`$ are introduced in our model to modify the LK model. A main consequence of our model is that the BH will be spun up rather than spun down provided that $`k`$ is greater than a critical value $`k^{bot}`$ with $`a_{}^{turn}<a_{}(0)<a_{}^{eq}`$. Our calculations show that more energy from a BH system with $`k`$ above $`k^{bot}`$ might be extracted in less time compared to a BH system with $`k`$ below $`k^{bot}`$.

However there are still some problems related to our model. First, the different values of the ratio $`k\mathrm{\Omega }_f/\mathrm{\Omega }_h`$ play a very important role in our model. Unfortunately, the determination of this ratio has remained one of the main unresolved problems in the BZ process. It was argued in a speculative way that the ratio might be regulated to about 0.5, if the charged particles conspired with the BH to have just the right amount of inertia for the impedance matching (MT82). However, it is difficult to understand how the load can conspire with the BH to have the same resistance and satisfy the matching condition, since the load is so far from the BH that it cannot be causally connected (Punsly & Coroniti 1990).

It was pointed out in the theory of BH magnetospheres that $`\mathrm{\Omega }_f`$ can be regarded as a function of magnetic flux $`\mathrm{\Psi }`$, i.e., $`\mathrm{\Omega }_f(\mathrm{\Psi })`$, since it is constant on the magnetic surface due to isorotation and axial symmetry (MT82). Very recently, Beskin and Kuznetsova (2000) discussed the stream equation describing magnetic surfaces $`\mathrm{\Psi }(r,\theta )`$, and found that $`\mathrm{\Omega }_f`$ can be either greater or less than $`\mathrm{\Omega }_h/2`$ for the different possibilities of the flow, with the electric current fixed by the pair creation region in the BH magnetosphere. However, we still have a long way to go for the determination of $`\mathrm{\Omega }_f`$ in a BH magnetosphere.

Second, we notice that the relation in equation(16) between $`\dot{M}_d`$ and $`B_h`$ is the key to affecting the evolutionary characteristics. To illustrate this, we replace equation (16) with another relation between $`\dot{M}_d`$ and $`B_h`$, which was proposed based on the balance between the pressure of the magnetic field on the horizon and the ram pressure of the innermost parts of the accretion flow (Moderskin et al. 1997), i.e.,

$$\frac{B_h^2}{8\pi }=P_{ram}\rho c^2\frac{\dot{M}_d}{4\pi r_h^2}.$$
(21)
From equation (21) we assume the relation to be

$$\dot{M}_d=\frac{r_h^2B_h^2}{2}.$$
(22)
Hereafter equation (22) is referred to as the MSL relation. Substituting the MSL relation into equation (13) and (15), we obtain the following evolutionary equations:

$$\frac{dM}{dt}=f_{_{MSL}}(a_{},k,\lambda )\dot{M}_d,$$
(23)
$$\frac{da_{}}{dt}=M^1g_{_{MSL}}(a_{},k,\lambda )\dot{M}_d,$$
(24)
where $`f_{_{MSL}}(a_{},k,\lambda )`$ and $`g_{_{MSL}}(a_{},k,\lambda )`$ are CFBHM and CFBHS corresponding to the MSL relation, respectively. In the same way, we obtain the parameter space consisting of $`a_{}`$ and $`k`$ as shown in Figure 15. The parameter space is very different from that depicted in Figure 2, and it is only divided into two regions by the boundary curve $`g_{_{MSL}}(a_{},k,\lambda )=0`$, corresponding to only two states of BH evolution. The characteristics of BH evolution related to the MSL relation are shown in Tabel 3.

In this paper an analytic model of BH evolution is proposed for GRBs. Although this model is rather particular and speculative , the potential application is still attractive: a little variation of the angular velocity of the magnetic field lines might affect the output energy and evolutionary timescales of a BH system remarkably. We also hope to improve and extend this model to AGNs and stellar-mass BHs in future work, since there is more evidence to strongly suggest the existence of fast-spinning BHs and possibly energy extraction via Blandford-Znajek-like processes(Krolik 2001).

The anonymous referees are thanked for their suggestions on the improvement of our manuscript. V. S. Beskin is thanked for his encouragement and the discussion of the determination of $`\mathrm{\Omega }_f`$ in a BH magnetosphere. This work is supported by the National Natural Science Foundation of China under grant No. 10173004.


format: text

text:
# Free energy determination of phase coexistence in model C60: A comprehensive Monte Carlo study

## I Introduction

We report an extensive investigation of the free energy characteristics of the Girifalco model of C60 fullerene . As is well known, this representation hinges on the fact that C60 molecules have almost spherical shape and freely rotate at sufficiently high temperatures . Under these conditions, the hollow molecular cages can be assimilated to spheres whose surface consists of a uniform distribution of carbon sites. The overall interaction between two fullerene particles is then obtained by an integral of the interaction between pairs of sites on different cages, eventually yielding an analytic central two-body potential. The latter is characterized by a harshly repulsive core at short distance, followed by a deep attractive well which rapidly decays with the interparticle distance .

The Girifalco model constitutes a prototype system in several respects; recent studies suggest that a similar “smeared out” spherical description can be attempted for Cn\>60 systems, although fullerene molecules with $`n>60`$ can have a sensibly non-spherical shape; the cases $`n=70`$, 76 and 84 have in particular been examined (see  and references therein). On the other hand, more refined calculations of the fullerene-fullerene interaction yield results very close to those predicted through the Girifalco model , and a similar representation has been used for the description of other hollow nanoparticles as carbon onions, or metal dichalcogenides (also termed inorganic fullerenes) as GaAs and CdSe . Moreover, a modification of the Girifalco model, suitable for the description of solid C60 at low temperatures, has been recently proposed ; this development seems of particular interest since fullerites, doped with organic molecules and upon the injection of electron (or holes), exhibit a superconducting behavior up to $`T=112`$ K ; an accurate description of such a simple model might prove useful for further studies on the lattice behavior upon impurity doping.

In our opinion, such a possible reference role of the Girifalco model for further studies on fullerenes and other systems, calls for a complete and confident determination of its phase diagram. With this purpose in mind, we have investigated the free energy characteristics of the C60 model for both the solid and the fluid phase through extensive Monte Carlo simulations, spanning the whole high-temperature region of the phase diagram.

A second general aspect of the Girifalco potential is related to its short-range nature. It is useful to recall in this respect that a ubiquitous definition of the range of interaction has been recently proposed in Ref. . As several studies have pointed out, the most apparent consequence of a reduced interaction length is the metastability of the liquid-vapor equilibrium, which is preempted by the fluid-solid coexistence . It is actually known that a stable liquid-vapor coexistence still survives for the model envisaged here, albeit restricted to a few tens degree temperature range . A tiny liquid pocket has been predicted also in the phase diagram of similar models of higher order fullerenes , although it is argued that even a modest reduction of the range of the forces might cause the disappearance of the liquid phase. The emerging borderline nature of the Girifalco model implies that the overall appearance of its phase diagram sensitively depends on the details of the interaction potential; in fact, the initial controversy around the existence of a stable liquid phase for this system  has been solved by taking into account on one side the full role of the attractive part of the interaction , and on the other, a conveniently large simulation sample . More generally, it has been recognized in Ref.  that the physical behavior of a wide class of systems characterized by short-range interactions, and in particular the onset of freezing, are substantially affected by the “perturbative” part (with respect to the repulsive core) of the potential, rather than being dominated by excluded-volume and packing (i.e. by entropic) effects, as is the case in the currently accepted van der Waals picture of simple liquids . We here investigate, as a further key purpose of this work, the same issue for the C60 model and try to reconcile the predictions coming from different structural indicators — usually related to the freezing threshold of several simple liquids  and adopted in early calculations for this model  — with the results of full free energy calculations.

The paper is organized as follows: we present in Section II the model and the simulation strategies; the results are reported and discussed in Sect. III, while Sect. IV is devoted to the conclusions and future perspectives of our investigation. A preliminary account of this work has been presented elsewhere .

## II Model and simulation strategies

The Girifalco potential between two C60 molecules, is written as :

$`v(r)`$ $`=`$ $`\alpha _1\left[{\displaystyle \frac{1}{s(s1)^3}}+{\displaystyle \frac{1}{s(s+1)^3}}{\displaystyle \frac{2}{s^4}}\right]`$ (2)
$`+\alpha _2\left[{\displaystyle \frac{1}{s(s1)^9}}+{\displaystyle \frac{1}{s(s+1)^9}}{\displaystyle \frac{2}{s^{10}}}\right],`$
where $`s=r/d`$, $`\alpha _1=N^2A/12d^6`$, and $`\alpha _2=N^2B/90d^{12}`$; $`N=60`$ and $`d=0.71`$ nm are the number of carbon atoms and the diameter, respectively, of the spherical particles; $`A=32\times 10^{60}`$ erg cm6 and $`B=55.77\times 10^{105}`$ erg cm12 are constants entering the Lennard-Jones 12-6 potential through which two carbon sites on different molecules are assumed to interact. The finite distance at which the potential (2) crosses zero and the minimum of the potential well depth are $`\sigma 0.959`$ nm and $`\epsilon 0.444\times 10^{12}`$ erg at $`r_{\mathrm{min}}=1.005`$ nm, respectively.

In order to investigate the coexistence properties of model (2), it is required the knowledge of the free energies of both the solid and the fluid phase. As a general strategy, the free energy of the system is first evaluated all along a supercritical isotherm at temperature $`\overline{T}`$ by integrating the pressure $`P`$ as a function of the density $`\rho `$ according to the formula :

$$\frac{\beta F(\rho ,\overline{T})}{N}=\frac{\beta F(\overline{\rho },\overline{T})}{N}+_{\overline{\rho }}^\rho \frac{\beta P(\rho ^{})}{\rho ^{}}\frac{d\rho ^{}}{\rho ^{}};$$
(3)
here $`F/N`$ is the Helmholtz free energy per particle, $`\beta =1/k_BT`$ is the inverse temperature, $`k_B`$ is the Boltzmann constant and $`(\overline{\rho },\overline{T})`$ is a thermodynamic state where the free energy is known (see below). The free energy at different temperatures is then calculated along isochoric paths as:

$$\frac{\beta F(\rho ,T)}{N}=\frac{\beta F(\rho ,\overline{T})}{N}_{\overline{T}}^T\frac{U(T^{})}{Nk_BT^{}}\frac{dT^{}}{T^{}},$$
(4)
where $`U/N`$ is the internal energy per particle of the system. The equilibrium conditions finally derive from the equality of the chemical potential and the pressure of the different phases.

All quantities entering Eqs. (3) and (4) have been determined through standard Monte Carlo simulations at constant volume or pressure. Simulations have been mostly carried out on a sample composed of $`N=864`$ C60 particles enclosed in a cubic box with periodic boundary conditions. The C60 interaction has been considered up to half the box length.

As far as the determination of the reference free energy $`F(\overline{\rho },\overline{T})`$ in Eq. (3) is concerned, we have used for the solid phase the Einstein crystal method described by Frenkel and Ladd . Namely, the C60 interaction in the solid fullerite is smoothly transformed, through a coupling parameter $`\lambda `$, into a corresponding harmonic potential so that the configurational energy $`U(\{𝐫\})`$ takes the following form :

$`U(\{𝐫\})`$ $`=`$ $`U_{\mathrm{C}_{60}}(\{𝐫_0\})+(1\lambda )\left[U_{\mathrm{C}_{60}}(\{𝐫\})U_{\mathrm{C}_{60}}(\{𝐫_0\})\right]`$ (6)
$`+\lambda {\displaystyle \underset{i=1}{\overset{N}{}}}\alpha (𝐫_i𝐫_{0,i})^2,`$
where $`𝐫_{0,i}`$ is the lattice position of atom $`i`$ and $`U_{\mathrm{C}_{60}}(\{𝐫_0\})`$ is the static contribution to the potential energy; $`\alpha `$ is the spring constant of the Einstein crystal. The free energy difference can then be written as:

$`{\displaystyle \frac{\beta F(\overline{\rho },\overline{T})}{N}}`$ $``$ $`{\displaystyle \frac{\beta F_{\mathrm{C}_{60}}}{N}}`$ (7)
$`=`$ $`{\displaystyle \frac{\beta F_{\mathrm{ein}}}{N}}+{\displaystyle _1^0}d\lambda {\displaystyle \frac{\beta }{N}}{\displaystyle \frac{U(\lambda )}{\lambda }}_\lambda ,`$ (8)
where the configurational free energy of the Einstein crystal is:

$`{\displaystyle \frac{\beta F_{\mathrm{ein}}}{N}}`$ $`=`$ $`{\displaystyle \frac{\beta U_{\mathrm{C}_{60}}(\{𝐫_0\})}{N}}{\displaystyle \frac{3(N1)}{2N}}\mathrm{ln}\left({\displaystyle \frac{\pi }{\alpha \beta }}\right)`$ (10)
$`+{\displaystyle \frac{3}{2N}}\mathrm{ln}N{\displaystyle \frac{\mathrm{ln}V}{N}}.`$
In the equation above, the last two terms account for the fixed center of mass constraint at which simulations are performed .

As for the absolute free energy of the fluid phase, the chemical potential $`\mu `$ has been estimated at several intermediate densities through the Widom test particle method .

In parallel with the free energy, the entropy of the fluid phase is also systematically analyzed, given that in several earlier papers the onset of freezing of the Girifalco model has been associated with the vanishing of the residual multiparticle entropy $`\mathrm{\Delta }s`$, namely

$$\mathrm{\Delta }ss_{\mathrm{ex}}s_2=0,$$
(11)
according to the one-phase freezing criterion originally proposed in Ref. . In Eq. (11), $`s_{\mathrm{ex}}`$ is the the excess entropy per particle of the system (in $`k_B`$ units) and $`s_2`$ is defined in terms of the radial distribution function $`g(r)`$ of the system as :

$$s_2=\frac{\rho }{2}\{g(r)\mathrm{ln}[g(r)]g(r)+1\}d𝐫.$$
(12)
In previous works $`\mathrm{\Delta }s`$ has been evaluated using liquid state integral equation theories. Here we report Monte Carlo data for the structural and thermodynamics quantities entering Eqs. (11) and (12), in order to perform a rigorous test of the theoretical predictions based on the criterion (11).

## III Results and discussion

The solid and liquid branches of the equation of state, to be integrated in Eq. (3), have been calculated at the supercritical isotherm $`\overline{T}=2100`$ K through both NVT and NPT Monte Carlo simulations, in order to crosscheck the predictions of the two algorithms. Five to six runs of 25 000 steps have been performed at each thermodynamic point investigated. The compressibility factor $`\beta P/\rho `$ is shown in Fig 1.

Results for the absolute free energy of the solid phase at $`\overline{T}=2100`$ K and $`\overline{\rho }=1.375`$ nm-3, as obtained through the Einstein crystal method, are reported in Fig. 2. Constant volume simulations have been carried out for several values of the switching parameter $`\lambda `$. The spring constant is set to a value $`\alpha /\epsilon =490`$, which makes the interactions in the pure Einstein crystal as close as possible to those of the original system, so to optmize the accuracy of the numerical integration scheme of Eq. (7. We have analyzed the free energy dependence on the system size; results with 256 and 2916 particles are shown in Fig. 2. It appears that the effect of $`N`$ on the estimate of the free energy is small but systematic; in the thermodynamic limit ($`N\mathrm{}`$) the smooth extrapolation in the bottom panel of Fig. 2 yields $`\beta F/N(\overline{\rho }=1.375\mathrm{nm}^3,\overline{T}=2100\mathrm{K})1.401`$.

As for the fluid phase, the excess chemical potential at $`\overline{T}=2100`$ K and $`\overline{\rho }=0.60`$ nm-3 has been calculated through the Widom test particle method as $`\beta \mu _{\mathrm{ex}}=1.523`$. Several tests at higher densities have also been conducted, in order to assess the results obtained via Eq. (3). The free energy at $`T=2100`$ K is shown in the top panel of Fig. 3, along with the common tangent construction, which determines the coexistence conditions. The $`\mu `$ vs $`P`$ behavior is displayed in the bottom panel of Fig. 3; it emerges that the thermodynamic integration of Eq. (3) fully agrees with the direct estimate of the chemical potential based on the Widom technique.

Starting from the knowledge of the free energy along the isotherm $`\overline{T}=2100`$ K, the free energy at different temperatures has been obtained through Eq. (4). We have examined the isochores $`\rho =1.25`$, 1.27, and 1.30 nm-3 in the solid phase and the density range $`\rho =[0.701.00]`$ nm-3 with steps $`\mathrm{\Delta }\rho =0.05`$ nm-3 in the liquid phase, descending down to $`T=1800`$ K, with temperature intervals $`\mathrm{\Delta }T=25`$ K. Simulations have been carried out with $`N=864`$ particles at constant density; four to eight cumulation runs of 10 000 steps at each state point are sufficient to yield accurate internal energies estimates (see Fig. 4), so to allow a smooth interpolation for the integration in Eq. (4). The pressure and the free energy along several isotherms are shown in Figs. 5 and 6, respectively; fully consistent results are obtained if we determine the free energy along an isochoric path first, and then integrating Eq. (3) at constant temperature. The accuracy of this global check is also evidenced in Fig. 6.

In order to analyze the thermodynamic properties of the system in the low density regime $`\rho =[0.050.20]`$ nm-3, we have calculated the chemical potential of the vapor phase along the isotherm $`T=1900`$ K through the Widom technique, and then estimated through Eq. (4) the free energy in the temperature range $`T=[18001900]`$ K.

As is visible in Fig. 7, where the chemical potential of the different phases is shown, liquid-solid equilibrium is stable at $`T=1900`$, while a solid-vapor coexistence takes place at $`T=1850`$ K. The intermediate temperature $`T=1875`$ K is characterized by almost a comparable value of the chemical potential of the various phases. We thus estimate that the triple point temperature is $`T_{\mathrm{tr}}1875`$ K, at the pressure $`P_{\mathrm{tr}}2.4`$ MPa; we then obtain from the equation of state $`\rho _{\mathrm{tr}}0.74`$ nm-3, in close agreement with the free energy results of Ref. , $`T=1880`$ K and $`\rho =0.74`$ nm-3. The features in the $`\mu `$ vs $`P`$ behavior which determine the narrow temperature width of the liquid pocket are clearly illustrated in Fig. 7, where it emerges that the solid and the liquid free energy branches have dissimilar slopes and considerably different spacing under equal variations of temperature. These two circumstances cause a fairly rapid shift of the liquid-solid intersection points and hence of the corresponding equilibrium parameters. Conversely, the vapor branch is hardly sensitive to temperature variations, and already at $`T1900`$ K tends to loose any further intersection with both the solid and the liquid branch. The only surviving liquid-solid equilibrium therefore fully characterizes the system behavior for $`T>1950`$ K.

The phase diagram of the system is displayed in Fig. 8, along with our previous Gibbs Ensemble Monte Carlo (GEMC) determination of the binodal line . Distinct liquid-vapor and liquid-solid equilibria take place at temperatures slightly higher than 1875 K; as far as the liquid-vapor coexistence is concerned, we observe a remarkable agreement between free energy calculations and GEMC results. At lower temperatures, the binodal points are metastable with respect to the freezing line; in this case as well, the GEMC approach fully reproduce the free energy data. We thus retain the GEMC results $`T_{\mathrm{cr}}1940`$ K, $`\rho _{\mathrm{cr}}0.42`$ nm-3, and $`P_{\mathrm{cr}}2.7`$ MPa, obtained in Ref.  with 1500 C60 particles, as a reliable estimate of the critical point parameters. We note in this context that an extended correponding-state behavior, which states a linear relationship between the critical temperature and the range $`R`$ of the interaction potential, has been recently outlined in Ref. . $`R`$ is defined by a map of the potential into an effective square well interaction with the same second virial coefficient. We obtain for the C60 model $`R=0.16`$, a value immediately over the minimum threshold of a stable liquid-vapor equilibrium, $`R=0.130.15`$ ; the extended rule’s estimate for the critical temperature is then $`T=1922`$ K, in fair agreement with the above $`T_{\mathrm{cr}}=1940`$ K result.

A comparison with the phase diagram determined by Hagen and coworkers  is reported in Fig 8. It has been conjectured, on the basis of a theoretical density-functional investigation , that the emerging discrepancy between Ref.  results and subsequent calculations might be due to an early cutoff of the C60 interactions, which substantially affects the location of the liquid-vapor binodal line. We argued on the other hand , that a fairly large simulation sample should be employed in this case, in order to take into account the peculiar density fluctuations in the GEMC simulations, which act to destabilize the liquid-vapor separation. We note that, while the overall fluid phase boundaries can sensitively depend on the interaction details (and hence on the cutoff), truncation effects play a minor role in the determination of the melting line; the latter appears indeed coincident with our estimate, and almost independent on temperature variations, ranging from $`\rho =1.26`$ nm-3 at $`T=2200`$ K to $`\rho =1.28`$ nm-3 at $`T=1800`$ K. Our results positively agree with the phase diagram obtained in Ref.  (also shown in Fig. 8), where a slightly higher critical temperature $`T=19541980`$ K is reported.

The phase diagram in the $`PT`$ representation is displayed in Fig. 9, where it appears that also the pressure range of the liquid phase is rather restricted, spanning only a few tens bar over the triple point pressure $`P_{\mathrm{tr}}=2.4`$ MPa.

We now turn to the indications on the freezing conditions obtained through the one-phase criterion expressed by Eq. (11). As is visible in Fig. 8, the $`\mathrm{\Delta }s=0`$ locus tends to overestimate the coexisting fluid density, thus affecting the location of the triple point in the phase diagram. We remark that a similar trend would emerge if the Hansen-Verlet prescription  for the height of the first peak of the structure factor were used . In fact, the close correspondence between the two interpretations of the freezing transition in terms of structural indicators, has been recently pointed out in Ref. . Similarly, the behavior of the internal energy, as well as the height of the first peak of the radial distribution function, show a non-monotonic behavior , suggesting that the system becomes unstable against the phase separation around the density $`\rho 1.0`$ nm-3 (see also Ref. ).

It thus appears that such indicators identify a very restricted range (if not a unique locus) of density vs temperature states, over which a structural reorganization of the fluid phase should be tendentially established, in order to satisfy purely excluded-volume, or equivalently entropic, demands. The almost vertical disposition of the $`\mathrm{\Delta }s=0`$ line (see Fig. 8), which means that the density keeps constant along the locus irrespective of the temperature, clearly reflects the substantial absence of any energy scale associated with such an indication. It is known at present that the vanishing of the residual multiparticle entropy accurately predicts the thermodynamic freezing threshold for a wide class of simple fluids, including hard-core models and the Lennard-Jones potential, in both three and two dimensions . In these models, the interparticle potential is dominated by steric effects, while the attractive forces can be treated as a perturbation to the inherent hard-sphere system, which essentially drives the liquid behavior, an approach commonly called the van der Waals picture of fluids . This is furthermore illustrated in Fig. 8, where it emerges that the reference hard-sphere system for the C60 model, obtained by splitting the potential into a repulsive and a perturbative part in the WCA fashion (after Weeks et al ), also freezes around the locus of vanishing residual multiparticle entropy, $`\rho (T)1.0`$ nm-3.

We argue that the discrepancy between one-phase indicators and full free energy calculations about the freezing transition of the system at issue, leads to a different scenario, where the strong attractive and rapidly decaying well in the interaction potential critically affects the conditions for the onset of the solid-liquid transition. The latter is driven in this case by “energetic” rather than by entropic effects, in agreement with several indications recently collected for systems with short-range interactions by A. A. Louis . It is shown in Ref.  that for this class of fluids most of the features exhibited by hard-sphere dominated systems do not arise, resulting in particular in the anticipation of the frezing threshold to lower densities than those predicted by solely structural conditions. In this respect, the lack of an accurate estimate of the freezing line of model C60 upon use of either the residual multiparticle entropy or the Hansen-Verlet prescriptions, must be interpreted less a “failure” of the criteria themselves, than a manifestation of their limited applicability in the present context, both indicators predicting the freezing threshold of the fluid on the basis of almost purely entropic requirements.

It appears in conclusion that the phase behavior of the Girifalco C60 model can be consistently understood in terms of general properties of systems interacting through short-range forces. Nevertheless, several issues are still open to further investigations and we refer firstly to the physical meaning to be associated to the fluid phase boundaries signalled by the above structural indicators. We note, on the other hand, the lack of clear indications of the freezing transition at a thermodynamic level, both in the internal energy and in the pressure, as documented in Figs. 4 and 5. More generally, the question of a detailed description of the microscopic behavior of the C60 model — as of other fluids interacting through short-range forces — already raised in Ref , is still unsolved. The formation of metastable clusters of strongly correlated particles in such fluids has been discussed in some recent papers , on the basis of the short-range nature of the interaction potential; however, in a recent molecular dynamics study  of quite a similar model, we have not been able to identify any net precursor, at a microscopic level, either of the fluid-solid threshold, or of the incipient crystallization of the liquid phase.

## IV Conclusions and perspectives

The free energy of the solid, liquid and vapor phases of the Girifalco model of C60 has been studied by means of extensive Monte Carlo simulations at constant density or pressure; the full phase diagram of the system has then been reconstructed on such a basis. It is confirmed by this comprehensive investigation that a stable liquid phase for this system exists, albeit confined to a rather restricted temperature range: we confidently estimate the triple and critical temperatures as $`T_{\mathrm{tr}}1875`$ K and $`T_{\mathrm{cr}}=1940`$ K, respectively. The pressure range of the liquid phase also appears rather narrow, spanning only a few tens bar over the triple point pressure. Our results illustrate how the interplay of various free energy branches determines the overall appearance of the phase diagram, and in particular the narrow extension of the liquid pocket.

The estimate of the freezing conditions, based on several structural indicators, is also critically discussed. It turns out that such indicators identify a practically unique thermodynamic locus, where the fluid phase becomes unfavoured due to entropic requirements. This locus almost coincides with the true thermodynamic freezing line for systems whose phase behavior is dominated by steric effects. For the model at issue, the solid-liquid transition is instead strongly affected by the deep, short-range attractive well in the interaction potential. As a result, the freezing transition of the fluid is driven to lower densities, mainly by energetic effects, in agreement with a scenario recently proposed by other authors.

As far as further studies are concerned, we are currently assessing, against the wide set of data produced in this work, the performances of several refined integral equation theories of the liquid state  and perturbation approaches, in order to describe on a full theoretical ground the phase diagram and the free energy properties of the Girifalco model. A refinement of our preliminary report on the phase diagram of higher order fullerenes Cn\>60  is also in progress. Finally, we plan to investigate the C60 properties, as well as the phase behavior of systems constituted by doped fullerites, in the low-temperature region of the solid phase.