text:
# THE FORMATION OF BROWN DWARFS AS EJECTED STELLAR EMBRYOS

## 1. INTRODUCTION

The past few years have seen a wealth of discoveries of brown dwarfs, whose bona fide nature as substellar objects have been confirmed by the presence of lithium in their spectra (Rebolo, Martín, Magazzù 1992; Martín, Basri, Zapatero Osorio 1999). Numerous brown dwarfs have been found in clusters like the Pleiades (e.g. Stauffer et al. 1998, Martín et al. 2000), some have been found in star forming regions like Taurus and Orion (e.g. Comerón, Neuhäuser & Kaas 2000; Luhman et al. 2000), and a few have been found as companions to main sequence or evolved stars (e.g. Nakajima et al. 1996). Brown dwarfs are now also rapidly showing up in infrared surveys of field stars like DENIS and 2MASS (e.g. Tinney, Delfosse & Forveille 1997; Kirkpatrick et al. 1999). The frequency of detection in both cluster and field star surveys is now so high that brown dwarfs may actually be as common as low mass stars. For a review, see e.g. Tinney (1999) and Basri (2000).

With such high abundances, it becomes even more important to understand the origin of brown dwarfs. It has been widely assumed that brown dwarfs form a natural extension to normal stars, i.e. large cores produce massive stars, and smaller cores produce lower-mass stars (e.g. Elmegreen 1999). Alternatively, it has been suggested that substellar objects may form through instabilities in circumstellar disks (e.g. Pickett et al. 2000). In a radical departure from these views, we here suggest that brown dwarfs are stellar embryos which have been ejected as a result of close dynamical interactions between small unstable groups of nascent stellar seeds, i.e. brown dwarfs differ from hydrogen burning stars only in that dynamical interactions deprived them from gaining further mass by prematurely cutting them off from their infalling gas reservoirs (Reipurth 2000). In this paper, we examine this suggestion in some detail, and propose various observational consequences and tests.

## 2. DISINTEGRATING MULTIPLE SYSTEMS

The dynamics of non-hierarchical triple or higher order systems has no analytical solutions and the chaotic motions of the members can only be followed numerically and analyzed statistically; for a review see Valtonen & Mikkola (1991). The motions can be divided into three categories: interplay, in which the members move chaotically among each other; close triple approach, when all members at the same time occupy a small volume of space; and ejection, in which one member departs the system after exchanging energy and momentum with the two others. A close triple approach is a necessary but not sufficient condition for ejection. Ejections often lead to escape, but can also result in the formation of a hierarchical triple system, with one body in an extended orbit. The remaining two members are bound closer to each other, forming a tighter and highly eccentric binary system. Most often, although not always, it is the lowest mass member that is ejected; the escape probability scales roughly as the inverse third power of the mass (e.g. Anosova 1986). The ejected member acquires a velocity $`v_{eject}15D_{ca}^{1/2}`$, where $`v_{eject}`$ is in km s-1 and the closest approach $`D_{ca}`$ is in AU (Armitage & Clarke 1997). Sterzik & Durisen (1995, 1998) have performed realistic numerical simulations of triple T Tauri systems and find ejection velocities of typically 3-4 km s-1 but with a higher velocity tail. The decay of a triple system occurs stochastically and can only be decribed in terms of the half-life of the process. Anosova (1986) finds that within about one hundred crossing times $`t_{cr}`$ almost all systems have decayed, where $`t_{cr}0.17(R^3/M)^{1/2}`$, and $`R`$ is a characteristic length scale for the system in AU, and $`M`$ is the total system mass in $`M_{}`$.

Three-body dynamical simulations including substellar members have been performed by Sterzik & Durisen (1999) and their dynamical results can be summarized as brown dwarfs are ejected because they are of low mass, in agreement with the general theory just outlined. This is in contrast to the evolutionary scenario advocated here which can be summed up as brown dwarfs are of low mass because they are ejected.

## 3. COMPETITION BETWEEN ACCRETION AND EJECTION

### 3.1. The General Problem

The mass of an individual star formed in a small multiple system is largely determined by the competition of two independent processes: on the one hand the accretion rate, which relates to the initial core properties, and on the other hand the decay of the multiple system, which leads to protostars being ejected and thus separated from their natal mass reservoir.

Consider a core containing a number of embryos whose initial mass is much less than the mass of distributed gas in the core. The timescale on which the bulk of the mass is accreted by the various embryos is of order the dynamical timescale of the parent core (Bonnell et al. 2001a), i.e. a few times $`10^5`$ years for canonical core properties. If the core mass divided by the number of embryos is larger than the hydrogen burning mass limit, then evidently if accretion went to completion and all embryos acquired their mass equally, no brown dwarfs would result. Brown dwarf production therefore requires at least one of the following conditions to be realized: (i) some embryos are ejected before the accretion process has gone to completion (i.e. the timescale for ejection is less than the timescale on which the embryo would otherwise acquire a mass equal to the hydrogen burning mass limit) and/or (ii) the embryos acquire their mass very inequitably.

A number of simulations, designed with rather different questions in mind, can be used to shed light on each of these scenarios.

In the first case, it is necessary for the ejection event to occur rapidly (i.e. on a timescale of less than the $`10^5`$ year timescale of the parent core). Therefore, the relevant embryos must form in a rather compact configuration (i.e. on a length scale that is much smaller than the core size). For example, Sterzik & Durisen (1995) showed that the half-life for ejection is around $`10^5`$ years if three embryos occupy a volume of diameter around 200 AU (i.e. on a size scale of around two orders of magnitude smaller than the typical size of the parent core.) The system then consists of a central mini-cluster of embryos which is fed mass from the infall of the surrounding envelope. \[To date, calculations of infall onto multiple systems have been restricted to co-planar hierarchical systems (Smith, Bonnell & Bate 1997), and so are not relevant to the current non-hierarchical situation\]. Numerical simulations of collapsing cores (e.g. Burkert & Bodenheimer 1993) show that the formation of a disk may be accompanied by the production of multiple fragments in non-hierarchical orbits, whereas star-disk and disk-disk interactions may likewise generate additional small fragments on the size scale of disks (Boffin et al. 1998, Watkins et al. 1998a,b). Fragmentation into even smaller volumes may be enhanced by magnetic fields (Boss 2000). Altogether, it would seem straightforward to achieve the initial conditions for the ejection model considered here.

In the second case, it is not necessary that the embryos form in a compact volume within the protostellar core, since an embryo does not have to be ejected in order to avoid reaching the hydrogen burning mass limit. \[One may instead envisage a situation similar to the ‘prompt initial fragmentation’ model proposed by Pringle (1989) in which fragments are seeded by existing structure in the core\]. In this case, it is enough if more gentle encounters perturb an embryo into an orbit where it does not intersect the densest part of the mass reservoir. Such a situation has an inbuilt feedback, since once an embryo has fallen behind in the race to acquire mass, it is more readily deflected away from the central densest regions. Likewise, embryos which acquire an early start in the mass race remain in the centre and preferentially intersect the infalling mass flow. This runaway to disparate masses is the basis of the competitive accretion scenario explored in a series of papers by Bonnell and collaborators (Bonnell et al. 1997, 2001a, 2001b). Of greatest relevance to the present discussion of small group dissolution are the simulations contained in Bonnell et al. (1997) which follow the evolution of cores containing 5-10 embryos initially distributed throughout the core volume. The trajectories and mass acquisition histories of the embryos clearly demonstrate the evolutionary path described above. Although these are only pilot simulations, and cannot be used to draw statistical inferences about the resulting mass distributions of stars, it is notable that the simulations do show a significant number of stars acquiring less than $`10\%`$ of the mean mass of stars in the clusters (corresponding roughly to the mass range of brown dwarfs: see, for example, Figure 2 of Bonnell et al. 1997 in which 2 out of the 10 objects end up in this category).

Although in the second case, an embryo can remain low mass without being ejected (on the mass acquisition timescale of $`10^5`$ years), the non-hierarchical situation is ultimately unstable, so that on timescales not much longer than this the system will break up into a system of binary and single stars. Thus in either case, the expected outcome is the production of a population of ‘stars’ whose mass is much less than the average mass produced by the core - i.e. objects we might loosely identify with brown dwarfs and very low mass stars. The critical difference between these two cases is the expected ejection velocity of the ”brown dwarfs”: for the case of the initially compact embryo configuration, the typical separations of close encounters are two orders of magnitude less than for the case that embryos are initially distributed throughout the core, and this translates into ejection velocities that are an order of magnitude greater (see Sects. 4.1 and 4.4).

The rough estimates presented here suggest that brown dwarfs could be common, maybe as common as normal stars, which is indeed what some recent surveys suggest. However, the only proper way to estimate the efficiency of this mechanism for producing brown dwarfs is through detailed numerical simulations of the dynamical evolution of small N-body systems in a time-variable potential and including the presence of circumstellar material. Though a number of simulations have attempted to capture aspects of this complex problem, there is clearly much more to be done in this area.

### 3.2. A Toy Model

To get a feeling for the time scales involved, consider the following order-of-magnitude estimates. Assume a flattened cloud is collapsing with a time-averaged mass infall rate $`\dot{M}_{infall}6\times 10^6(T/10K)^{3/2}`$ M yr-1, as derived by Hartmann, Calvet & Boss (1996). A 1 M cloud at 10 K would then take 1.7 $`\times `$ 105 yr to collapse. Assume that the collapse fragments into a number $`N_{mul}`$ of stellar embryos and that the infalling mass is distributed among them according to their mass ratios $`q_i=M_i/M_{prim}`$ (which may be time-dependent), where $`1iN_{mul}`$, and $`M_{prim}`$ is the time-dependent mass of the most massive embryo. Further assume that a part $`f`$ of the infalling gas is lost in outflow activity, where we adopt $`f`$$``$0.3. The growth rate of the $`i`$th embryo is then $`\dot{M}_i=\dot{M}_{infall}\times (1f)\times q_i/N_{mul}`$, and for the simplest case where $`N_{mul}=3`$ and $`q_i=1`$ it will then take a typical embryo about 6 $`\times `$ 104 yr to grow beyond a substellar mass.

If the three stellar embryos occupy a volume with diameter 200 AU, and if we consider a total mass in the range 0.06 to 0.24 M, then the characteristic crossing time $`t_{cr}0.17(R^3/M)^{1/2}`$ is between 1000 yr and 2000 yr. If, on the other hand, the configuration is tighter, occupying a volume with a diameter of only 20 AU, then $`T_{cr}`$ is merely 30-60 yr.

The usual decay equation is $`n_t/n_o=exp(0.693t/\tau )`$, where $`\tau `$ is the half-life of the decay. In their numerical simulations, Sterzik & Durisen (1995) found that about 95% of all systems had decayed after about 100 crossing times. This allows us to link the crossing time and the half life of the decay, $`\tau =23.1t_{cr}`$, so that to first order $`n_t/n_o=exp(0.18\dot{M}_{infall}^{1/2}(1f)^{1/2}N_{mul}^{1/2}R^{3/2}t^{3/2})`$.

While the stellar embryos are still very small, the crossing time, and therefore the half life, is extremely long, i.e. the embryos do not effectively start to interact until a time Ti when they have a certain mass. When the embryos are very small, their awareness of each other is limited because they are surrounded by the massive infalling envelope, and this further adds to $`T_i`$. But as the envelope material thins out and the embryo masses increase, their dynamical interactions become important. A more precise value of $`T_i`$ can only be determined by numerical experiments.

We will say that there is a reasonable chance that a brown dwarf is ejected if the half-life $`\tau `$ of the decay is less than $`T_{}T_i`$, where $`T_{}`$ is the free-fall time required to build up a multiple system of objects of average mass 0.08 M, i.e. $`T_{}=(N_{mul}\times 0.08)/\dot{M}_{infall}`$. If the ejection occurs at a time later than $`T_{}`$, then the ejected object is a star. In either case it has been assumed that the amount of circumstellar material the ejected object brings with it is small.

Figure 1 shows a simple, schematic presentation of these processes in a growth vs decay diagram, based on the values already discussed. The stapled line represents the growth of an embryo, which reaches the stellar mass threshold after $`6\times 10^4`$ yr. The solid curve shows $`n_t/n_o`$, which is a measure of the probability that the multiple system has not yet decayed. To be conservative, we have assumed the embryos occupy a volume of diameter 200 AU, resulting in a large $`t_{cr}`$. The time when interactions begin to take place, $`T_i`$, has arbitrarily been set to $`3\times 10^4`$ yr. In the particular example shown, about one third of a population of triple systems would have decayed before the expulsed member reached the stellar mass limit.

An important issue concerns the rate at which mass is accreted. The mass accretion rate adopted in Figure 1 is constant, in accordance with Shu (1977), but studies of shock-induced collapse suggest that infall might start with a powerful burst, briefly reaching values as high as 10-4 M$`{}_{}{}^{}yr^1`$, and subsequently decreasing (e.g. Boss 1995). If so, an embryo might much more quickly reach and exceed the stellar mass limit, leaving less time for the ejection of brown dwarfs.

This is, however, not a problem for the ejection scenario, because each of the following four physical processes will allow more time for ejection. Firstly, the infalling mass is distributed among the stellar seeds, and if $`N_{mul}`$ is higher than three, $`T_{}`$ becomes correspondingly higher. Secondly, in the calculations above we have for simplicity assumed that the components of newborn multiple systems have equal access to growth from infalling material, but as discussed in Sect. 3.1, competitive accretion is likely to be a very important process (Bonnell et al. 1997), leading to $`q_i`$ different from 1, and this would leave one or more embryos as substellar objects for extended periods of time. Third, if embryos are born in a more compact rather than a wide configuration (Sect. 3.1), then the crossing time decreases dramatically, leading to enhanced probability of very rapid ejection, on timescales of just a few thousand years. Finally, it is conceivable that stellar embryos may be formed over a period of time, with late-forming embryos being rapidly ejected by their more evolved siblings.

## 4. OBSERVATIONAL CONSEQUENCES OF THE EJECTION MODEL

In the following, we outline a number of issues relevant to current studies of brown dwarfs as seen from the perspective of the ejection model.

### 4.1. Detecting Brown Dwarfs near Embedded Sources

All stars, as they form initially as stellar embryos and rapidly gain mass through infall, must at some early point pass from substellar to stellar mass objects. Therefore, if isolated brown dwarfs are the result of the disintegration of a system of multiple stellar embryos, the moment of decay must take place very early in the infall process, and it thus appears likely that brown dwarfs are ejected during the Class 0 phase.

It follows that in order to find the very youngest brown dwarfs one should search near Class 0 sources. The number of Class 0 sources known is limited, and suggests that any attempt to compare decay calculations with observations will be severely limited by small number statistics. To get a practical sense of the observed separation between a nascent brown dwarf and its siblings, assume that it is observed at a time $`t`$ \[yr\] after ejection and moving with the space velocity $`v`$ \[km s-1\] at an angle $`\alpha `$ to the line-of-sight in a star forming region at a distance $`d`$ \[pc\]. Then the projected separation $`s`$ in arcsec is $`s=0.21vtd^1sin\alpha .`$

As an example, assume that a brown dwarf is moving out of a nearby ($`d`$ $``$ 130 pc) cloud at an angle of 60o to the line-of-sight with a space velocity of 3 km s-1. At a time of only 104 yr after ejection, the brown dwarf is already 42 arcsec from its site of birth, and after 105 yr, when the embedded phase of the original Class 0 source is about to end, the brown dwarf is 7 arcmin away, and its origin is soon lost in the mist of time.

Half of all ejected brown dwarfs will move into the cloud from which they formed. Such objects will be detectable only as highly extincted and weak infrared sources.

The rather high ejection velocities quoted above correspond to brown dwarfs ejected from a compact configuration (see Sterzik & Durisen 1998) and one sees that one has little chance of associating such ejected objects with the Class 0 star at an age comparable with the average age of Class 0 sources. If, however, brown dwarfs are ejected by the softer interactions envisaged in the competitive accretion scenario (see Section 3) then their ejection velocities are typically smaller, due both to the dissipative effects of gas drag and to the larger mean encounter distances in this case. The ejection velocities have not been quantified in this case, but rough estimates suggest that they may be only of order a few times $`0.1`$ km s-1. In this case, therefore, one would expect to detect one or more brown dwarfs around any given Class 0 object.

### 4.2. Binary Brown Dwarfs

As more brown dwarfs become known and more detailed studies of individual objects are undertaken, the number of known brown dwarfs paired in binaries is certain to increase (e.g. Martín et al. 2000; Reid et al. 2001). Among the first identified bound pairs of brown dwarfs is PPl 15, which has an orbital period of six days (Basri & Martín 1999), and other pairs have been recently reported (e.g. Martín, Brandner & Basri 1999).

Intriguingly, there is gathering evidence that whereas close pairs may be common (comparable with the binary fraction for solar type stars), there is a striking absence of wider pairs: to date, no brown dwarf binaries have been detected with separations larger than 23 AU, despite a lack of any obvious selection bias against the detection of such binaries (Martín et al. 2000). The fact that this apparent dearth of wider brown dwarf pairs is found both in the field and in clusters, hints that its origin is not due to the larger scale environment, but may reflect smaller scale processes that are common to both environments. In the context of the present paper, we would argue that the small scale clustering is ubiquitous (regardless of whether the parent regions survive as large scale bound clusters) and that any observational signatures of this origin would be expected in both field and cluster environments. \[Note that although some young clusters such as the Orion Nebula Cluster show no hint of primordial sub-clustering (Bate, Clarke & McCaughrean 1999), there is ample time for such clustering to have dissolved over the 2 Myr age of the cluster (Scally, in preparation)\]

It is clear, however, that purely point mass dynamical interactions will not produce brown dwarf binaries. Such interactions give rise to a dynamical biasing (McDonald & Clarke 1993), such that the resultant binary is usually composed of the most massive two members of the mini-cluster (van Albada 1968). McDonald and Clarke showed that under this assumption binaries comprising two brown dwarfs would essentially never be formed, a result confirmed by the more recent numerical simulations of Sterzik & Durisen (1998).

The situation is, however, different in the case that cluster dissolution is associated with dissipative interactions between the stars (embryos). An obvious source of such dissipation is provided by the presence of extended disks around each object, which exert a drag on neighboring stars whose orbits intersect the disk (Clarke & Pringle 1991a, 1991b, Hall, Clarke & Pringle 1995). McDonald & Clarke (1995) modeled the disintegration of small multiple systems in which the dynamics was crudely modified by parameterized star-disk interactions. They showed that the dynamical bias is considerably weakened in this case: dissipative interactions harden the resulting binaries and allow for the survival of low mass pairs that would otherwise have been broken up in dissipationless environments. Nevertheless, the predicted binary fraction for brown dwarfs in these simulations is still low (around $`5\%`$). More realistic calculations (Delgado, in preparation) are required in order to investigate whether this represents a serious objection to the ejection model.

Whatever the origin of brown dwarf binaries in mini-clusters, they clearly have to be rather close in order to survive the ejection. More specifically, Sterzik & Durisen (1998) found that the binaries that survive had separations not exceeding roughly one third of the mini-cluster radius. The apparent observed dearth of binary brown dwarfs with separation greater than $`2030`$ AU would then certainly be consistent with the scenario that brown dwarfs were ejected from a compact configuration (i.e. within a volume of radius $`100`$ AU).

A further criterion necessary for the survival of brown dwarf binaries is that after formation they do not continue to accrete beyond the hydrogen burning mass limit. In order for this to be the case, it is probably necessary for the interaction giving rise to the binary to eject it in the process. Again, further simulations will clarify if this is feasible.

### 4.3. Brown Dwarfs as Companions to Normal Stars

Many of the earliest searches for brown dwarfs were done towards main sequence stars or white dwarfs, and their mostly negative outcome added to the early speculations that brown dwarfs are rare. We now know, of course, that brown dwarfs are not rare, only as companions to other stars, the “brown dwarf desert” (Marcy & Butler 1998). Complete references to these searches are given by Basri (2000), who notes that the incidence of brown dwarf companions to stars with masses of 0.5 M or more does not exceed 1%.

If brown dwarfs are formed as stellar embryos which were prematurely ejected, this observational result is readily explained. If these stellar embryos had not been ejected, they would have continued to accrete matter from the infalling gas, and eventually would have passed from the substellar to the stellar regime. Indeed, in the ejection model, the problem is not to explain the absence of brown dwarfs around more evolved stars, but rather to explain how any brown dwarf could be still bound in such a system, because continued accretion on to both of the binary components would raise the brown dwarf companion above the hydrogen burning mass limit. There are two different approaches to this.

In the first scenario, we must assume that binaries consisting of a brown dwarf and a main sequence star must have been ejected very early from the mini-cluster, in precisely the same manner discussed in Sect. 4.2 for pairs of brown dwarfs, except that in this case one of the components happens to have already exceeded the hydrogen burning limit by the time of ejection. Since statistically the lightest members of a mini-cluster are more likely to be ejected, it follows that binaries with a main sequence star and a brown dwarf are more likely to have a very low mass star like an M dwarf as a member than a more massive star like a G dwarf. Within the limited statistics available so far, this indeed seems to be the case (Rebolo et al. 1998).

In the second scenario, competitive accretion in a mini-cluster creates a group of objects with a wide range in masses. As the small N-body system dissolves through dynamical interactions, the formation and subsequent ejection of binary systems may be influenced by a dynamical bias against binaries with very low mass companions. In the case of purely point mass gravitational encounters, brown dwarfs should rarely be companions to normal stars (McDonald & Clarke 1993), since the relatively fragile nature of such extreme mass ratio binaries makes them vulnerable to exchange reactions with other more massive stars in the newborn multiple system. The observation of at least a small number of such pairs might imply a role for dissipative interactions which harden the resulting binary and permit its survival: in the simulations of McDonald & Clarke (1995), which employed very massive and extended disks, large numbers of low mass stars or brown dwarfs ended up as companions to solar type stars. The fact that, observationally, such pairs are not plentiful suggests a more limited role for dissipative interactions in real systems.

Because the lowest mass members of a multiple system have a much increased chance of being ejected, it follows that if giant planets orbiting other stars were formed as stellar embryos, then they should be exceedingly rare. This is obviously contrary to the observation of numerous giant planets around stars, and it suggests that brown dwarfs and giant planets must form in two different processes.

### 4.4. Kinematics of Brown Dwarfs

If brown dwarfs are ejected stellar embryos, they should carry kinematic evidence reflecting their origins in small multiple systems. To first order, the ejected member from such a system acquires a velocity comparable to the velocity attained at pericenter in the close triple encounter. Low mass ejected members will therefore generally have higher velocities than higher mass objects.

Whether such an ejection history is likely to be detectable, through radial velocity and proper motion measurements of brown dwarfs, depends on the magnitude of the ejection velocities (which themselves depend on the closeness of triple encounters) and the general level of stellar motions in the star forming complex. For example, in clusters such as Orion or the Pleiades, where the velocity dispersion of stars in the large scale cluster potential is several km s-1, brown dwarfs ejected at such speeds from primordial mini-clusters would leave no kinematic imprint. Note that the very existence of brown dwarfs in the Pleiades means that retained brown dwarfs must have been ejected at less than the cluster’s escape velocity; in practice, even in the compact simulations of Sterzik & Durisen (1998), most brown dwarfs would be retained in the Pleiades. There may, however, be a significant tail of brown dwarfs ejected at higher velocities, which would populate the outer reaches of the Pleiades. Numerical simulations of the Pleiades that do not include primordial sub-clustering (Fuente Marcos & Fuente Marcos 2000) indicate that very few brown dwarfs are expected to be ejected by subsequent encounters within the Pleiades. Hence the detection of a brown dwarf halo around the Pleiades (or objects with discrepant proper motions) would provide unambiguous evidence of their origin in very compact mini-clusters. Further brown dwarf searches in the vicinity of the Pleiades (e.g. Bouvier et al., in preparation) are therefore likely to be particularly valuable.

Further constraints on the kinematic history of brown dwarfs are provided by their spatial distributions within star formation regions. For example, Hillenbrand & Carpenter (2000) have shown that in the core of the Orion Nebula Cluster, the mass function turns over at around 0.15 M. Given that there is evidence for mass segregation amongst more massive stars in the ONC (Hillenbrand & Hartmann 1998), it is obviously of great interest to determine if the substellar mass function has a similar form at larger radii in the cluster. In the case of Taurus (where the stars are not mutually bound) it is tempting to speculate that the paucity of brown dwarfs recently reported by Luhman et al. (2000) may be due to their higher ejection velocities, so that they have already left the region within which the bulk of T Tauri stars are currently concentrated.

It should be emphasized that brown dwarfs in very massive and/or very old clusters have velocities that are higher because their kinematics have been dominated by two-body relaxation. Therefore, such clusters should also be surrounded by a halo of brown dwarfs, but not because of their high ejection velocities. A proper test of the ejection scenario is made only by finding a halo of brown dwarfs around clusters small enough or young enough that relaxation has not yet dominated their kinematics. With regard to old brown dwarfs in the field, these have not been dynamically relaxed, but the majority must have evaporated from large clusters, in the process losing their original kinematic history (Spitzer & Harm 1958).

### 4.5. Substellar Equivalents to the Classical T Tauri Stars

Classical T Tauri stars show emission line characteristics and veiling which identify them as young accreting objects (Herbig 1962). The source of this accretion is a circumstellar disk. If brown dwarfs are ejected stellar embryos, then their disks are likely to be pruned during the close interactions that ejected them. Armitage & Clarke (1997) showed that very close encounters, which truncate the disks to radii of a few AU, could promote the rapid decline of classical T Tauri star characteristics thereafter and one might on these grounds expect that all but the very youngest brown dwarfs are unlikely to show any of the spectroscopic signatures of infall and outflow which we have come to associate with extremely young stellar objects. This appears to be consistent with current observations of young brown dwarfs (Martín 2001). In the case of ejected brown dwarfs, an empirical upper limit on expected encounter separations may be provided by the observed lack of brown dwarf binaries at separations greater than $`2030`$ AU. Further modeling is required in order to quantify how quickly accretion characteristics would be expected to decline following brown dwarf ejection: we, however, point out that because of the draining of the pruned disk, and because ejected stars have no source of disk replenishment, it may well be difficult to find the substellar equivalents of classical T Tauri stars, except perhaps as newly ejected objects located in the vicinity of Class 0 sources.

A further issue is the effect of ejection and disk pruning on the rotation of brown dwarfs, since in many models the disk plays the major role in braking young stars (Königl 1991, Cameron & Campbell 1993, Armitage & Clarke 1996). One might thus expect that ejected brown dwarfs in star forming regions should be rapidly rotating (i.e. at speeds close to break up). However, it turns out that extremely close encounters, even at the sub-AU level, are required in order to significantly affect disk braking, at least in the case of more massive stars (see the discussion in Clarke & Bouvier 2000). Further calculations using brown dwarf parameters and plausible enounter distances, are required before one can make firm predictions of the importance of this effect, although the sign of the prediction (a tendency towards more rapid rotation) is clear.

### 4.6. The Low-Mass End of the IMF

Numerous studies have explored the origin of the initial mass function (see e.g. Kroupa 1995; Larson 1999; Elmegreen 2001). Such analyses have been greatly helped by the increasingly accurate determinations of the mass function reaching into the very low mass and substellar regimes (e.g. Reid et al. 1999, Luhman et al. 2000, Najita et al. 2000, Hillenbrand & Carpenter 2000, Muench, Lada & Lada 2000). These studies have firmly established that the standard Salpeter power law slope of the IMF at higher masses gives way to a much shallower slope at lower masses or even a turn-over. While the main power law slope probably is related to the hierarchical structure of interstellar clouds, the flat part of the IMF at lower masses may result from a second, unrelated, physical process (e.g. Elmegreen 2000).

The ejection model can naturally produce the observed turn in the IMF . With reference to Figure 1, the shape of the very low mass end of the IMF is significantly influenced by the difference between $`T_i`$ and $`T_{}`$. If $`T_i`$ should be greater than $`T_{}`$ all objects ejected would be low mass stars. If interactions occur even just slightly before $`T_{}`$, the steep decline of $`n_t/n_o`$ ensures that a population of brown dwarfs will come into existence. And if $`T_i<<T_{}`$, brown dwarfs will be abundant.

## 5. DO LOW MASS STARS FORM IN SMALL MULTIPLE SYSTEMS?

In recent years it has been increasingly well documented that binarity is common among pre-main sequence stars and probably higher than among main sequence stars (e.g. Reipurth & Zinnecker 1993, Ghez, Neugebauer & Matthews 1993, Köhler & Leinert 1998). More recently, an analysis of deeply embedded outflow sources has revealed an observed binary frequency of around 80%, with half in higher order systems (Reipurth 2000).

On the theoretical side, there is growing consensus that fragmentation during the protostellar collapse phase is the principal mechanism for forming binaries and multiple systems (e.g. Bodenheimer et al. 2000, Bonnell 2001, Whitworth 2001), although the formation of close binary stars may require additional steps, as for example dynamical interactions between the nascent stars.

For brown dwarfs to be formed as ejected stellar embryos, close triple approaches must occur very early in the infall process. This is greatly facilitated if from the outset stellar seeds are formed close to each other. In this context it is of interest that the latest three-dimensional calculations of fragmentation in a collapsing core including magnetic fields show that magnetic tension helps avoiding a central density singularity (Boss 2000). These models suggest that fragmentation leads to a transient quadruple system. This is the situation which is ideal for dynamical interactions resulting in the early expulsion of a stellar embryo.

Dynamical interactions between members of small multiple systems are well studied and it is incontrovertible that they sometimes must play a role during the star formation process. The issue is therefore not whether stellar embryos can be ejected from a nascent multiple system, but rather how often this process occurs. Larson (1972, 1995) has suggested that all stars may be born in binary or multiple systems, and the high observed multiplicity frequency of Class 0/I sources near $``$80% would seem to support such a scenario (Reipurth 2000). The detection and statistical study of newborn brown dwarfs in star forming regions may therefore cast light on the general question of how low mass stars form.

## 6. CONCLUSIONS

We have argued here that the primary route for brown dwarf formation is not through the collapse of low mass cores, but that rather they are ejected as a result of encounters in multiple systems comprising a small number of stellar embryos. A number of simulations give plausibility to this picture, though the expected ejection velocities (and other properties such as the presence of accretion diagnostics, the range of binary separations and brown dwarf rotation rates) will depend on whether the embryos form on size scales comparable with the parent core ($`10^4`$ AU; ‘prompt initial fragmentation’) or on size scales of protostellar disks ($`100`$ AU).

To date, simulations have tended to address the problem in piecemeal fashion, focusing separately on the roles of point mass dynamics, star-disk interactions and competitive accretion. A clear goal for theorists is more realistic simulations that integrate all these features, a goal that is becoming increasingly realizable with the aid of advances in massively parallel computer architecture. In the absence of such simulations, the predictions given above have necessarily been laced with numerous caveats. Rather than repeating these here, we set out the key observations in the coming years that will provide the datasets necessary to test the models:

1. Brown dwarf searches in the vicinity of Class 0 sources

A null result in this area would not contradict the ejection model, but would imply high ejection velocities (several km s-1) consistent with embryos being ejected from a compact configuration ($`100`$ AU). A positive result would provide strong evidence in favour of the hypothesis, but, if common, would imply rather slow ejection velocities.

2. Brown dwarf binaries

The apparent upper limit of $`2030`$ AU on the separation of brown dwarf binaries may contain an important fossil record of a history of close encounters. Hence the placing of this result on a firmer statistical basis is an important goal. Likewise, it is important to establish whether close brown dwarf pairs are indeed as common as current estimates suggest. Although the survival of such pairs is unproblematical in the ejection model, it is not clear whether it is easy to form such low mass pairs starting from non-hierarchical initial conditions.

3. Brown dwarfs as companions to normal stars

The result that such companions are rare, at least in the case of solar type stars, is now well established. An important goal for simulations will be to discover whether they can simultaneously reproduce this property together with a high fraction of brown dwarfs as binary primaries (point 2. above).

4. Brown dwarf kinematics and distributions around open clusters

The presence of brown dwarfs in open clusters places an upper limit of several km s-1 on the ejection velocities of brown dwarfs that are retained in the cluster. Brown dwarf searches in the vicinity of the cluster will, however, explore a possible tail of higher velocity ejections. The existence of such a brown dwarf halo around the Pleiades, for example, would provide strong evidence in favour of the ejection hypothesis, since simulations starting from smooth initial conditions (no mini-clusters) do not eject brown dwarfs to large distances. On the other hand, a null result would not rule out the ejection hypothesis but would constrain the velocity distribution of ejected objects.

5. Circumstellar diagnostics and rotation rates of young brown dwarfs

Ejected brown dwarfs have suffered encounters that will have pruned their circumstellar disks. In consequence, one might expect weaker accretion diagnostics and faster rotation in young brown dwarfs compared with young stars. Further work is required to quantify the predicted magnitude of these effects, whilst considerably more observational data is needed in this area.

We thank John Bally, Jerome Bouvier, Nick Gnedin, Richard Larson, Eduardo Martín, and Michael Sterzik for discussions and helpful comments, Ka Chun Yu for preparing the figure, and Alan Boss for a helpful referee report.


format: text

text:
# Gamma Ray Bursts statistical properties and limitations on the physical model

## 1 Introduction

Cosmological origin of GRB had been first suggested in \[Prilutsky and Usov(1975)\] soon after their discovery. The present model of cosmological GRB based on production of gamma quanta from neutrino collisions $`\nu +\overline{\nu }e^++e^{}`$ was first considered in \[Berezinsky and Prilutsky(1987)\]. The efficiency of transformation of the neutrino flux $`W_{\nu \overline{\nu }}6\times 10^{53}`$ ergs into gamma quanta was estimated as $`6\times 10^6`$, giving a pulse $`W_\gamma 3\times 10^{48}`$ ergs. It could explain the cosmological GRB only at rather narrow pulse beam. In the giant GRB 990123 the isotropic energy production is very large \[Kulkarni et al.(1999)\] in gamma $`W_\gamma 2.3\times 10^{54}`$ ergs, and in optics $`W_{\mathrm{𝐨𝐩𝐭}}10^{51}`$ ergs. Simultaneous strong beaming in gamma and optical bands is rather unplausible. Strong beaming would modify the observed smooth optical light curve in presence of a source rotation. Some problems in GRB interpretation and modelling are discussed.

## 2 Statistics and restrictions to the model

BATSE data start to deviate from the uniform distribution with 3/2 slope at rather large fluences, for which KONUS data are well defined. Analysis of KONUS data had been done in \[Higdon and Schmidt(1990)\]. Taking into account selection effects, the resulting value $`V/V_{max}=0.45\pm 0.03`$ was obtained. KONUS data had been obtained in conditions of constant background. Similar analysis \[Schmidt(1999)\] of BATSE data, obtained in conditions of substantially variable background, gave resulting $`V/V_{max}=0.334\pm 0.008`$. These two results seems to be in contradiction, because KONUS sensitivity was only 3 times less than that of BATSE, where deviations from the uniform distribution $`V/V_{max}=0.5`$ in BATSE data are still rather large \[Fishman and Meegan(1995)\].

In presence of a threshold deviations of $`V/V_{max}`$ from its uniform Euclidean value 0.5 may be connected with the errors in determination of the burst peak luminosity or total fluence \[Bisnovatyi-Kogan(1997)\]. Such errors may be connected with spectral differences, variable sensitivity of detectors for bursts coming from different directions, variable background. All these reasons lead to underestimation of the burst luminosity, and decrease the slope of the curve $`\mathrm{log}N(\mathrm{log}S)`$. There is no angular correlation between GRB and sample of any other objects in the universe.

From the energy conservation law it follows $`W<Mc^2`$, where $`M`$ is a mass of the source. A proper account of physical laws put much stronger restrictions to the energy output. Calculations of ns-ns collisions gave energy output in $`(X,\gamma )`$ region not exceeding $`10^{50}`$ ergs \[Ruffert and Janka(1998)\], and similar results characterize ns-bh collision \[Ruffert and Janka(1999)\]. Magnetorotational explosion does not give larger energy output in $`(X,\gamma )`$ region, transforming about 5% of the rotational energy into a kinetic one \[Ardeljan et al.(2000)\]. The problems with vaguely defined ”hypernova” model had been discussed in \[Blinnikov and Postnov(1998)\].

The largest $`\gamma `$-ray production efficiency, close to 100%$`Mc^2`$ may be expected, if GRB originate from matter-antimatter star collisions. That arises a problem of antistar creation in the early universe.

Simultaneous $`\gamma ,X`$ and optical observations in GRB, accompanied by spectral and polarization experiments are very important. Search of hard $`X`$-ray lines and of annihilation 0.511 keV, line declared by KONUS, remain to be a puzzle which should be solved. Cosmological GRB explosion in a dense molecular cloud would lead to a specific optical light curve \[Bisnovatyi-Kogan and Timokhin(1997)\], which discovery would reveal conditions in the region of cosmological GRB explosion. Study of hard $`\gamma `$ afterglow, similar to the one observed by EGRET \[Schneid et al.(1992)\] is expected in a near future.

## 3 Short GRB and SGR

All afterglows had been measured only for long bursts. It cannot be excluded that short bursts could have another, may be Galactic origin. There is also a possibility, that short bursts are connected with a giant bursts observed in 3 soft gamma repeaters (from 4 known). At larger distances only giant bursts, appeared as short GRB could be observed. If we accept the present interpretation of SRG, as galactic and LMC sources at distances 10-50 kpc, than only giant bursts should be visible in the nearest galaxies as weak (about few $`10^7`$) short GRB. Taking into account that Andromeda is $`4`$ times more massive than our Galaxy \[Vorontsov-Velyaminov(1972)\], we should expect \[Bisnovatyi-Kogan(1999)\] to see about 10 short GRB in its direction during the observation time, while no one was yet observed. Another large galaxy in the local group of galaxies Maffei IC 1805 is also more massive than the Galaxy, and short GRB from it are also expected.

Presently SRG are interpreted as young neutron stars with very strong magnetic field - ”magnetars” \[Duncan and Thompson(1992)\]. The estimation of the distance and, consequently, the luminosity is based on SGR identification with supernovae remnants (SNR), leading to large energy losses. This interpretation has several theoretical objections \[Bisnovatyi-Kogan(1999)\].

1. Hard gamma pulsars observed in 3 SGR have luminosities strongly exceeding the critical Eddington luminosity. At such luminosity a strong mass loss should smear out the pulses.

2. Rotational energy losses estimated from the period increase rate are much smaller than the observed gamma and $`X`$-ray luminosity even in a quiescent state. In the magnetar model the energy comes from the annihilation of magnetic field. Such annihilation should be accompanied by creation of energetic electrons and radio-emission. The radio-emission of SGR is very weak, its discovery is very difficult, and still not firmly established.

3. Giant bursts observed in 3 SGR at present interpretation are accompanies by a huge energy production, part of which should go into particle acceleration and kinetic energy outbursts. It should influence the near-by SNR, and produce a visible changes in radio and optics, similar to those produced by pulsar glitches in the Crab nebula, when much smaller amount of energy is released. No such changes had been reported up to now, probably because they have not been present there.

Another interpretation of SGR, free of these contradictions needs a smaller distance to SGR, what is possible if the connection with SNR would not be confirmed. Note, that all SGR are situated at the very edge, or even outside of the SNR envelope, requiring very high 1000-3000 km/s speed of the neutron star. Refusing this connection and suggesting $`1030`$ times smaller distance to SGR would remove the upper objections. The even smaller distances are less probable, because most SGR are situated in the galactic disk, and so should be situated at distances larger than this disk thickness. Existence of one SGR outside the galactic disk direction could indicate to its big age during which it could leave the galactic disk. Discovery of big population of neutron stars in the globular clusters and in the galactic bulge, as recycled pulsars, indicate to existence of neutron stars in the whole volume of the Galaxy.

## 4 Acknowledgments

This work was partially supported by RFBR grant 99-02-18180 and INTAS-ESA grant 120. The authors are grateful to Prof. J.C.Wheeler and the Organizing Committee for support and hospitality, and to B.Schaefer for discussion.


text:
# Flux dynamics and vortex phase diagram of the new superconductor 𝑀𝑔𝐵₂

\[

## Abstract

Magnetic critical current density and relaxation rate have been measured on $`MgB_2`$ bulks from 1.6 K to $`T_c`$ at magnetic fields up to 8 Tesla. A vortex phase diagram is depicted based on these measurement. Two phase boundaries $`H_{irr}^{bulk}(T)`$ and $`H_{irr}^g(T)`$ characterizing different irreversible flux motions are found. The $`H_{irr}^{bulk}(T)`$ is characterized by the appearance of the linear resistivity and is attributed to quantum vortex melting induced by quantum fluctuation of vortices in the rather clean system. The second boundary $`H_{irr}^g(T)`$ reflects the irreversible flux motion in some local regions due to either very strong pinning or the surface barrier on the tiny grains.

\]

The recently discovered new superconductor $`MgB_2`$ generates enormous interests in the field of superconductivity. Many important thermodynamic parameters have already been derived, such as the upper critical field $`H_{c2}(0)`$ = 13 - 20.4 T, the Ginzburg-Landau parameter $`\kappa `$ 26, and the bulk critical superconducting current density $`j_c2\times 10^6A/cm^2`$ at 4.2 K and 1 T in bulk samples and $`j_c8\times 10^4A/cm^2`$ at 4.2 K and 12 T in thin films. One big issue concerns however how fast the critical current will decay under a magnetic field and in which region on the field-temperature ( H-T ) phase diagram the superconductor can carry a large critical current density ( $`j_c`$ ). This $`j_c`$ is controlled by the mobility of the magnetic vortices, and vanishes at the melting point between the vortex solid and liquid. A finite linear resistivity $`\rho =(E/j)_{j0}`$ will appear and the relaxation rate will reach 100% at this melting point showing the starting of the reversible flux motion. In this Letter we present an extensive investigation on the flux dynamics by magnetic relaxation method. A vortex phase diagram will be depicted based on the magnetic relaxation data.

Samples investigated here were fabricated by both high pressure ( HP ) ( P = 6 GPa at 950 C for 0.5 hour ) and ambient pressure ( AP ) synthesis which was described very clearly in a recent publication. High pressure synthesis is a good technique for producing the $`MgB_2`$ superconductor since it can make the sample more dense and uniform ( in sub-micron scale ) and also prevent the oxidization of Mg element during the solid reaction. Our HP samples are very dense and look like metals with shiny surfaces after polishing. Scanning electron microscopy ( SEM ) shows that the HP sample is uniform in sub-micron scale but some disordered fine structures are found in 10 nanometer scale, being similar to the internal structure of large grains seen in the AP sample. X-ray diffraction ( XRD ) analysis on both type of samples show that they are nearly in a single phase with the second phase ( probably $`MgO`$ or $`MgB_4`$ ) less than 1 wt.%. For simplicity we present in this Letter only the results from the HP sample.

Resistive transition was measured by the standard four-probe technique and the magnetic measurements were carried out by a superconducting quantum interference device ( SQUID, Quantum Design MPMS 5.5 T ) and a vibrating sample magnetometer ( VSM 8T, Oxford 3001 ). To precisely calculate the critical current density $`j_c`$ the sample has been cut with a diamond saw into a rectangular shape with sizes of 3.2 mm ( length ) $`\times `$ 2.7 mm ( width ) $`\times `$ 0.4 mm ( thickness ). Fig.1 shows the diamagnetic transition of one of the samples measured in the field-cooled ( FC ) and zero-field-cooled ( ZFC ) process. In the FC process, the temperature was lowered from above $`T_c`$ to a desired temperature below $`T_c`$ under a magnetic field, and the data are collected in warming up process with field. Its signal generally describes the surface shielding current and the internal frozen magnetic flux profile. In the ZFC process, the temperature was lowered from above $`T_c`$ to a desired temperature below $`T_c`$ at a zero field and the data are collected in the warming up process with a field. Its signal generally describes the internal magnetic flux profile which is ultimately related to the flux motion. Resistive and diamagnetic measurement show that the transition is very sharp with a perfect diamagnetic signal and $`T_c^{onset}`$ = 40.3 K with the transition width only 0.4 K.

In Fig.1 we show the magnetization hysteresis loops ( MHL ) measured at temperatures ranging from 2 K to 38 K. The symmetric MHLs observed at temperatures up to 38 K indicate the dominance of the bulk current instead of the surface shielding current. The MHLs measured at low temperatures, such as 2 K, 4K show quite strong flux jump which will be discussed elsewhere. A small kink appearing in the field ascending process in low field region may be attributed to the Bean-Livingston surface barrier and will be presented separately. It is easy to see that the MHLs measured at low temperatures ( e.g., 2 K to 10 K ) are too close to be distinguishable. From these MHLs one can calculate $`j_c`$ via $`j_c=20\mathrm{\Delta }M/Va(1a/3b)`$ based on the Bean critical state model, where $`\mathrm{\Delta }M`$ is the width of the MHL, V, a and b are the volume, length and width of the sample, respectively. The result of $`j_c`$ is

shown in Fig.2. It is clear that the bulk critical current density $`j_c`$ of our sample is rather high. For example, at T = 18 K and $`\mu _0H`$ = 1 T, we have $`j_c=1.1\times 10^5A/cm^2`$, which is among the highest values that have been reported in bulk samples. The $`j_c(H)`$ curves have been measured with three different field sweeping rate 200Oe/s, 100 Oe/s and 50 Oe/s. It is interesting to note that for all MHLs, there are small tails ( shown in the inset to Fig.2 ) in high field region. Accordingly the $`j_c(H)`$ curves show also a small tail in high field region. Since $`j_c`$ value in the small tail region is almost $`4\times 10^4`$ times smaller than that in the low field region, therefore it is safe to conclude that this small tail is due to some secondary effect, such as some local regions with very strong pinning or the surface pinning by the tiny grains. From the contour of $`j_c`$ vs. H shown in Fig.2 one can see that there are two ways to determine the so-called irreversibility line $`H_{irr}`$. The first way is to take a criterion for $`j_c`$ = 30 $`A/cm^2`$ just before the appearance of the tail. This is named as the bulk irreversibility line $`H_{irr}^{bulk}(T)`$ which signals the sharp drop of $`j_c`$ before the tail and reflects the irreversible flux motion in clean bulk samples. The second way is just to follow the small tail to a criterion of $`3A/cm^2`$. In this way one can determine a higher irreversibility line $`H_{irr}^g`$ which marks a boundary above which a complete reversible flux motion occurs. As mentioned above, the small tail in our present samples is due to a secondary effect, therefore the major part of the vortex system melt at $`H_{irr}^{bulk}`$. Thus we may conclude that the $`H_{irr}^{bulk}`$ reflects the melting of the vortex matter in the clean limit, for example, in single crystals. If more pinning centers are introduced into the sample, one should be possible to push the IL more close to the upper critical

field $`H_{c2}(T)`$. The ” irreversibility lines ” by following these two different methods are shown in Fig.3. It is clear that the lower phase boundary $`H_{irr}^{bulk}`$ terminates at about 8 T at 2 K, while the extrapolated values ( $`H_{irr}^g`$ ) reaches about 12 T at 2 K. There is a large separation between these two phase boundaries. It is important to note that Larbalestier et al. have found the similar feature where they regard the lower phase line $`H_{irr}^{bulk}`$ as the Kramer line. Fig.3 shows also the upper critical field $`H_{c2}(T)`$ determined from the temperature dependent magnetization by defining the $`H_{c2}(T)`$ as the point at which the magnetization start to deviate from the normal state linear background and that determined from the resistive measurement by Takano et al. on HP samples. A striking result is that the $`H_{irr}^{bulk}(T)`$ extrapolates to a rather low field, here for example, $`H_{irr}^{bulk}(0)`$ 8 T, while the $`H_{c2}(T)`$ extrapolates to a much higher value ( $``$ 15 T ) at zero K. There is a large separation between the two fields $`H_{c2}(0)`$ and $`H_{irr}(0)`$. If following the hypothesis of the vortex liquid above $`H_{irr}`$, we would conclude that there is a large region for the existence of a quantum vortex liquid at zero K. This can be attributed to a quantum fluctuation effect of vortices in bulk $`MgB_2`$. Although the lowest temperature in our present experiment is 1.6 K, however, from the experimental data one cannot find any tendency for $`H_{irr}^{bulk}(T)`$ to turn upward rapidly to meet the $`H_{c2}(0)`$ at zero K. One may argue that the $`H_{irr}(T)`$ probably can

be increased to higher values by introducing more pinning centers into this sample. This is basically correct since recently a higher irreversibility line $`H_{irr}(T)`$ has been found in thin films and bulk samples irradiated by heavy ions. This indicates that the quantum vortex melting in our present bulk sample is due to the strong quantum fluctuation which smears the perfect vortex lattice leading to the vanishing of the shear module $`C_{66}`$ of the vortex matter ( probably within grains ). Dense disorders will strengthen the shear module and thus enhance the bulk irreversibility line. The irradiation of heavy ions by Perkins et al.  did not suppress but strongly increase $`j_c`$ at a high field would suggest that the low value of $`H_{irr}^{bulk}(T)`$ measured in unirradiated bulks ( like our samples here ) is not due to the weak links since otherwise the $`j_c`$ value would drop even faster with increasing the magnetic field after the irradiation. This indicates that the rather low $`H_{irr}^{bulk}(T)`$ observed in our present work reflects the melting of the vortex matter in sample with relatively pure structural details rather than the breaking of the Josephson couplings between grains. In addition, as mentioned above, the $`H_{irr}^{bulk}(T)`$ found in present work is very close ( or identical ) to those by other authors found in bulks, therefore it shows a more intrinsic feature corresponding to the melting of vortex matter in the rather pure system. All these imply a strongly linked current flow between grains.

In order to investigate the flux dynamics in the vortex solid state below $`H_{irr}^{bulk}`$ we have carried out the dynamical

relaxation measurement. According to Schnack et al. and Jirsa et al., in a field sweeping process, if the field sweeping rate is high enough, the quantity $`Q=dln\mathrm{\Delta }M/dln(dH/dt)`$ is identical to the relaxation rate $`S=dlnM/dlnt`$ measured in the conventional relaxation method, where Q is called as the dynamical relaxation rate, $`\mathrm{\Delta }M`$ the width of the MHL, $`dH/dt`$ the field sweeping rate. The raw data with three different sweeping rate ( 200 Oe/s, 100 Oe/s, and 50 Oe / s ) are shown in Fig.2. The Q values vs. field for different temperatures are determined and shown in Fig.4. It is clear that the relaxation rate increases monotonically with the external magnetic field and extrapolates to 100% at the bulk melting point $`H_{irr}^{bulk}`$. A closer inspection to the data finds that the rising of the relaxation rate Q vs. field is slowered down when it runs into the small tail region. Since the small tail is a secondary effect as mentioned before, we leave it to a future investigation. Here we concentrate on the flux dynamics before the setting in of the small tail. At 2 K it is found from the Q(H) data that the melting field is about 8 T, being very close to $`H_{irr}^{bulk}(T=2K)`$ determined from the $`j_c(H)`$ curve. It is known that the $`H_{irr}^{bulk}(T)`$ is rather stable in low temperature region, therefore we can anticipate a rather low value of $`H_{irr}^{bulk}(0)`$ which is below 9 T being much lower than $`H_{c2}(0)`$. The large separation between $`H_{irr}^{bulk}(0)`$ and $`H_{c2}(0)`$ may manifest the existence of the quantum vortex liquid due to strong quantum fluctuation of vortices in the pure system of $`MgB_2`$. Fig.5 shows the temperature dependence of the relaxation rate Q. The arrows point at the irreversibility temperatures at the corresponding fields $`H_{irr}^{bulk}(T)`$. It is clear that the relaxation rate in wide temperature region keeps rather stable against the thermal activation and fluctuation. However, when the bulk melting point $`H_{irr}^{bulk}(T)`$ is approached the relaxation

rate will quickly jump to 100%. This may indicate that the thermal fluctuation is not the dominant process for the flux depinning in the new superconductor $`MgB_2`$. It shows high possibility for the vortex quantum melting even at a finite temperature. Worthy to note is that the quantum tunneling rate is extremely low at an intermediate field, such as Q = 0.7% at 2 K and 2 T, but is rather high at a high field, for example Q = 33% at 2 K and 6 T. This may imply that the field will greatly enhance the vortex quantum fluctuation and tunneling.

In conclusion, in rather pure samples of $`MgB_2`$, the irreversibility field is rather low comparing to the upper critical field in low temperature region showing the possible existence of the quantum vortex liquid due to strong quantum fluctuation. The temperature and field dependence of the relaxation rate may further suggest that the vortex melting at a finite temperature is also induced by strong quantum fluctuation in pure systems, such as single crystals and bulks. A small tail following the $`j_c(H)`$ curve in high field region is attributed to a secondary effect induced by strong pinning or surface barrier on some tiny grains. This opens a way for the higher irreversibility line when more pinning centers are introduced into the pure sample.

###### Acknowledgements.

This work is supported by the National Science Foundation of China (NSFC 19825111) and the Ministry of Science and Technology of China ( project: NKBRSF-G1999064602 ). HHW gratefully acknowledges Prof. B. Ivlev and Dr. A. F. Th. Hoekstra for fruitful discussions, and continuing financial support from the Alexander von Humboldt foundation, Germany.


text:
# Cluster Emission and Phase Transition Behavior in the Nuclear Disassembly

## I Introduction

Phase transition and critical phenomenon is an extensive debating subject in nature sciences. Recently, the same concept was introduced into the microscopic systems, such as in atomic cluster Bert?  and nuclei Finn82 . The break-up of nuclei due to violent collisions into several intermediate mass fragments, can be viewed as critical phenomenon as observed in fluid, atomic, and other systems. It prompts the possible signature on the liquid gas phase transition in the nuclear system. One fact is that the sudden opening of the multifragmentation Biza93  and vaporization Rive96  channels can be interpreted as the signature of the boundaries of phase mixture Gros90 . In addition, the fact that the nuclear caloric curve Poch95 ; Haug96 ; Mayg97 ; Serf98 ; Tsang97 ; Xi97  in a certain excitation energy range shows a plateau gives a possible indication of a first order phase transition Poch95  as predicted in the framework of statistical equilibrium models Bond98 . On the other hand, the extraction of critical exponents in the charge or mass distribution of the multifragmentation system Gilk94  can be explained as an evidence of the phase transition.

Several theoretical models have been developed to treat such a phase transition in nuclear disassembly. These models are percolation model Camp88 ; Baue88 , lattice gas model (LGM) and statistical multifragmentation model Bond95  etc. In this work, what we are interested in is the lattice gas model. It is a simple short range interaction model Yang52 , but it can be successfully applied to nuclear systems with isospin symmetry and asymmetry Jpan95 ; Jpan96 ; Camp97 ; Carmona ; Gupta97 ; Sray97 ; Biro86 ; Sama ; Mull97 ; Gulm98 ; Ma99 . LGM is carried assuming the system is in a ”freeze-out” density $`\rho _f`$ with thermal equilibrium at temperature $`T`$. Previous calculations Gupta97  with LGM showed that there exists a phase transition for the finite nuclear systems by studying the effective power law parameter ($`\tau `$) of cluster mass or charge distribution Fish67 ; Lit93 , their second moments ($`S_2`$) Camp88 , and the specific heat. More recently, we proposed two novel criteria, namely multiplicity information entropy ($`H`$) and nuclear Zipf’s law to diagnose the onset of liquid gas phase transition in the framework of the isospin dependent LGM (I-LGM) and isospin dependent classical molecular dynamics (I-CMD) model Maprl .

In this work, we will show that the emission rate of clusters is a useful tool to diagnose the nuclear liquid gas phase transition, and while the largest fluctuation of cluster multiplicities is simultaneously revealed at the point of the phase transition by investigating the features of light particles, charged particles, intermediate mass fragments and the largest fragment of disassembling source in the frameworks of I-LGM and I-CMD.

The paper is organized as following. Sec.II gives the descriptions of I-LGM and I-CMD. Results and discussions are shown in Sec. III where the multiplicities of cluster emissions, their slopes and their fluctuations are investigated. The influence of the ”freeze-out” density on cluster emission is also presented in the framework of I-LGM, and comparisons of the results between I-LGM and I-CMD at a given ”freeze-out” density is followed. Finally the conclusion is given in Sec.IV.

## II Description of Models

### II.1 Isospin dependent lattice gas model

The lattice gas model of Lee and Yang Yang52 , where the grand canonical partition function of a gas with one type of atoms is mapped into the canonical ensemble of an Ising model for spin 1/2 particles, has successfully described the liquid-gas phase transition for atomic system. The same model has already been applied to microscopic nuclear system. Detailed descriptions for the features of LGM can be found in some literatures, eg. papers of Pan and Das Gupta et al. Jpan95 ; Jpan96 , of Campi and Krivine Camp97 , of Carmona, Richert and Tarancon Carmona  and of Chomaz and Gulminelli Gulm98  etc. Here we will adopt the model developed by Pan and Das Gupta Jpan95 ; Jpan96 . In order to better understand the context of this work, the models are described below.

In the LGM, $`A`$ (=$`N+Z`$) nucleons with an occupation number $`s`$ (i.e. ”spin”) which is defined as $`s`$ = 1 (-1) for a proton (neutron) or $`s`$ = 0 for a vacancy, are placed in the $`L`$ sites of 3 dimension lattice. Nucleons in the nearest neighboring sites interact with an energy $`ϵ_{s_is_j}`$. The hamiltonian of the system is written as

$$H=\underset{i=1}{\overset{A}{}}\frac{P_i^2}{2m}\underset{i<j}{}ϵ_{s_is_j}s_is_j.$$
(1)
The interaction constant $`ϵ_{s_is_j}`$ is related to the binding energy of the nuclei. In order to incorporate the isospin effect in the lattice gas model, the short-range interaction constant $`ϵ_{s_is_j}`$ is chosen to be different between the nearest neighboring like nucleons and unlike nucleons:

$`ϵ_{nn}`$ $`=`$ $`ϵ_{pp}=0.MeV,`$
$`ϵ_{pn}`$ $`=`$ $`5.33MeV,`$ (2)
which indicates the repulsion between the nearest-neighboring like nucleons and the attraction between the nearest-neighboring unlike nucleons. This kind of isospin dependent interaction incorporates, to a certain extent, Pauli exclusion principle and effectively avoids producing unreasonable clusters, such as di-proton and di-neutron clusters, etc. The disassembly of the system is to be calculated at an assumed ”freeze-out” density $`\rho _f`$ = $`(A/L)\rho _0`$, where $`\rho _0`$ is the normal nucleon density, beyond $`\rho _f`$ nucleons are too far apart to interact.

In the model, $`N+Z`$ nucleons are put in $`L`$ sites by Monte Carlo sampling using the canonical Metropolis algorithm Mull97 ; Metr53 . As pointed out in Ref. Carmona ; Shida , however, one has to be careful treating the process of Metropolis sampling in order to satisfy the detailed balance principle and therefore warrant the correct equilibrium distribution in the final state. Speaking in detail, in this work, first an initial configuration with $`N+Z`$ nucleons is established. Second, for each event, sufficient number of ”spin”-exchange steps is tested, eg. 20000 steps to let the system generate states with a probability proportional to the Boltzmann probability distribution with Metropolis algorithm. In each ”spin”-exchange step, a random trial change on the basis of the previous configuration is made. For instance, we choose a nucleon at random and attempt to exchange it with one of its neighboring nucleons or vacancies at random regardless of the sign of its ”spin” (Kawasaki-like spin-exchange dynamics Kawasaki ), then compute the change $`\mathrm{\Delta }E`$ in the energy of the system due to the trial change. If $`\mathrm{\Delta }E`$ is less than or equal to zero, accept the new configuration and repeat the next ”spin”-exchange step. If $`\mathrm{\Delta }E`$ is positive, compute the ”transition probability” $`W`$ = $`e^{\mathrm{\Delta }E/T}`$ and compare it with a random number $`r`$ in the interval $`[0,1]`$. If $`r`$ $``$ $`W`$, accept the new configuration; otherwise retain the previous configuration. 20000 ”spin”-exchange steps are performed to assure we get the equilibrium state. Third, once the nucleons have been placed stably on the cubic lattice after 20000 ”spin”-exchange steps for each event, their momenta are generated by Monte Carlo sampling of Maxwell Boltzmann distribution. Thus various observables based on phase space can be calculated in a straightforward fashion for each event. One important point of such Monte-Carlo Metropolis computations is that the above ”spin”-exchange approach between the nearest neighbors, independently of their ”spin”, is evidenced to be satisfied by the detailed balance condition as noted in Ref. Carmona ; Shida . In other words, this sampling method will guarantee that the generated microscopic states form an equilibrium canonical ensemble.

Once this is done the LGM immediately gives the cluster distribution by using the rule that two nucleons are part of the same cluster if their relative kinetic energy is insufficient to overcome the attractive bond Jpan95 :

$$P_r^2/2\mu ϵ_{s_is_j}s_is_j<0.$$
(3)
This method has been proved to be similar to the so-called Coniglio-Klein’s prescription in the condensed matter physics Coni80  and was shown to be valid in LGM Jpan95 ; Jpan96 ; Camp97 ; Gulm98 ; Ma99 .

### II.2 Isospin dependent classical molecular dynamics model

Since the lattice gas model is a model of the nearest neighboring interaction, a long range Coulomb force is not amenable to lattice gas type calculation. Pan and Das Gupta provided a prescription, based on simple physical reasoning, to decide if two nucleons occupying neighboring sites, form part of the same cluster or not Jpan96 . They first try to map the lattice gas model calculation to a molecular dynamics type prediction, both first done without any Coulomb interaction. If the calculations match quiet faithfully then they can study the effects of Coulomb interaction by adding that to the molecular dynamics calculation. Here we adopt the same prescription to use the molecular dynamics and therefore investigate Coulomb effect. The results and conclusions can now be compared and checked between the LGM and the CMD. Obviously, here we do not do any $`abinitio`$ molecular dynamics calculation but only use it for simulating the nuclear disassembly starting from a thermally equilibrated source which has been produced by the above I-LGM: i.e. the nucleons have been initialized at their lattice sites with the Metropolis sampling and have their initial momenta with Maxewell Boltzmann sampling. From this starting point we switch the calculation to a classical molecular dynamics evolution under the influence of a chosen force. Noting that in this case $`\rho _f`$ is, strictly speaking, not a ”freeze-out” density for molecular dynamics calculation but merely defines the starting point for time evolution. However since classical evolution of a many particle system is entirely deterministic, the initialization does have in it all the information of the asympototic cluster distribution, we will continue to call $`\rho _f`$ the ”freeze-out” density.

The form of the force in the CMD is chosen to be also isospin dependent in order to compare with the results of I-LGM. The potential for unlike nucleons is expressed as

$`v_{\mathrm{np}}(r)(r/r0<a)`$ $`=`$ $`C\left[B(r_0/r)^p(r_0/r)^q\right]exp({\displaystyle \frac{1}{(r/r_0)a}}),`$
$`v_{\mathrm{np}}(r)(r/r_0>a)`$ $`=`$ $`0.`$ (4)
where $`r_0=1.842fm`$ is the distance between the centers of two adjacent cubes. The parameters of the potentials are $`p`$ = 2, $`q`$ = 1, $`a`$ = 1.3, $`B`$ = 0.924, and $`C`$ = 1966 MeV. With these parameters the potential is minimum at $`r_0`$ with the value -5.33 MeV, is zero when the nucleons are more than 1.3$`r_0`$ apart and becomes strongly repulsive when $`r`$ is significantly less than $`r_0`$. We now turn to the nuclear potential between like nucleons. Although we take $`ϵ_{pp}`$ = $`ϵ_{nn}`$ = 0 in I-LGM, the fact that we don’t put two like nucleons in the same cube would suggest that there is short range repulsion between them. We have taken the nuclear force between two like nucleons to be the same expressions as above plus 5.33 MeV up to $`r`$ = 1.842 fm and zero afterwards:

$`v_{\mathrm{pp}}(r)(r<r0)`$ $`=`$ $`v_{\mathrm{np}}(r)v_{\mathrm{np}}(r_0),`$
$`v_{\mathrm{pp}}(r)(r>r0)`$ $`=`$ $`0.`$ (5)
The system evolves with the above potential. The evolution equations for each nucleon are, as usual, given by

$`\stackrel{}{p_i}/t`$ $`=`$ $`\mathrm{\Sigma }_{ji}_iv(r_{ij}),`$
$`\stackrel{}{r_i}/t`$ $`=`$ $`\stackrel{}{p_i}/m.`$ (6)
Numerically, the particles are propagated in the phase space by integrating Newton’s equations of motion through a ”leap-frog” algorithm. At asymptotic times, for instance, the original blob of matter has expanded to 64 times its volume in the initialization, the clusters are easily recognized: nucleons which stay together after arbitrarily long time are part of the same cluster. The observables based on cluster distribution in both models can now be compared. While they can be also compared by switching on/off Coulomb interaction within the molecular dynamics.

## III Results and Discussions

We choose a medium-size nucleus 129Xe as an example. The input parameters are temperature $`T`$ and ”freeze-out” density $`\rho _f`$ in the model calculations. In this work $`T`$ is mostly limited to the range of 3 to 7 MeV and the ”freeze-out” density $`\rho _f`$ is chosen to be the range of 0.18$`\rho _0`$ to 0.60$`\rho _0`$. In most cases, $`\rho _f`$ is chosen to be 0.38 $`\rho _0`$, since the experimental data can be best fitted by a $`\rho _f`$ between 0.3$`\rho _0`$ and 0.4$`\rho _0`$ in previous LGM calculations Jpan95 ; Beau96 , which corresponds to $`7\times 7\times 7`$ cubic lattice for Xe. There is also good support from experiment that the value of $`\rho _f`$ is significantly below 0.5$`\rho _0`$ Agos96 . In addition, 0.18$`\rho _0`$, corresponding to $`9\times 9\times 9`$ cubic lattice and 0.60$`\rho _0`$, corresponding to $`6\times 6\times 6`$ cubic lattice of ”freeze-out” densities for 129Xe are also taken to compare and check the results with different ”freeze-out” density. 1000 events are simulated for each combination of $`T`$ and $`\rho _f`$ which ensures good statistics for results.

### III.1 I-LGM in different ”freeze-out” density

Fig.1 shows that the mean multiplicities of emitted neutrons, protons, charged particles and intermediate mass fragments and the mean mass for the largest fragment evolve with temperature at different ”freeze-out” densities in the I-LGM calculation. At a fixed ”freeze-out” density, average neutron multiplicity ($`N_n`$), proton multiplicity ($`N_p`$), charged particle multiplicity ($`N_{cp}`$) and the largest fragment mass ($`A_{max}`$) display monotonously increasing or decreasing with temperature as expectedly. But the multiplicity ($`N_{imf}`$) of $`IMF`$ shows a rise and fall with temperature Ogil91 ; Tsang93 ; Mayg95 , when the system probably crosses the phase transition boundary. With the decreasing of the ”freeze-out” density, $`N_p`$, $`N_n`$, $`N_{cp}`$ and $`A_{max}`$ decreases since larger space separation among nucleons at smaller ”freeze-out” density makes the clusters less bound and therefore decreases the sizes of free clusters and then increases the cluster multiplicities. The situation of $`N_{imf}`$ is slightly complicate, i.e. it increases with the decreasing of ”freeze-out” density in lower temperature branch but it is contrary in high temperature branch.

It seems to be impossible to discover the possibility of phase transition of nuclei if we only see these mean quantities as shown above ($`N_{imf}`$ is an exception). However, when we focus on their slopes to temperature: figure 2, the sharp changes are observed at the nearly same temperature at each fixed ”freeze-out” density: for instance, 4 MeV at 0.18$`\rho _0`$, 5 MeV at 0.38$`\rho _0`$ and 6 MeV at 0.60$`\rho _0`$. At such a transition point, (1) the multiplicities of emitted clusters increase rapidly and after that the emission rate slows down; (2) the decrease of the largest fragment size reaches to a valley for such a finite system. Physically the largest fragment is simply related to the order parameter $`\rho _l`$ \- $`\rho _g`$ (the difference of density in nuclear ”liquid” and ”gas” phases). In the infinite matter, the infinite cluster exists only on the ”liquid” side of critical point. In a finite matter, the largest cluster is present on both sides of phase transition point. In this calculation, a valley for the slope of $`A_{max}`$ to temperature may correspond to a sudden disappearance of infinite cluster (”bulk liquid”) near the phase transition temperature. It is not occasional producing such waves of the slopes, it should reflect the onset of phase transition there. This idea can be supported by surveying the other phase transition observables: such as the effective power law parameter, the Campi’s second moment from the mass or charge distribution of fragment and the information entropy $`H`$ of event multiplicity distribution Ma99 ; Maprl . Figure 3 depicts these results. The minimum of $`\tau `$ and the maxima of $`S_2`$ and $`H`$ appear around respective phase transition temperatures at different values of ”freeze-out” density, $`i.e.`$ about 4 - 4.25 MeV at 0.18$`\rho _0`$, 5 - 5.5 MeV at 0.38$`\rho _0`$ and 6 - 6.5 MeV at 0.60$`\rho _0`$. These temperatures are consistent with the ones extracted from the choppy position of the above slopes. It indicates that the above slopes (emission rate) can be taken as a probe of a liquid gas phase transition of nuclei.

Furthermore, the largest fluctuations of cluster multiplicities are found around the phase transition point in the same calculation. Figure 4 illustrates that RMS width ($`\sigma `$) of the multiplicity distributions of neutrons, protons, charged particles and intermediate mass fragments, and the distribution of the largest fragment masses. These variances generally show peaks at the same phase transition temperatures as the ones extracted from the above observables for each fixed ”freeze-out” density. Noting that the fluctuation of $`A_{max}`$ is related to the compressibility of the system. These features are also consistent to the one of phase transition behaviors, i.e. there exists the largest fluctuation at the phase transition point Chas96 . This fluctuation represents internal feature of the disassembling system not numerical fluctuation.

We would like to point that the ”freeze-out” density dependent phase transition temperature, extracting from $`\tau `$ and $`S_2`$, was also observed in a previous work Jpan95 , but this $`\rho _f`$ dependent phenomenon will vanish when excitation energy is used as a variable Maepja . In other words, the excitation energy is perhaps a good correspondence of the critical temperature because of only one critical point and hence only one critical temperature for a system.

### III.2 Comparison of I-LGM to I-CMD in the same ”freeze-out” density

Considering the absence of long range Coulomb force in the LGM, we will adopt the CMD to investigate the Coulomb effect and check the features of cluster emission and its relation to the phase transition behavior. The same ”freeze-out” density of 0.38$`\rho _0`$ is used to make a comparison for the results of I-CMD and I-LGM. Figure 5 shows that $`N_p`$, $`N_n`$, $`N_{cp}`$, $`N_{imf}`$ and $`A_{max}`$ changes with temperature within both models. The multiplicities of clusters and the largest fragment mass are close each other between I-LGM and I-CMD except for the multiplicities of neutrons and protons, illustrating that the I-LGM is, in general, a good tool to describe the fragmentation if Coulomb interaction can be ignored. When Coulomb interaction is switching on, $`N_p`$, $`N_{cp}`$ and $`N_{imf}`$ increase due to the repulsive role among protons and while $`A_{max}`$ decreases. Meanwhile, $`N_n`$ does not change because of no Coulomb interaction. $`N_{imf}`$ also shows a rise and fall with temperature in the I-CMD cases. Coulomb force makes the turning temperature of $`N_{imf}`$ smaller due to the long range repulsion.

The slopes of multiplicities of emitted clusters and of mean mass of the largest fragment have been plotted as a function of temperature in case of I-CMD in figure 6. The definite peaks of slopes are found as in the I-LGM case. The corresponding temperature at peaks locates about 4 MeV in the I-CMD case with Coulomb interaction, and about 5 MeV in the I-CMD case without Coulomb interaction. This turning temperature should also reflect the onset of phase transition there. If we investigate $`\tau `$, $`S_2`$ and $`H`$: Figure 7, we will find that there are a minimum of $`\tau `$ and the maxima of $`S_2`$ and $`H`$ around phase transition temperatures, $`i.e.`$ about 4 - 4.25 Mev for I-CMD with Coulomb, around 5 MeV for I-CMD without Coulomb and around 5.5 MeV for I-LGM, respectively. Obviously, these phase transition temperatures are consistent with the ones extracted from the choppy position of the slopes of figure 6.

Finally the RMS widths of the multiplicity distributions of the clusters and of the largest fragment mass are checked in the I-CMD case in figure 8. The widths of the multiplicity distributions of neutrons and protons tend to be saturate at higher temperature, and while the ones for $`CP`$, $`IMF`$ and $`A_{max}`$ demonstrate peaks at a certain fixed temperature, i.e. around 4 MeV for the case I-CMD with Coulomb, 5 MeV for the case of I-CMD without the Coulomb, which is similar to the I-LGM case. These turning temperatures are also consistent to the phase transition temperature as shown in figure 6 and 7 in the I-CMD cases. Again, the largest fluctuation simultaneously appears in the point of phase transition.

Overall, the phase transition temperature seems to rely on some ingredients: such as the ”freeze-out” density (or pressure), the model and its interaction potential, but the rule that emission rate and fluctuation of cluster multiplicity as a probe of phase transition is not changed.

## IV Conclusions

In conclusion, features of the emissions of light particles, charged particles, intermediate mass fragments and the largest fragments are investigated for a medium size nucleus $`{}_{}{}^{129}Xe`$ in the frameworks of isospin dependent lattice gas model and classical molecular dynamics model. $`N_n`$, $`N_p`$, $`N_{cp}`$ and $`A_{max}`$ show monotonously increasing or decreasing but $`N_{imf}`$ shows a rise and fall with temperature. Slopes of these observables go through extrema at the same temperature where the largest fluctuation of cluster multiplicity distributions is also observed. This temperature is consistent to the phase transition temperature extracted from the extreme values of $`\tau `$, $`S_2`$ and $`H`$. It gives an indication that the cluster emission rate can be taken as a probe of the phase transition of nuclei and while the largest fluctuation is simultaneously accompanied when the onset of phase transition occurs. In addition, the systematical comparison of I-LGM and I-CMD shows that LGM is a rather good tool to study the nuclear disassembly when the system is not large where the Coulomb interaction can be ignored. In light of this work, we think that the experimental study about the cluster emission is rather meaningful, especially to measure the excitation function of the multiplicities, their slopes and variances, from which some signals for phase transition could be found.

ACKNOWLEDGEMENT

It’s my pleasure to thank Prof. S. Das Gupta and Dr. J. Pan for kindly providing the original codes and Prof. B. Tamain for helps. This work was partly supported by the NSFC for Distinguished Young Scholar under Grant No. 19725521, the NSFC under Grant No. 19705012, the Presidential Foundation of Chinese Academy of Sciences, the Science and Technology Development Foundation of Shanghai under Grant No. 97QA14038.