text:
# Modelling Self-Interacting CDM Haloes with a Cosmological Boltzmann Code

## 1 Introduction

The current Cosmology based on the growth of small fluctuations through gravitational instability in a universe with cold dark matter (CDM) defines a physically well-motivated theory, capable to explain most of the properties of the large-scale structures of the universe. In this scenario, where a large fraction of the matter is supposed non-dissipative, cold and collisionless, the dark haloes form hierarchically via gravitational collapse of primordial density fluctuations. The density profiles of these dark haloes have been studied systematically using N-body simulations (Navarro, Frenk $`\&`$ 1996, 1997; hereafter NFW). It was found that a simple analytical density profile with one free parameter is able to describe the typical structure of the virialized CDM haloes. On its own, this free parameter, which can be the concentration, depends only on the halo mass.

The NFW density profile goes as $`\rho r^1`$ at small radii. According to recent high-resolution simulations, the inner density profile of the haloes is even steeper than $`r^1`$ ($`\rho r^{1.5}`$) (Moore et al. 1999; Jing & Suto 2000). This prediction of the CDM scenario is in potential conflict with observations: the halo inner density profiles inferred from the HI rotation curves of dwarf and low surface brightness (LSB) galaxies and from a mass map constructed by strong lensing techniques for the galaxy cluster CL0024+1654 seem to be much shallower than $`r^1`$ (Moore 1994; Flores & Primack 1995; Burkert 1995; de Blok & McGaugh 1997; Tyson et al. 1998). Recently, the rotation curves of some of the dwarf and LSB galaxies were re-obtained using high-resolution $`H\alpha `$ observations, resulting the inner density profiles more concentrated than those obtained from the HI observations, but on average still shallower than $`r^1`$ and with a large scatter (Swaters et al. 2000; van den Bosch & Swaters 2000; Dalcanton & Bernstein 2000). Besides, as Firmani et al. 2000b have shown, the halo inner density profiles of most of the observed LSB galaxies could be shallower than reported if the baryon contraction is taken into account (see also Firmani & Avila-Reese 2000).

An interesting property to be confirmed with more observational studies, in particular for galaxy clusters, is that the central density of dark matter haloes is roughly independent of the halo mass with a value around 0.02 M pc-3 (Firmani et al. 2000a,b). Not only the predictions of the CDM scenario, but also of alternative scenarios as the warm dark matter one are unable to predict such a scale invariance of the central density of dark haloes (Avila-Reese et al. 1998; Colín et al. 2000).

In a recent burst of activity, several alternatives to the well motivated CDM scenario were proposed in order to overcome the difficulties mentioned above and others as the excess of satellites and a little angular momentum in the galaxy disc predicted by numerical simulations. These alternative theories try to suppress the excess of central mass of the galactic haloes and may be summarized into two classes: 1) a manipulation of the power spectrum with a suppression of the small scales as an intrinsic property of a dark matter: Avila-Reese et al. 1998; Moore et al. 1999; Firmani et al. 2000a,b; Peebles 2000; Colín et al. 2000. These works make different assumptions concerning thermal energy at the beginning of the mass aggregation history (MAH). 2) Assuming for dark matter particles a different nature from collisionless: self-interacting (Spergel $`\&`$ Steinhardt 2000); repulsive (Goodman 2000); decaying (Cen 2000); with a limiting phase-space density (Hogan $`\&`$ Dalcanton 2000) and annihilating (Kaplinghat et al. 2000). Intriguingly, only the assumption of a self-interacting or annihilating cold dark matter seems to be able to reproduce the observed central density scale invariance ranging from galaxy to galaxy cluster size.

In a previous work, adopting the self-interacting CDM scenario suggested by Spergel & Steinhardt, we have proposed a thermodynamical model able to describe the gravothermal expansion of the core in the haloes (Firmani et al. 2000a); with this model, the scale invariance of the halo central densities was predicted. The NFW density profile shows a velocity dispersion raising outwards with radius (Cole $`\&`$ Lacey 1996; Fukushige $`\&`$ Makino 1997). Therefore, if particles are self-interacting, the inward heat transfer induces a thermalization process able to avoid the formation of a cuspy core. The influence of weak self-interaction on the structure of the CDM haloes has been recently explored using cosmological numerical simulations (Yoshida et al. 2000; Davé et al. 2000). These works, where constant cross sections were used, have led to interesting results, that we will discuss in comparison with ours, later.

A cross section inversely proportional to the particle velocity dispersion seems to be better motivated by observations (Firmani et al. 2000a,b; Wyithe et al. 2000). In this paper we focus on detailed dynamical predictions of the density profiles of self-interacting spherical symmetric dark matter haloes ranging from galaxy to galaxy cluster mass. The cosmological model we assume is flat with vacuum density $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and expansion rate $`H_0=65`$ km s-1 Mpc-1.

## 2 The cosmological Boltzmann code

Our dynamical approach aimed to investigate the cosmological halo density profiles is based on a solution of the collisional Boltzmann equation for dark matter. The virialized dark halo is described by the distribution function $`N(E,J)`$, where $`N(E,J)\mathrm{\Delta }E\mathrm{\Delta }J`$ gives the number of particles within the intervals $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }J`$ centred at the total energy E and the angular momentum J. The particle population is assumed to be uniform with respect to the directions of J. The mass aggregation supplies new particles, each one with a total energy and angular momentum. The computation of each particle orbit is based on a gravitational field derived from the mass distribution in a self-consistent way. The change in time of the gravitational field due to the mass aggregation modifies the particle orbit in agreement with its adiabatic invariants. The particle mass is distributed along its orbit according to the fraction of the period spent in any interval. For collisionless particles this approach is similar to the secondary infall (Avila-Reese et al. 1998) and allows to calculate a virialized halo density profile for mass and angular momentum aggregation histories in a specific CDM universe. For a given universe, the mass aggregation history (MAH) is obtained making use of the extended Press-Schechter formalism based on the conditional probabilities for a Gaussian random field. For a given mass, we generate a set of MAHs through Monte Carlo simulations and we calculate the average MAH. If the perihelion to aphelion ratio of the orbits (angular momentum) is opportunely selected, the NFW profile is obtained on galaxy and galaxy cluster scales.

The effects of the collisions between particles is described by the collisional term of the Boltzmann equation and induce a further evolution on the distribution function. Our scheme will provide an accurate solution only when the collision time is sufficiently larger than the orbital period of the particles. The collision cross section multiplied the relative velocity of the colliding particles is assumed to be constant in agreement with the previous discussion in Firmani et al. 2000b. This hypothesis makes particularly efficient the Monte Carlo method implemented to describe the collisions of the particles of each discrete state (E,J). By each trial we move a large amount of particles from one state to another according to the probability of the process. The results show an enormous space resolution obtained with a minimum computational expense, with each run taking approximately a couple of hours on a 450MHz pentium processor.

The efficiency of our method is counterbalanced by the practical difficulty to explore solutions with collision times lesser than the particle orbital periods (high cross section), as well as halo structures which deviate from the spherical symmetry. For these cases complementary studies with N-body techniques have been planned.

## 3 Self-interacting dark matter halo evolution

We have calculated the evolution of haloes in a $`\mathrm{\Lambda }`$CDM cosmology assuming a self-interaction cross section per unit mass ($`\sigma `$) given by: $`\sigma v_{100}=10^{24}`$ cm2/GeV, where $`v_{100}`$ is the relative velocity of the colliding particles in units of 100 km s-1. Fig. 1 shows the present-day density profiles of haloes with masses of 6 1011 and $`\mathrm{6\; 10}^{15}`$ M (solid lines). For comparison, the corresponding NFW profiles are also shown (dotted lines). Remarkably, the self-interacting CDM haloes appear to be thermalized in the centre resulting in a non-singular isothermal sphere, while in the outer parts, the density profile is well described by the NFW shape which is representative of the hierarchical merger history. The non-singular isothermal profile obtained with our code confirms the halo profile (dashed line) we derived in a previous work in order to describe observations (Firmani et al. 2000b). In that case, the hydrostatic equilibrium was assumed in order to pass from a NFW profile for the outer regions to an isothermal sphere in the centre. For that “naive” model the core-radius was obtained once assigned the central density.

With the dynamical model based on the solution of the transfer equation we presented here, we can now estimate the typical number of collisions in the halo core over a Hubble time ($`t_H`$): $`N_{coll}=\rho \sigma vt_H`$. For a 6 1011 M halo and the assumed cross section, we find $`N_{\mathrm{coll}}4`$.

Fig. 2 shows the evolution of the density profile of an halo of 6 1011M. Since early epochs ($`z40`$) a soft core of few pc size is formed. At $`z=20`$ the central density is about $`4M_{}`$ pc-3 and the core size is already of $`40`$ pc. The central density is high enough to allow collisions on time scales of only a few hundred million years. Due to the mass aggregation, the inner radial velocity dispersion gradient is preserved. Then, due to heat transfers inwards, the soft core continues growing while the central density decreases. At any time, the core density is determined by a self-regulating mechanism which makes the collision efficiency able to transfer enough energy to expand the core. The central density decreases and the core size rises, both reaching values comparable with the present observations, i.e., $`0.03M_{}`$ pc-3 and 6 kpc, respectively.

Of great interest is the invariance of the central density, shown in Fig.1, in passing from an halo of galactic size to an halo with mass four order of magnitude larger, representative of a cluster of galaxies. In Fig. 3a and Fig.3b we compare the central densities and core radii predicted by the models with those inferred from observations of dwarfs (filled squares), LSB galaxies (empty squares) and the cluster CL0024+1654 (filled circle) (see Firmani et al. 2000b for details). Big empty squares are for those LSB galaxies with $`H\alpha `$ rotation curves (Swaters et al. 2000). The models were calculated for two cross sections: $`\sigma v_{100}=10^{24}`$ and $`\mathrm{2\; 10}^{24}`$ cm2 GeV-1 (dashed and dotted lines, respectively).

The agreement of the self-interacting CDM model with observations is rather remarkable. In a range of scales from dwarfs to galaxy clusters, the haloes have (1) a nearly constant central density, and (2) a core radius roughly proportional to the maximum circular velocity. We stress that the only way to obtain halo central densities independent of their masses is assuming $`\sigma 1/v`$.

It is interesting to note that constraining the density for the cluster to the value of $`\rho _c=0.02M_{}`$ pc-3 as inferred by the observation of CL0024+1654, we derive an estimate of the cross section: $`\sigma v_{100}10^{24}`$ cm2/GeV.

## 4 Comparison with previous works

In Fig.4 we show a grid of models for two different masses: $`\mathrm{6\; 10}^{11}M_{}`$ and $`\mathrm{6\; 10}^{15}M_{}`$ and wide range of cross sections. In order to scale the models with the mass for a given value of $`\sigma v_{100}`$, Fig.3 shows that $`\rho _c`$ is roughly scale invariant while $`R_c`$ is proportional to $`V_m`$.

With the aim to compare our models, derived with the assumption $`\sigma v^1`$ to models with $`\sigma =const`$, we will assume the latter velocity dispersion as representative of the collisional velocity. From Fig.2 of Yoshida and co-workers (2000), extrapolating the profile towards the center for the model S1W-a, with the same pattern of S1W-b and S1W-c, we read a central density of $`\rho _c=2\mathrm{3\; 10}^{16}M_{}`$/h2/Mpc3 and a core radius $`R_c=3550`$ kpc/h for a halo of $`\mathrm{7.4\; 10}^{14}M_{}`$/h. We have simulated a cluster of $`M=\mathrm{6\; 10}^{15}M_{}`$. Rescaling the mass of the Yoshida et al. (2000) halo ($`\mathrm{7.4\; 10}^{14}M_{}`$/h and $`v=1500`$ km/s) to the cluster mass used here ($`\mathrm{6\; 10}^{15}M_{}`$ and $`v=2400`$ km/s), it is easy to obtain $`\rho _c=0.0130.008M_{}`$/pc3 and $`R_c=85120`$ kpc. In our model, for the cross section of S1W-a, $`\sigma =0.1`$ cm2/g, we find $`\rho _c=0.01M_{}`$/pc3 and $`R_c=100`$ kpc. Let us now analyze the case S1W-b of Yoshida et al. (2000). In this model, which resolution is higher than the S1W-a case, the central density is $`\rho _c=\mathrm{4\; 10}^{15}M_{}`$h2 /Mpc3 and the core radius $`R_c=150`$ kpc/h. Making the same mass rescaling we obtain $`\rho _c=\mathrm{1.7\; 10}^3M_{}`$/pc3 and $`R_c=370`$ kpc. This case is for a cross section $`\sigma =1`$ cm2/g. In Fig.4 we see that our model for the same $`\sigma `$ predicts $`\rho _c=\mathrm{2\; 10}^3M_{}`$/pc3 and $`R_c=355`$ kpc. Taking into account the different approaches the agreement is really outstanding.

Regarding Davé et al (2000) work, the only case we have simulated is the first top panel from left to right of Fig.1 of the same paper. For $`\sigma =10^{24}`$ cm2 GeV-1 and for an halo of $`M=\mathrm{6.17\; 10}^{11}M_{}`$ they find $`\rho _c=0.03M_{}`$/pc3 in agreement with our value of central density of $`0.03M_{}`$/pc3.

## 5 Summary and conclusions

We have calculated the evolution of the density profiles of self-interacting CDM haloes using a numerical approach based on the collisional Boltzmann equation. This study presents two original points: (i) we were able to explore the halo density profiles covering about five orders of magnitude in the radius and for haloes from galaxy to galaxy cluster scales; (ii) we used a self-interacting cross section inversely proportional to the collision velocity of the particles. Analysing the results obtained with our code we conclude:

* The interior of dark haloes may be deeply affected by self-interaction. A modest cross section value ($`\sigma v_{100}10^{24}`$ cm2 GeV-1) is enough to produce soft cores.
* Collisions between dark particles induce a thermalization process with the consequent appearance of a soft core. Thermal equilibrium in the core is established with a few collisions in a Hubble time ($`4`$).
* Soft cores appear already at the very early dynamical evolution ($`z40`$); with time the central density decreases and the core size increases.
* The observed scale invariance of the central density from galactic to galaxy cluster scales is reproduced only if the cross section is inversely proportional to the collision velocity.
* The core radius predicted by the model increases proportional to the halo maximum rotation velocity in agreement with the observed trend.
* The observed central density of CL0024+1654 is consistent with assuming a self-interacting dark matter cross section given by: $`\sigma v_{100}10^{24}`$ cm2 GeV-1.

The code is limited in analysing dark halo properties which depart from spherical symmetry. Further studies concerning dark halo triaxiality and the evolution of substructures have been planned using N-body techniques.

## 6 Acknowledgments

We are grateful to Vladimir Avila-Reese for a critical reading of an early draft of the paper and to Carmelo Guzman, Liliana Hernández and Gilberto Zavala for computing assistance. ED thanks Fondazione Cariplo for financial support and is grateful to UNAM in Mexico City for its hospitality during the preparation of this work.


format: text

text:
# An Approach to the Implementation of Overlapping Rules in Standard ML

## 1 Introduction

Rule-based systems111In this paper, we focus exclusively on production rule systems. Related systems, such as those based on a notion of term rewrite rule, have not been considered at this point. have had a long history in AI and powerful implementations have been developed. The most problematic aspect of this work has always been that of integrating the rule-based approach to a general-purpose language for application support. As an example along those lines, in , Crawford *et al* describe R++, a rule-based extension to C++ completely integrated with the object-oriented features of the language, implemented as a rewrite of R++ into C++ code.

One common use of rules, and indeed the primary motivation for this work, is to help write rule-based reactive systems, reactive in the sense of having the system react to changes in the environment. A change in the environment should enable a certain number of appropriate rules that can be fired to react to the environment change, possibly effecting new changes to the environment that will enable other rules.

One aspect of rule-based reactive systems we focus on is that of long-acting overlapping rules. In most rule-based systems, a single rule is active at any given time: when the system determines that a rule is enable, the rule is fired, the environment is updated, and a new rule can be selected. This is not so great is some of the rules are computationally intensive: when such a rule is fired, it will take time to execute and if the system relies on other rules to ensure say responsiveness of the interface, the system will not respond to the user until the rule has finished executing. In general one, may want rules that span multiple other rules firing. For example, one may pre-fire a rule when certain conditions are met, perhaps pre-computing some values, and get ready for the “real” firing once the “real” conditions are in place. Or one may want a rule that when fired will perform a given action at every subsequent rule firing until maybe a condition occurs that stops this behavior. It is of course possible to achieve these effects in certain rule-based systems, by a process of chaining rules (at the end of one rule, setting up the firing conditions of the next rule), along with a suitable notion of concurrently fired rules.

In this paper, we describe an approach to the integration of rule-based programming in Standard ML (SML) , with a strong focus on overlapping rules to achieve the effects described above. The resulting framework is suitable for the design of domain-specific abstractions for various rule systems. In contrast with other work, we do not worry about efficiency issues in this paper, but rather concentrate on expressiveness and applicability of the framework. We assume throughout the paper a basic knowledge of SML, as described in various introductory material such as .

After reviewing the basic notions of rule-based programming in Section 2, we discuss the framework in Sections 34. We give an application of the flexibility of the framework in Section 5, where we design domain-specific abstractions for controlling a goal-directed agent. Section 6 focuses on implementation details, including a reference implementation in terms of an existing reactive library for programming reactive systems in SML, introduced in , and based on the reactive approach of Boussinot .

## 2 OPS-style production rules

In this section, we establish the terminology and the model of rules we are interested in, namely OPS-style production rules . The *production-system* model of computation is a paradigm on the same footing as the procedural paradigm, the functional paradigm or the object-oriented paradigm: it is a view of what a computation ought to be to best achieve a given goal.

A *production-system program* is an unordered collection of basic units of computation called *production rules* (henceforth simply called rules). Each rule has a condition part and an action part. An *inference engine* is used to execute the rules: it determines which rules are *enabled* by checking which conditions are true, and select rules to execute or *fire* from the enable rules. A rule is fired by executing its action part, which typically will have a side effect of performing input and output or computing a value and updating some date in memory. Thus rule-based systems fundamentally based on the notion of side effects.

Many mechanisms can be used to select which enabled rule to fire. In the literature, the term *conflict set* is often used to name the set of rules which are enabled, reflecting the intuition that somehow these rules can be conflicting, that is update the store in different incompatible ways. To ensure such problems do not occur, production systems will typically select a single rule to fire, by methods involving various notions such as priority or probabilities, along with notions such as the best matching of the conditions and so on.

In a rule-based system, control is *data-driven*, that is the data determines which part of the program will execute — furthermore, communication between different units is done solely through the use of data. There is no concept of a subroutine call to another rule, or anything of that sort. Rule-based systems allow a cleaner separation of knowledge (in the form of rules) from the control (encapsulated in the inference engine). This makes rule-based systems well-suited to program expert systems for analysis problems, and for programs for which the exact flow of control is not known. In this paper, we will make the further refined point that with the appropriate extensions, rule-based systems are well-suited for the compositional development of reactive systems.

## 3 A framework for rule-based programming

We begin by considering a general framework for the handling of simple OPS-style rules in SML, where actions are executed atomically and terminate before a new rule can be fired. We discuss the framework abstractly, that is in terms its interface.

The first notion of importance is that of a set of rules, as an abstract type with a single basic operation that creates an empty set of rules. The reason for keeping the set of rules abstract is to allow for different implementations, some possibly aimed at optimizing the evaluation of the conditions.

* type rule\_set

val newSet : unit -\> rule\_set

A rule is simply defined as a pair of a condition and an action, as in OPS. Since we want the condition to be dynamically evaluated, it is implemented as a function (the standard way to delay evaluation in an eager language such as SML). A condition evaluates to a positive integer (a word), which we call the *fitness* of the condition, indicating the degree to which the condition is satisfied. There is no *a priori* semantics or range associated with those, they are left to the discretion of the programmer. The only restriction is that a fitness of zero is used to indicate that a rule is not enabled. One can implement simple boolean conditions as values 0 and 1, if need be. The only operation on rules is to add them to rule sets. Note that because functions are first-class in SML, rules become first-class as well. That is, they can be passed as argument to functions, and returned from functions.

* type rule = {cond : unit -\> word,

action : unit -\> unit}

val addRule : rule * rule\_set -\> rule\_set

val mkSet : rule list -\> rule\_set

The final ingredient of the system is the inference engine, that selects which of the rules will be fired. At this point, the issue of *when* the rules should be fired must be addressed. Many systems tie the firing of the rules, or at least the evaluation of the conditions, to a change to variables that affect the conditions. We choose a much more fundamental approach using an explicit call that monitors the conditions, from which we can derive a *change-of-state* trigger. While inefficient, this approach has the advantage of being general. The other issue this does not address is that of conflict resolution, which rules to fire of all those that are enabled. We provide a selection of conflict resolution strategies. A function monitor is used to select rules to fire in a given set of rules. It takes as argument a conflict resolution flag determining the resolution strategy to use:

* datatype conflict\_res = AllBest

| RandBest

| AllDownTo of word

| RandDownTo of word

val monitor : conflict\_res -\> rule\_set -\> unit

AllBest fires all the rules that qualify as the best rule to apply, sorted according to fitness. RandBest randomly picks one of the rules that qualify as the best rule. AllDownTo and RandDownTo perform similarly, but consider all the rules whose fitness is at least the given value. One can easily extend the framework to allow for custom conflict resolution, which we do not pursue in this paper for simplicity.

Our notion of fitness is general. As we noted, we can imagine a binary fitness (0-1) for boolean firing, but also a fitness based on how close the conditions are to being completely satisfied, or even so far as how many conditions are actually satisfied (if we allow firing based on partially satisfied conditions). An easy extension to the framework would be to pass information from the condition to the action. We can mimic this easily by using a reference cell which can also be hidden in a closure of the rule, as follows:

* let

val r = ref 0

in

{cond = compute fitness, store something in r,

action = some action using the value in r}

end

In striking difference with other systems, rules are not persistent in this framework: once a rule fires and executes, it is removed from the rule set. To make a rule persistent, we can use the function persistent. This allows us to dispense with a function to remove rules from rule sets. Moreover, persistent comes for free given our extension for managing overlapping rules, as we will see in the next section.

* val persistent : rule -\> rule

As an example of how to use the framework, consider a simple rule-based program to compute the greatest common divisor of two integers, the classic Euclid’s algorithm. The example is artificial (it is easily implemented in SML without rules), but serves well to illustrate the basics. A more complete example is presented in Section 5. The program can be expressed as follows in Dijkstra’s language of guarded commands:

$$\text{do}X>YX:=XY|Y>XY:=YX\text{od}$$

The corresponding code in our framework is more verbose, but essentially similar:

* fun gcd (x,y) = let

val rx = ref x

val ry = ref y

fun r (r1,r2) = persistent

{cond = fn () =\> if (!r1) \> (!r2)

then 1 else 0,

action = fn () =\> r1 := (!r1) - (!r2)}

val rs = mkSet \[r (rx,ry),r (ry,rx)\]

fun loop () = if ((!rx) = (!ry)) then (!rx)

else (monitor AllBest rs; loop ())

in

loop ()

end

The above example is interesting because it shows how the fact that rules are first-class in the framework allow for parametrized rules: the rule r in the above code is parametrized over the reference cells containing the two arguments, parameterization which nicely showcases the symmetry of the rules.

As a final remark, we note that the rules as we have presented them in this section are simpler than they are in the actual framework. The rules in the implemented framework contain an extra field generically called *data* whose type is a parameter to the structure implementing the rules (in effect, the library is implemented as a functor). The monitoring function takes as an extra argument a function to compute a fitness both from the result of the evaluation of the condition and the data (which for example can contain notions such as rule priority and so on). Since describing this explicitly would require us to go into the details of both the module system and the type system of SML, we punt on these issues in this paper.

## 4 Managing overlapping rules

The main point of this work was to introduce *overlapping rules*, that is rules that can span multiple invocations and be performed in parallel with other rules. From an interface point of view, the framework only requires the addition of a single primitive, which we call wait, to add the desired functionality:

* val wait : (unit -\> word) option222a primitive SML type defined as datatype ’a option = NONE | SOME of ’a. -\> unit

The semantics of wait is simple. Fundamentally, wait interrupts the action of the current rule, as if the rule was finished executing, except that at the next time the monitoring function is invoked to fire a rule, the rules that were interrupted are allowed to resume while the new rule fires. Hence the term *overlapping*. In fact, the wait primitive takes two forms. The form wait (NONE) behaves as an unconditional interruption. Execution continues the next time the monitoring function is invoked. The form wait (SOME (f)) with f as a condition (that is, a unit $``$ word function) also interrupts, but only resumes the rule the next time the monitoring function is invoked with the condition f being satisfied. In effect, wait (SOME (f)) interrupts the rule and conceptually replaces it with a new rule containing the remainder of the interrupted rule, with a condition f.

As we mentioned in the previous section, rules are not persistent: once they are fired and execute, they are removed from the rule set.We can implement the persistent function using wait:

* fun persistent {cond,action} =

{cond = cond,

action = let

fun loop () = (action ();

wait (SOME (cond));

loop ())

in

loop

end}

This nicely shows the power of first-class rules.

## 5 An application: goal-directed agent control

We describe in this section one of the motivating applications for the development of the framework in the first place, that of controlling goal-directed agents. The architecture we have in mind is inspired by Hap, a reactive, goal-driven architecture for controlling agents in the Oz virtual environment . The main structure in Hap is an *active behavior tree* (ABT), which represents all the goals and behaviors an agent is pursuing at any given point . An agent chooses the next step to perform by selecting one of the leaves of its ABT. Three types of actions can be performed depending on the type of the node selected.

1. Primitive physical action: an action sent to the action server, which can either succeed or fail depending on the state of the world.
2. Primitive mental action: an action that simply performs a computation (possibly with side effects) and which always succeeds.
3. Subgoal: an action corresponding to a new subgoal; an appropriate behavior is selected that matches that subgoal, and the ABT is expanded by adding the steps specified by the behavior to the tree as children of the subgoal selected.

(For our purposes, we drop the distinction between mental and physical actions, since we can model mental actions as physical actions that always succeed). Programming an agent reduces to programming behaviors for various goals. Goals have no intrinsic meaning, they are simply names on top of which the programmer can attach any semantics she desires. Behaviors are defined by specifying to which goal they apply, a pre-condition for the application of that behavior (simply a predicate over the state of the world), and the steps that the behavior prescribes (physical actions, mental actions, subgoals). At this point, many details enter the description, to provide control over managing goals and behaviors. There are three kind of behaviors:

1. Sequential: the steps are performed in order, and failure of any step signifies the failure of the corresponding goal to which the behavior is attached; success of the last step signifies success of the corresponding goal.
2. Concurrent: the steps are performed in any order, but again failure of any step signifies the failure of the corresponding goal.
3. Collection: the steps are performed in any order, but success or failure of the steps are irrelevant. When all steps have succeeded or failed, the corresponding goal succeeds.

Subgoals steps in behaviors can moreover be annotated as *persistent*, that is when they succeed or fail, they are not removed, but rather persist as a continuing goal. Typically, top-level goals are persistent. Conversely, subgoals can be annotated with a *success test*, a predicate over the state of the world, which gets tested every time the ABT is activated. If the success test of a subgoal is true, the subgoal automatically succeeds.

Choosing a step to perform is by default done at random over all the applicable steps in an ABT, that is all the leaves that can be either executed right away or subgoals that can be expanded because a behavior applies (no behavior may apply because either none has been defined or no pre-condition is satisfied). Similarly, choosing a behavior to perform once a subgoal step has been chosen is by default done at random over all applicable behaviors. One can modify this default by assigning *priorities* to various subgoals.

All of this is meant to evoke the kind of structure we would like to express in our framework. Since Hap revolves around the notion of goals, we abstractly provide a notion of goal to the framework of the previous sections, where goals are for simplicity represented as strings.

* datatype goal\_status = Success | Failure

| Active | Available

| NoSuch

val goalSet : string -\> unit

val goalSucceed : string -\> unit

val goalFail : string -\> unit

val goalStatus : string -\> goal\_status

val goalClear : string -\> unit

where goalSet enables the given goal, such that it is to be pursued by the agent (it becomes available). The functions goalSucceed and goalFail are used to record that a goal has succeeded or failed. The function goalStatus returns the status of the given goal. The status of a goal is either Success or Failure if the goal has been recorded as such, or Active if a behavior is actively pursuing the goal, but is not done with it yet. A status of Available indicates that the goal is enabled, but that no behavior is pursuing it, while a status of NoSuch indicates that no such goal exists. The function goalClear removes a goal from the active list of goals. We define the boolean-valued helper functions isAvailable and isDone to check the status of a goal to be respectively Available or Success/Failure.

We interpret an ABT behavior as a rule, triggered both by the pre-condition of the behavior (if present) and the apparition as Available of the goal the behavior is meant to pursue. We do not worry about either sequential, collection or concurrent annotations, choosing rather to let the programmer manage the steps of the behavior explicitly. Patterns quickly emerge. For instance, a behavior for goal g triggered by a condition c and sequentially performing subgoal g1, action a and subgoal g2 can be interpreted as a rule:

* val beh1 =

{cond = fn () =\> if isAvailable (gandalso c

then 1 else 0,

action = fn () =\>

(goalSet (g1);

wait (SOME (fn () =\> isDone (g1)));

case (goalStatus (g1))

of Success =\> (goalClear (g1);

a;

wait (NONE);

goalSet (g2);

wait (SOME (fn () =\> isDone(g2)));

case (goalStatus (g2))

of Success =\> (goalClear (g2);

goalSucceed (g))

| _ =\> (goalClear (g2);

goalFail (g)))

| _ =\> (goalClear (g1);

goalFail (g)))}

Similarly, the previous behavior can be implemented concurrently by setting all the goals at once and waiting for all the goals to be done.

* val beh2 =

{cond = fn () =\> if isAvailable (gandalso c

then 1 else 0,

action = fn () =\>

(goalSet (g1);

goalSet (g2);

a;

wait (SOME (fn () =\> isDone (g1andalso isDone (g2)));

case (goalStatus (g1),goalStatus (g2))

of (Success,Success) =\> (goalClear (g1);

goalClear (g2);

goalSucceed (g))

| (\_,\_) =\> (goalClear (g1);

goalClear (g2);

goalFail (g)))}

By virtue of the andalso in the first waiting condition of beh2, the conditions must all be true for the system to proceed at that point. Although it does sequentializes the testing of the conditions, this is not an issue given our current framework since the rule is resumed when all the conditions are satisfied. It may however become an issue if we attempt to optimize the satisfaction of the conditions. Persistence of goals can be implemented by clearing them and setting them again, while goal success tests can be wrapped inside the wait condition for that goal.

One difficulty of this approach, immediately noticeable from the above code, is that it is very error-prone, even if it is much more flexible than ABTs. In effect, we have to implement the handling of the goals explicitly, for every single behavior. However, the flexibility of first-class rules and first-class functions allows us to easily generate such rules from a declarative description of the intended behavior. Consider a type behavior that describes a behavior declaratively:

* type behavior = {goal : string,

precond : (unit -\> bool) option,

kind : behavior\_kind,

steps : behavior\_step list}

datatype behavior\_kind = Sequential | Concurrent

datatype behavior\_step = Subgoal of string

| Action of unit -\> bool

(For simplicity, we drop the collective kind of behavior, its handling similar enough to the concurrent one to not cause a problem). We can describe the previous two behaviors beh1 and beh2 as follow (this time using behavior parameterization!):

* fun beh (k) = {goal = g,

precond = SOME (fn () =\> c),

kind = k,

steps = \[Subgoal (g1),

Action (fn () =\> a),

Subgoal (g2)\]}

val beh1 = beh (Sequential)

val beh2 = beh (Concurrent)

We can then interpret such descriptions in our framework, by a function behaviorRule, which takes a description of type behavior and returns a rule of type rule. The implementation of behaviorRule is simply a question of writing a rule whose action is an interpreter for lists of behavior\_step. For the truly interested, the code for behaviorRule is given in Figure 1.

## 6 Implementation

The framework we have described has been implemented for the Standard ML of New Jersey compiler  using the reactive library described in . One advantage of this approach is that the semantics of the system is easily derivable from the one in . Before discussing the implementation, let us give an overview of the reactive library.

The library defines a type rexp of reactive expressions, which are expressions that define control points. A reactive expression is created through a function rexp that expects a unit $``$ unit function as argument. The argument function calls the function stop to define a control point. The function react is used to take a reactive expression and to evaluate the code starting from the last control point reached until the next control point is reached. This is called *activating* a reactive expression. When a reactive expressions evaluates to a value without reaching a control point, it is said to *terminate*. Interesting combinators can be defined to take reactive expressions and combining them. The most important of such combinators is the merge combinator, which takes two reactive expression $`e_1`$ and $`e_2`$ and creates a new reactive expression $`e`$ that behaves as follows: when $`e`$ is activated, $`e_1`$ and $`e_2`$ are activated, one after the other. In effect, this interleaves the execution of $`e_1`$ and $`e_2`$. The combinator loop takes a reactive expressions $`e_1`$ and creates a new reactive expression $`e`$ that behaves as follows: when $`e`$ is activated, $`e_1`$ is activated. If $`e_1`$ terminates, it is reset (the reactive expression is re-instanciated) and activated again. The reactive expression nothing simply terminates immediately.

A more fine-grained notion of control point is also available. A reactive expression can call suspend to suspend its current execution. A suspend acts as a stop, except that a special combinator close is available. Given a reactive expression $`e_1`$, close returns a new reactive expression which behaves as follows: when $`e`$ is activated, it repeatedly activates $`e_1`$ until all its reactive subexpressions have reached stop-defined control points. This allows finer control over the order of execution of the reactive subexpressions, an example of which we will see in this section.

It is clear, given this description, how the library may be useful to implement the details of overlapping rules. For brevity, we assume the reactive library has been bound to a structure R. The implementation of the framework is rather simple, although it is complicated by technical details and some mismatches with the underlying reactive library. We define a rule set as a reactive expression, a *merge* of all the relevant rules.

* type rule\_set = R.rexp

fun newSet () = R.nothing

A rule is defined as before, while adding a rule to a set of rules consists of merging the reactive expression corresponding to the rule in the merge of the rule set.

* type rule = {cond : unit -\> word,

action : unit -\> unit}

fun addRule rs {cond,action} = let

val r = R.exp (fn () =\> (condition (cond);

action ()))

in

R.merge (rs,r)

end

val mkSet = foldr addRule (newSet ())

The function condition terminates immediately if monitoring determines that it should be fired (according to the fitness of the condition), otherwise it stops to wait for another instant where the condition is deemed fit. To bypass a limitation of the current reactive library, which is more geared towards locally determining whether a given reactive expression is allowed to continue rather than being selected through a global check, we introduce a global variable to hold a list of all computed fitnesses and allow the system to select the ones that will execute 333This does make the library non-reentrant. This could be corrected by an appropriate change to the reactive library (implementing reaction-specific data, for instance), or by including a notion of execution context to bind the use of the variable to a given context..

* val globalFitnesses = ref (\[\]) : (unit ref * word) list ref

We uniquely tag each condition being computed using a value of type unit ref, a trick commonly used in SML to get unique identifiers that can be quickly compared for identity. When the function condition is encountered, the current fitness is computed and stored along with a unique identifier for the condition. The reactive expression is then *suspended* (not stopped). This gives the other reactive expressions running *in parallel* a chance to evaluate their conditions. Upon resumption of the suspension, each condition checks if it is allowed to continue by seeing if it is listed in a list of *allowed-to-continue* conditions, stored in the previous globalFitnesses variable.

* fun condition (f) = let

val r = ref ()

fun loop () = (globalFitnesses := (r,f)::(!globalFitnesses)

suspend ();

if (List.exists (fn (r’,\_) =\> r=r’) (!globalFitnesses))

then ()

else (stop (); loop ()))

in

loop ()

end

The function wait that implements an interruption of the execution of the rule is simple:

* fun wait (NONE) = (stop (); suspend ())

| wait (SOME (f)) = (stop (); condition (f))

(The suspend in wait (NONE) is a technicality, needed to prevent the firing of stopped rules with no conditions until all the conditions of the other rules have been checked).

The core of the work is done in monitor, which is in charge of activating the reactive expression corresponding to the rule set, gather the results, compute which conditions are satisfied, and resume the suspensions.

* fun monitor c rs = let

val clearVar = R.loop (R.rexp (fn () =\> (globalFitnesses := \[\];

stop ())))

val r = R.loop (R.rexp (fn () =\> (computeEnabled (c);

stop ())))

in

R.react (R.close (R.merge (clearVar,R.merge (rs,r))))

end

The above code encompasses the control flow of the monitoring process, and relies heavily on the fact that merges are deterministic. A non-deterministic implementation could play with suspend calls to achieve the right order of execution: we need to make sure at every instant that clearVar is executed first, and computeEnabled is activated after all suspensions. The actual rule selection is performed by computeEnabled:

* fun computeEnabled (AllBest) = let

fun max (curr,\[\]) = curr

| max (curr,(\_,x)::xs) = if x\>curr then max (x,xs) else max (curr,xs)

val bestVal = max (1,!globalFitnesses)

val lst = List.filter (fn (\_,a) =\> a=bestVal) (!globalFitnesses)

in

globalFitnesses := lst

end

| computeEnabled (RandBest) = (as AllBest, but return random element)

| computeEnabled (AllDownTo (v)) = let

val lst = List.filter (fn (\_,a) =\> (a\>=v)) (!globalFitnesses)

in

globalFitnesses := lst

end

| computeEnable (RandDownTo (v)) = (as AllDownTo, but return random element)

## 7 Conclusion

So what have we done? We have developed a general framework for rule-based programming in SML, flexible enough to handle standard OPS-style rules, as well as overlapping rules. The fact that rules are first-class in the framework gave us enough flexibility to express Loyall’s active behavior trees through an interpreter of declarative behaviors.

We have not worried about the efficiency of the framework. Some rather standard optimizations as found in modern rule-based systems can easily be applied, although optimizations of the conditions may require them to be lifted into a more structured type, such as for example

* datatype condition = Basic of unit -\> word

| And of condition list

| Or of condition list

| Not of condition

| ...

in order to allow for such things as caching of condition evaluations and so on. On another note, evaluating a condition is only required if the references the condition refers to have been changed, so an optimized framework should maintain a list of the references used by conditions and record changes accordingly.

The tradeoff between generality and conciseness is not new. If we are willing to restrict flexibility, we can design an appropriate surface language which we can translate into this framework for execution. This is what we did for the implementation of behaviors in Section 5. This makes our framework a target language for the interpretation of domain-specific rule-based languages. In a similar way, the reactive library of  was designed as a target language for the interpretation of domain-specific reactive languages, which is the way it is being used in this paper.

Using the SugarCubes framework described in , we could move most of the framework to Java, but SML has distinct advantages, at least at first brush: we have seen how first-class rules and first-class functions helped design suitable domain-specific abstractions; the module system, although not used in this paper, plays a crucial role in the generalization of the framework to general rules that can carry arbitrary data (not only conditions and actions). The whole approach is also helped by the fact that SML already implements state-based programming through the use of explicit references. It will be interesting to see how much of this can be carried over to Java.

This frameworks, one of the first implemented using the library in , also points to some conclusions about the reactive approach. If the order in which parallel reactive expressions are activated is important, then we need to be careful, appropriately using suspend calls to get the order right; this makes the resulting system brittle and the code hard to see correct. A better approach may be to allow one to explicitly control the order of execution of the branches of a merge. On a related note, the reactive library is geared towards determining reactive expression activation locally, while we have encountered in this paper a reasonable instance where the activation decision is taken on a global level. It would be interesting to find an extension of the reactive library to do that cleanly, without resorting to a list of unique identifiers indicating which expression is allowed to resume.


text:
# Amino acid classes and the protein folding problem

The principal theme of this paper is to address the issue of determining the minimum number of distinct amino acids that are needed to make protein-like sequences with folds similar to those found in nature. A general answer to this question would have important ramifications in the design of proteins and in the origin of life. Our analysis is carried out within the framework of the Miyazawa-Jernigan (MJ)  interaction matrix which was derived by them using a statistical approach and is commonly used as a measure of the coarse-grained interactions between amino acids in a protein.

Experiments  and theoretical studies  have shown the critical role played by two kinds of residues: hydrophobic and polar. There is an effective attraction between hydrophobic amino acids that arises from their aversion to the solvent and lead to such amino acids forming the core in the protein native state. Recent protein engineering experiments suggest that maybe not two but certainly several amino acids can effectively substitute for the full 20-alphabet set. It has been shown  that helical bundles can be built with the set of three amino acids: hydrophobic L (leucine), polar E (glutamic acid), and polar K (lysine). Helical bundles have also been found  to be viable when hydrophobic sites are filled from the set \[F (phenylalanine), L, I (isoleucine), M (methionine), V (valine)\] and the polar sites from the set \[E, D (aspartic acid), K, N (aspargine), Q (glutamine), H (histidine)\]. Encoding the $`\beta `$-sheet SH3 domain, however, requires five amino acids: hydrophobic I (isoleucine), K, E, A (alanine), and G (glycine) . This suggests that the 20 amino acids can be grouped into five distinct clusters with the members of each group having quite similar properties. Riddle et al.  and Wolynes  have presented persuasive arguments that five groups are needed to provide enough specificity to form a folding funnel and generate few traps in the energy landscape.

Wang and Wang  have suggested that a justification for the five group clustering scheme can be provided by minimizing the mismatch between the reduced and complete interaction matrices. Specifically, they considered the MJ  matrix and deduced a clustering scheme, shown at the top of Figure 1, in which the representative amino acids were IKEAG – precisely as in the experiments on SH3 . The computational scheme for the clustering is stochastic in nature and does not permit a single entry group due to technical reasons associated with the computation.

Here, we implement a much simpler and deterministic clustering scheme which is based on considering the ”distances” between the amino acids. The groupings we get are different from those obtained by Wang and Wang . We have carried out detailed tests of two characteristics of real proteins – thermodynamic stability and kinetic accessibility – on well-designed sequences within the context of a lattice model in three dimensions. When the individual amino acids are substituted by the representatives of their groups, the thermodynamic stability test is passed successfully both by the Wang-Wang and our schemes. However, there is a qualitative difference in the test results on rapid folding into the native state, with the Wang-Wang scheme performing significantly worse than our physically-based approach. A straightforward eigenvalue analysis of the kind carried out by Li et al.  and Chan  lead to results in perfect accord with our clustering scheme. Furthermore, when one works with the bare MJ matrix  (without subtracting off the mean value as done by Li et al. ), one obtains the simple result that the MJ matrix may be represented as

$$M_{ij}^1=\lambda _1v_iv_j,$$
(1)
where $`\lambda _1`$=69.68 and $`v`$ has rank-ordered components 0.333, 0.332, 0.309, 0.282, 0.276, 0.275, 0.255, 0.249, 0.198, 0.196, 0.177, 0.171, 0.169, 0.163, 0.161, 0.158, 0.152, 0.144, 0.144, 0.125 in the order LFIMVWCYHATGPRQSNEDK respectively. The representation given in eq. 1 is similar in spirit to the reprezentation  in terms of “charges”, $`q_i`$, which yields $`M_{ij}^{}=C_0+C_1(q_i+q_j)+C_2q_iq_j`$. Eq. (1) results in the same level of accuracy (an average error of about 6%), even though it involves only one constant, $`\lambda _1`$, instead of three, $`C_0`$, $`C_1`$, and $`C_2`$.

The interactions of a given amino acid $`i`$ with each of the 20 amino acids form a vector with 20 components. There are 20 such vectors and the Euclidean distance, $`R_{ij}`$, between amino acids $`i`$ and $`j`$ is defined through $`R_{ij}^2=_k(M_{ik}M_{jk})^2`$. Similarly, the Manhattan distance involves the sum of the absolute values of $`M_{ik}M_{jk}`$. $`R_{ij}`$ is a measure of the fidelity of substitution of one amino acid by the other. As in the construction of optimal paths in a strongly disordered medium , we select two amino acids that are separated by the shortest distance and combine them into one group. The effective couplings of a group with other groups or individual amino acids or indeed with itself is simply obtained as an arithmetic average over the individual amino acid interactions. More generally, a weighted average would be appropriate if the weights, which could depend on the frequency of occurrence of a given amino acid and other factors, were known. The procedure is now iterated resulting in fewer and fewer groups. The leader of the group is determined as the amino acid whose original couplings deviate the least from the effective couplings associated with the group. For a two-member group the choice of a leader is ambiguous. At each stage, the advancement of clustering is characterized by measuring the smallest distance between the remaining individuals or groups, $`R_{min}`$. Figure 2 shows the behavior $`R_{min}`$ as a function of $`N_a`$, the number of groups at any stage of the iteration process. For either distance metric, the number $`N_c=5`$ is a special number of groups, below which $`R_{min}`$ goes up strongly. This sharp increase is a direct reflection of the forced clustering of incompatible amino acids.

Our analysis can also be used to probe the nature of other theoretically derived potentials of interaction between amino acids. For example, for the matrix of interactions proposed by Kolinski et al. , the low $`N_a`$ growth is much weaker but $`N_c`$ is still around 5. For the matrix introduced by Betancourt and Thirumalai , $`R_{min}`$ varies very weakly with $`N_a`$ and the sharp increase takes place at $`N_c`$ equal to 3.

Our results for clustering of the Miyazawa-Jernigan amino acids are shown in Figure 1. Our clustering into five groups separates the single hydrophobic group of Wang and Wang into two groups and it isolates K as a single element group. The membership of the five groups does not depend on the distance metric, even though different amino acids are selected as “leaders”, or the most characteristic representatives of the groups. The Euclidean choice proposes FVASK as the best whereas the choice of the Manhattan metric yields LWASK. These groupings are consistent with the experimental results except that Riddle et al.  take two choices, E and G, from the TGPRQSNED group and no choice from the second hydrophobic group MVWCY. It is heartening that the ultimate division into hydrophobic and polar groups is quite robust and yields the same group members in both schemes and for any distance metric. In order to test the five group clustering, we have considered a 27-monomer lattice model with the Miyazawa-Jernigan couplings. The native states are maximally compact, i.e. they fit the $`3\times 3\times 3`$ lattice. We generate a bank of 94 well designed sequences. The first step in the design procedure  involved the generation of a random sequence with a uniform composition. The second step was to perform random permutations in the placement of the amonoacids and to select only those sequences which have energy gaps (the energy difference between the two lowest energy maximally compact conformations) larger than 6.1 (in units of the MJ interaction matrix). This bank of sequences was then tested against substitution of the amino acids by the leaders corresponding to the five-group division. For both the distance-based and mismatch based schemes we find identical performance – 95% of the resulting effective sequences continue to have nondegenerate native states which coincide with the native state conformations of the original sequences. (A lower performance rate of 80%, quoted by Wang and Wang , is obtained on considering less stringently designed sequences.)

A more stringent test is provided by the kinetics. Figure 3 shows the folding behavior of four of the well-designed sequences. The first passage time in 21 runs starting from randomly chosen unfolded conformations are shown for all cases which have a folding time less than $`5\times 10^8\tau _0`$. The dynamics are based on a Monte Carlo process which satisfies detailed balance . The folding is studied at the folding transition temperature, $`T_f`$, at which the probability of being in the native conformation equals $`\frac{1}{2}`$. This temperature was located using a long unfolding Monte Carlo process, which explores the entire space of conformations and is not restricted merely to maximally compact conformations. The distance based substitution either raises $`T_f`$ or leaves it intact and preserves the range of values of the folding times. On the other hand, the mismatch based substitution lowers all four values of $`T_f`$ (some of them significantly) and extends folding times substantially. The simple origin of this failure may be traced to the fact that $`K`$ which is a relatively inert amino acid is chosen as one of the leaders in the mismatch scheme and sequences with many amino acids belonging to the group with $`K`$ are all substituted by an innocuous representative.

The clustering schemes illustrated in Figure 1 correspond to the partitioning of a string of all amino acids into segments provided that the amino acids are first arranged into a particular order. The order corresponding to the mismatch based scheme is almost the same as shown in Figure 1 for the distance based scheme except that K is forced to be placed earlier along the string. What is it that determines this optimal order? Following Li et al.  and Chan , we address this question through an eigenanalysis of the MJ matrix. As shown by Chan, the MJ matrix is clearly dominated by one mode and lends itself to an exceedingly simple representation. Figure 4 shows that this one mode description of the MJ matrix is fairly accurate (the largest errors are located at the smallest magnitude entries). In fact, it is as accurate as the two mode description of Li et al. . Their two mode picture emerges because of the subtraction of the mean value from each matrix element. This mean value has physical significance because it determines the degree of compactness and indeed the degree of aggregation of a protein . It is interesting to note that the dominant eigenvector varies across the amino acids in a smooth way except near the transition from Y to H which marks a transition from the hydrophobic to polar groups. The order in which the amino acids are plotted on the x axis in Figure 4 and in (most of) Figure 1 is thus determined by the weight with which a particular amino acid contributes to the dominant eigenvector. In summary, we have shown that a simple distance-based scheme of clustering leads to a powerful and simple representation of groupings of amino acids making up a protein. Detailed tests within the context of a MJ matrix and a lattice model show that both the thermodynamic stability and the folding kinetics of protein-like sequences are preserved by the substitution of the full 20 amino acid alphabet by merely 5 groups. Our scheme is found to be in good physical accord with that obtained using eigenanalysis and does not have the defects associated with a more complex mismatch based scheme.

We are grateful to Ruxandra Dima for fruitful discussions. This work was supported by KBN (Poland - grant 2P03B-146-18), MURST (Italy), NASA, IGERT-NSF and the Petroleum Research Fund administered by the American Chemical Society.


text:
# Observations of Mkn 421 with the CELESTE Experiment

## Introduction

CELESTE has been in stable operation with a 40 heliostat array since October 1999. Our main targets during this initial period have been the three most reliably detected sources in the TeV sky; the Crab nebula, Markarian 421 and Markarian 501. Our detection of the Crab nebula is described elsewhere in these proceedings crabceleste . This report is concerned with the observations of Mkn 421.

Mkn 421 is the closest X-ray selected BL Lac object (z=0.030) and was the first extragalactic object detected at TeV energies punch . The emission from this source is extremely variable over timescales as short as 15 minutes gaidos  and multiwavelength observations have shown that the TeV variability is correlated with the X-ray emission buckley . Observations in the 50 GeV region, not possible until now, will help to improve our understanding of the emission processes at work.

The CELESTE experiment is described in detail in the experimental proposal proposal  and in NIM ; berry . Briefly, CELESTE uses forty, $`54m^2`$ heliostats of a former solar electrical plant at the Themis site in the French Pyrénées (N. $`42.50^{}`$, E. $`1.97^{}`$, altitude 1650 m). The light from the seperate heliostats is resolved using secondary optics at the top of a 100 m tall tower. At the secondary mirror focus is one PMT for each heliostat. A solid Winston cone at the face of each PMT defines the field of view of each tube to be 10 mrad.

The trigger is constructed by summing eight PMT signals in five groups and sending each sum to a discriminator with a threshold set to 4.5 photo-electrons (PE) per heliostat. A logic coincidence of at least three of the five groups triggers the experiment. Each PMT signal is sent to a 1 GHz flash ADC circuit which digitizes the signal. When a trigger occurs, digitization stops and a window of 100 ns centered at the nominal Čerenkov pulse arrival time is read out.

## Analysis

CELESTE, with its large mirror area per PMT, is very susceptible to systematic effects in the data caused by differences in the night sky background light between the ON source observation and OFF source control data. In the case of Mkn 421, a star of magnitude 6.1 located at the source position causes a difference of $``$15% in the PMT anode currents between the two fields. To compensate, we have adapted the software padding method proposed by Cawley cawley  for use with our FADC data crab . We also apply a software trigger by using the FADC data to reconstruct the five trigger group sums and then demanding:$``$4 groups $`>`$5.0 PE per heliostat. This increases the energy threshold, but provides us with comparable background data in the ON and OFF source fields. We also require $`>10`$ Čerenkov peaks with an amplitude $`>25`$ digital counts to ensure the shower parameters are measurable.

The rest of the analysis is concerned with improving the gamma ray signal to hadron background ratio using two parameters. The first, $`\sigma _{grp}`$, is a measurement of the homogeneity of the light distribution between the trigger groups and should be the most efficient background rejection cut according to our simulations, rejecting 80% of the remaining hadrons for only 30% of the gammas. The second parameter, $`\theta `$, is the shower axis angle relative to the pointing direction reconstructed using two points; the barycentre of the light distribution at ground level and the position of the shower core relative to the pointing direction at the altitude of maximum emission deduced from a sphere fit to the Čerenkov pulse arrival times at the heliostats.

## Observations and Results

Mkn 421 was observed between November 1999 and March 2000. After rejecting those observations with bad weather or equipment problems, the dataset consists of 80 ON/OFF pairs with an ON source exposure of 22.5 hours.

The CAT imaging telescope barrau , situated on the same site as CELESTE, made observations of Mkn 421 during the same period. In December 1999, no evidence for a signal was reported by CAT. The left-hand plot of Fig. 1 shows the distribution of the difference in the shower orientation parameter, $`\theta `$, between the ON and OFF source observations as measured by CELESTE during December. No evidence for a signal at low $`\theta `$ is seen (for 1.5 hours of ON source exposure). Table 1 shows the number of events remaining at each stage of the analysis.

In January/February 2000 CAT reported that Mkn 421 was in an active state. In particular, two bright flares were observed, with emission levels $`>`$2 Crab. The right-hand plot of Fig. 1 shows the difference distribution in $`\theta `$ for CELESTE observations during these flares (5.2 hours of ON source exposure). A strong signal is apparent at low values of $`\theta `$, as expected for gamma rays. The significance of the excess is 8.1 $`\sigma `$. Table 2 shows the improvement in the signal/background ratio at each stage of the analysis after the initial cleaning cuts. The significance before cuts is not very meaningful as the excess, while undoubtably largely due to gamma rays, is contaminated by systematic effects due to sky noise differences.

The whole Mkn 421 dataset has been analysed and Fig. 2 shows the lightcurve as observed by both CAT and CELESTE on a night by night basis. Emission levels up to $`3`$ Crab ($`7\gamma \mathrm{min}^1`$) were observed, detectable by CELESTE in a single night. Also shown in the figure is a plot of the gamma ray flux measured by CELESTE against that measured by CAT showing reasonably good correlation.

## Conclusions

The TeV blazar Mkn 421 has been detected by CELESTE for the first time at an energy of $`50`$ GeV. Evidence for night by night variablility has been observed, and the variablity measurements by CELESTE are correlated with those of CAT. Importantly, a non-detection of the source in a quiescent state, followed by a very significant detection in a flaring state demonstrates that the CELESTE analysis is working, and that we have been able to correct for systematic effects due to night sky light differences between the ON and OFF source fields.


text:
# Self-Gravitating N-Body Systems out of Equilibrium

## 1. Introduction

An energy-flow through a system is related to an entropy-flow. If the entropy-flow leaving the system is larger than the entering one, the system evacuates its internal, by irreversible processes produced entropy to the outer world. Consequently order is created inside the system, diverging from the classical thermodynamical equilibrium of a closed system. Far from equilibrium the system is no longer characterized by an extremum principle, thus losing its stability. Therefore perturbations can lead to long range order, through which the system acts as a whole. Such a behavior is well known in laboratory hydrodynamics and chemistry. The underlying concepts such as “dissipative structures” and “self-organization” were extensively studied (see e.g. Nicolis & Prigogine 1977). But in spite of their popularity they are until now only little studied in the context of non-equilibrium structures in self-gravitating astrophysical systems. Therefore we take up some ideas of these concepts and build a simple model of an open self-gravitating system. With this model we want to check if an energy-flow can maintain a self-gravitating system in an statistically stable state, out of thermodynamical equilibrium and if gravitation in combination with an energy flow can create structures with a higher degree of order.

## 2. Dissipative N-body Systems

Taking into account the highly clumpy nature of the interstellar medium, we use in our model dissipative particles, representing dense cloud fragments, to simulate cosmic gas (Pfenniger 1998). Moreover, with such a realization we can check some thermodynamic results of self-gravitating systems.

## 3. Model

In order to prevent gravitationally unbound particles from dissolution, we confine the particles in a spherical potential well. This prevents matter flow. But the system is subjected to an energy flow, maintaining the system out of equilibrium. This flow is sustained by energy injection (heating) on large scales and local dissipation. The energy injection is due to time and position dependent potential perturbations. If the system represents a molecular cloud, then these potential perturbations can stem from star clusters, clouds or other high mass objects passing in the vicinity. Indeed such stochastic encounters must be quite frequent in galactic discs and we assume, $`1/f_{enc}t_{dyn}`$, where $`f_{enc}`$ is the mean frequency of the encounters. Thus these encounters can provide a continuing low frequency energy injection on large scales. The dissipation is due to “inelastic particle encounters”. Therefore we add friction forces to the equation of motion, depending on the relative particle velocities and positions.

If we switch off the heating process and use instead the local dissipation scheme (“inelastic particle encounters”) a scheme that dissipates the energy globally, then we can maintain the system during its evolution nearly in thermodynamic equilibrium. Thus we can perform some simple “thermodynamic experiments” of self-gravitating systems with softened potentials.

## 4. Results

Potential perturbations caused by astrophysical objects passing in the vicinity of a system can compensate the energy loss due to dissipation, thus prevent the system from collapsing and maintain a statistical equilibrium. This statistically steady state consists of a cold core moving in a hotter halo. The energy-flow leads however not to a phase transition to higher ordered structures.

The gravo-thermal experiments lead to the following results: A Plummer softened potential yields the same range of negative specific heat as found by Follana & Laliena. The negative specific heat range is related to a phase transition, separating a high energy homogeneous phase from a collapsed phase. As long as cooling processes are at work the collapsed phase shows no core-halo structure. After the cooling has stopped such a structure is formed on the relaxation time scale.

## References

Follana E., Laliena V., cond-mat/9911107

Nicolis G., Prigogine I., 1977, Self-organization in non-equilibrium systems, Wiley, New York

Pfenniger D., 1998, in: $`H_2`$ in the Early Universe, eds. Palla F., Corbelli E., Galli D., Memorie Della Societa Astronomica Italiana, 429


text:
# Kinetic equation with exact charge conservation

## I Introduction

Relativistic statistical thermodynamics has long been used as a tool to describe particle production in heavy-ion and in high-energy particle collisions . Recent analyses have shown that the statistical models can indeed give a satisfactory description of the multiplicities of most hadrons measured in A-A collisions at AGS and SPS energies. However, the dynamics of particle equilibration, and in particular chemical equilibration, is still not well understood . In this Letter, we shall address the problem of chemical equilibration within the statistical kinetic approach.

Within this approach, particle production is commonly described using the grand canonical ensemble, where event-averaged multiplicities are controlled by chemical potentials. In this description the net value of a given $`U(1)`$ charge (e.g., electric, baryon, strangeness, charm, etc.) fluctuates. These fluctuations can be neglected only if the particles carrying the charge in question are abundant. In this case the charge will be conserved on the average. In the opposite limit of rare particle production, conservation laws must be implemented locally on an event-by-event basis ; i.e.,a canonical ensemble must be used.

The local conservation of quantum numbers in the canonical approach severely reduces the phase space available for particle production . Recently, it has been shown that the canonical statistical model provides a good description of particle yields measured in low-energy heavy-ion  and high-energy hadron-nucleus, hadron-hadron and $`e^+e^{}`$ reactions .

In view of these results it is of particular importance to formulate a kinetic theory for the time evolution of particle production in order to investigate the approach to the canonical chemical equilibrium. Obviously, if particles are abundantly produced, the equilibrium result for the particle multiplicity should coincide with the grand canonical one. For rare processes, however, the particle production is strongly correlated and the canonical equilibrium result is expected. The kinetic formulation for the production of strongly correlated particles was studied in the literature . However, no complete solution has been obtained.

In this letter we consider the time evolution of the multiplicity of particles that carry the charges corresponding to an $`U(1)`$ internal symmetry. We formulate a kinetic master equation valid for arbitrary particle multiplicity. It reproduces the canonical equilibrium solution for rare particle production, and reduces to the standard grand canonical rate equation for abundant particle production.

## II Rate equation

In the standard formulation , the rate equation for a binary process $`a_1a_2b_1b_2`$ with $`ab`$ is described by the following population equation:

$`{\displaystyle \frac{dN_{b_1}}{d\tau }}={\displaystyle \frac{G}{V}}N_{a_1}N_{a_2}{\displaystyle \frac{L}{V}}N_{b_1}N_{b_2},`$ (1)
where $`G\sigma _Gv`$ and $`L\sigma _Lv`$ give the momentum-averaged cross sections for the gain process $`a_1a_2b_1b_2`$ and the loss process $`b_1b_2a_1a_2`$, respectively. $`N_k`$ represents the total number of particles $`k`$, and $`V`$ is the proper volume. Among such processes, a typical example is the kaon production/annihilation via $`\pi ^+\pi ^{}K^+K^{}`$. The above rate equation, however, cannot be applied to the situation where particle production is rare and is strongly correlated by exact charge conservation.

To account for the correlation between the production/annihilation of particles $`b_1`$ and $`b_2`$, let us define $`P_{i,j}`$ as the probability to find a number $`i`$ of particle $`b_1`$ and a number $`j`$ of particle $`b_2`$ in an event. We also denote as $`P_i`$ the probability to find a number $`i`$ of particle $`b`$ in an event. The average number of particles $`b`$ per event is then defined as:

$`N_b={\displaystyle \underset{i=0}{\overset{\mathrm{}}{}}}iP_i.`$ (2)
We can now write the following general rate equation for the average particle multiplicity:

$`{\displaystyle \frac{dN_{b_1}}{d\tau }}={\displaystyle \frac{G}{V}}N_{a_1}N_{a_2}{\displaystyle \frac{L}{V}}{\displaystyle \underset{i,j}{}}ijP_{i,j}.`$ (3)
Further, let particles $`b_1`$ and $`b_2`$ carry opposite units of a charge, corresponding to an $`U(1)`$ internal symmetry (strangeness in the case of kaons). Then the $`U(1)`$ charge neutrality of the system gives $`NN_{b_1}=N_{b_2}`$. We have then,

$`P_{i,j}`$ $`=`$ $`P_i\delta _{ij},`$ (4)
$`{\displaystyle \underset{i,j}{}}ijP_{i,j}`$ $`=`$ $`{\displaystyle \underset{i}{}}i^2P_iN^2=N^2+\delta N^2,`$ (5)
where $`\delta N^2`$ represents the event-by-event fluctuation of the number of $`b_1b_2`$ pairs. Note that we always consider $`a_1`$ and $`a_2`$ particles abundant (such as, e.g., pions) so that we can neglect the event-by-event fluctuations of their multiplicity and the change of their number due to the considered processes.

Following Eqs.(3) and (5), the general rate equation for the average number of $`b_1b_2`$ pairs can be written as

$`{\displaystyle \frac{dN}{d\tau }}={\displaystyle \frac{G}{V}}N_{a_1}N_{a_2}{\displaystyle \frac{L}{V}}N^2.`$ (6)
For abundant production of $`b_1b_2`$ pairs, where $`N1`$,

$`N^2N^2,`$ (7)
and Eq.(6) obviously reduces to the standard form, i.e.,

$`{\displaystyle \frac{dN}{d\tau }}{\displaystyle \frac{G}{V}}N_{a_1}N_{a_2}{\displaystyle \frac{L}{V}}N^2.`$ (8)
However, for rare production of $`b_1b_2`$ pairs, where $`N1`$, the rate equations (1) and (8) are no longer valid. We have instead

$`N^2N,`$ (9)
which reduces Eq.(6) to the following form:

$`{\displaystyle \frac{dN}{d\tau }}{\displaystyle \frac{G}{V}}N_{a_1}N_{a_2}{\displaystyle \frac{L}{V}}N.`$ (10)
Thus, in the limit $`N1`$, the absorption term depends on the pair number only linearly, instead of quadratically for the limit $`N1`$.

## III Equilibrium multiplicity and relaxation time

To illustrate the differences in the time evolution of particle abundance, we consider the two limiting cases, $`N1`$ and $`N1`$, and present their equilibrium values and relaxation times for the production of $`b_1b_2`$ pairs. As an example, we consider a system at fixed temperature and volume and with no initial $`b_1b_2`$ pairs, i.e., $`N(\tau =0)=0`$.

In the limit when $`N1`$, the standard Eq.(8) is valid and has the following well-known solution:

$`N^{\mathrm{GC}}(\tau )=N_{\mathrm{eq}}^{\mathrm{GC}}\mathrm{tanh}\left(\tau /\tau _0^{\mathrm{GC}}\right),`$ (11)
where the equilibrium value for the number of $`b_1b_2`$ pairs $`N_{\mathrm{eq}}^{\mathrm{GC}}`$ and the relaxation time constant $`\tau _0^{\mathrm{GC}}`$ are given by

$`N_{\mathrm{eq}}^{\mathrm{GC}}=\sqrt{ϵ},\tau _0^{\mathrm{GC}}={\displaystyle \frac{V}{L\sqrt{ϵ}}},`$ (12)
respectively, with $`ϵGN_{a_1}N_{a_2}/L`$.

In the particular case where particle momentum distributions are thermal, the gain ($`G`$) and loss ($`L`$) terms just represent the thermal averages of the production and absorption cross sections with

$`{\displaystyle \frac{G}{L}}={\displaystyle \frac{d_{b_1}\alpha _{b_1}^2K_2(\alpha _{b_1})d_{b_2}\alpha _{b_2}^2K_2(\alpha _{b_2})}{d_{a_1}\alpha _{a_1}^2K_2(\alpha _{a_1})d_{a_2}\alpha _{a_2}^2K_2(\alpha _{a_2})}},`$ (13)
where $`d_k`$’s denote the degeneracy factors, and $`\alpha _km_k/T`$. The equilibrium value for the number of $`b_1b_2`$ pairs in Eq.(12) now reads as

$`N_{\mathrm{eq}}^{\mathrm{GC}}={\displaystyle \frac{d_{b_1}}{2\pi ^2}}VT^3\alpha _{b_1}^2K_2(\alpha _{b_1}).`$ (14)
Thus it is described by the Grand Canonical (GC) result with vanishing chemical potential due to our requirement of the (average) $`U(1)`$ charge neutrality of the system.

In the opposite limit where $`N1`$, the time evolution is described by Eq.(10), which has the following solution:

$`N^\mathrm{C}(\tau )=N_{\mathrm{eq}}^\mathrm{C}\left(1e^{\tau /\tau _0^\mathrm{C}}\right),`$ (15)
with the equilibrium value and relaxation time given by

$`N_{\mathrm{eq}}^\mathrm{C}=ϵ,\tau _0^\mathrm{C}={\displaystyle \frac{V}{L}}.`$ (16)
With a thermal momentum distribution the equilibrium value of $`b_1b_2`$ pair multiplicity has the following form:

$`N_{\mathrm{eq}}^\mathrm{C}=\left[{\displaystyle \frac{d_{b_1}}{2\pi ^2}}VT^3\alpha _{b_1}^2K_2(\alpha _{b_1})\right]\left[{\displaystyle \frac{d_{b_2}}{2\pi ^2}}VT^3\alpha _{b_2}^2K_2(\alpha _{b_2})\right].`$ (17)
This equation demonstrates the locality of the $`U(1)`$ charge conservation. With each particle $`b_1`$, a particle $`b_2`$ with the opposite charge is produced in the same event in order to conserve charge locally. This is the result expected from the Canonical (C) formulation of conservation laws .

We note that Eq.(17) is just the leading term in the expansion of the canonical result for the multiplicity of particles carrying the $`U(1)`$ charges. The general expression is known to have the following form :

$`N_{\mathrm{eq}}^\mathrm{C}=N_{\mathrm{eq}}^{\mathrm{GC}}{\displaystyle \frac{I_1(2N_{\mathrm{eq}}^{\mathrm{GC}})}{I_0(2N_{\mathrm{eq}}^{\mathrm{GC}})}},`$ (18)
where $`N_{\mathrm{eq}}^{\mathrm{GC}}`$ is given by Eq.(14) and $`I_i`$’s are the modified Bessel functions.

Comparing Eq.(12) and Eq.(16), we first find that, for $`N1`$, the equilibrium multiplicity in the canonical formulation is much lower than what is expected from the grand canonical result,

$`N_{\mathrm{eq}}^\mathrm{C}=(N_{\mathrm{eq}}^{\mathrm{GC}})^2N_{\mathrm{eq}}^{\mathrm{GC}}.`$ (19)
This shows the importance of the canonical description of charge conservation when the multiplicity of charged particles is small. We also note that the volume dependence in the two cases is different. The particle density in the GC (abundant) limit is independent of $`V`$, whereas in the opposite canonical (rare) limit the density scales linearly with $`V`$.

Secondly, the relaxation time for a canonical system is far shorter than what is expected from the grand canonical result,

$`\tau _0^\mathrm{C}=\tau _0^{\mathrm{GC}}N_{\mathrm{eq}}^{\mathrm{GC}}\tau _0^{\mathrm{GC}},`$ (20)
due to small number of particles ($`N_{\mathrm{eq}}^{\mathrm{GC}}1`$). For example, the total number of produced kaons in $`Au+Au`$ collisions at 1 GeV/$`A`$ is of the order of 0.02. Thus the canonical relaxation time is a factor of 50 shorter than what is expected from the grand canonical formulation.

We note from Eq.(6) that these two limits are essentially determined by the size of $`\delta N^2`$, the event-by-event fluctuation of the number of $`b_1b_2`$ pairs. The grand canonical results correspond to small fluctuations, i.e., $`\delta N^2/N^2<<1`$, while the canonical description is necessary in the opposite limit.

## IV Master equation

In this section, we formulate the general evolution equation which is valid for an arbitrary value of $`N`$. It is a master equation for $`P_n(\tau )`$, the probability of finding $`n`$ pairs of $`b_1b_2`$ at time $`\tau `$. This probability increases with time due to transitions from $`n1`$ and $`n+1`$ states to the $`n`$ state, while it also decreases due to transitions from the $`n`$ state to $`n1`$ and $`n+1`$ states. The transition probability $`nn+1`$ per unit time due to pair creation is $`GN_{a_1}N_{a_2}/V`$ and the transition probability $`nn1`$ due to pair annihilation is $`n^2L/V`$. Therefore, the master equation set has the form:

$`{\displaystyle \frac{dP_n}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{G}{V}}N_{a_1}N_{a_2}\left[P_{n1}P_n\right]`$ (21)
$``$ $`{\displaystyle \frac{L}{V}}\left[n^2P_n(n+1)^2P_{n+1}\right],`$ (22)
where $`n=0,1,2,3,\mathrm{}`$. Multiplying the above equation by $`n`$ and summing over $`n`$, one recovers Eq.(6), the general rate equation for the time evolution of the average number of $`b_1b_2`$ pairs. However, the master equation (22) contains much more information than the rate equation (6). It contains enough information to solve for the evolution of $`N`$ (and all moments of $`N`$) for arbitrary $`N`$. For example, one can also obtain the time evolution of particle fluctuations, $`\delta N^2`$, which are of physical importance for rarely produced particles.

We can convert the iterative equations (22) for $`P_n`$’s into a partial differential equation for the generating function

$$g(x,\tau )=\underset{n=0}{\overset{\mathrm{}}{}}x^nP_n(\tau ).$$
(23)
Multiplying Eq.(22) by $`x^n`$ and summing over $`n`$, we find

$`{\displaystyle \frac{g(x,\tau )}{\tau }}={\displaystyle \frac{L}{V}}(1x)\left(xg^{\prime \prime }+g^{}ϵg\right),`$ (24)
where $`g^{}g/x`$. Note that the $`g(1,\tau )`$ does not change with time, which is equivalent to the conservation of total probability evident in Eq.(22).

The equilibrium solution, $`g_{\mathrm{eq}}(x)`$, thus, obeys the following equation:

$`xg_{\mathrm{eq}}^{\prime \prime }+g_{\mathrm{eq}}^{}ϵg_{\mathrm{eq}}=0.`$ (25)
By variable substitution ($`x=y^2/(4ϵ)`$) this equation can be reduced to the Bessel equation. The solution that is regular at $`x=0`$ (since $`g(0)=P_01`$) is given by

$`g_{\mathrm{eq}}(x)={\displaystyle \frac{1}{I_0(2\sqrt{ϵ})}}I_0(2\sqrt{ϵx}),`$ (26)
where the normalization is fixed by $`g(1)=P_n=1`$.

The equilibrium probability distribution $`P_n`$ can now be found from Eqs.(23,26) as

$`P_{n,\mathrm{eq}}={\displaystyle \frac{ϵ^n}{I_0(2\sqrt{ϵ})(n!)^2}},`$ (27)
and the value of the average number of $`b_1b_2`$ pairs per event at equilibrium is given by

$`N_{\mathrm{eq}}=g^{}(1)=\sqrt{ϵ}{\displaystyle \frac{I_1(2\sqrt{ϵ})}{I_0(2\sqrt{ϵ})}},`$ (28)
which obviously coincides with the expected result for particle multiplicity in the canonical ensemble given by Eq.(18) and also reduces to the results for the two limiting cases in Sec.III.

## V Conclusions

We have formulated the kinetic master equation for strongly correlated production of particles, where the correlation is due to the local charge conservation required by an $`U(1)`$ internal symmetry. Our general rate equation is valid for arbitrary value of $`N`$, thus it reduces to the grand canonical results for large $`N`$ and to the canonical results for small $`N`$. Therefore, our equation provides a generalization of the standard rate equation beyond the grand canonical limit. We have shown that for rare particle production the equilibrium multiplicity is much smaller and the relaxation time is much shorter than expected from the standard rate equation. For abundant particle production, where the standard rate equation applies to first order, one can use the general rate equations to study finite-size corrections to the grand canonical results.

Our results could be of importance in the description and understanding of equilibration phenomena and equilibrium properties of partonic or hadronic matter created in heavy-ion collisions. For example, it could provide insights into the equilibration time of strange particles produced at SIS, or open charm and charmonium productions at SPS and higher energies. It may also be meaningful for transport model studies of rare particle production, especially the perturbative procedure for rare processes , where the local charge conservation requires that the production probability should be assigned to the pair instead to each particle separately.

## Acknowledgments

We acknowledge stimulating discussions with M. Gyulassy, S. Jeon and L.D. McLerran. K.R. also acknowledges comments from C. Greiner, J. Knoll, and H. Feldmeier. M.S. and Z.L. are grateful to the nuclear theory group at LBNL for supporting the summer program where this work was started. This work was supported by the Director, Office of Energy Research, Office of High Energy and Nuclear Physics, Division of Nuclear Physics, the Office of Basic Energy Science, Division of Nuclear Science, of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098, the Committee for Scientific Research (KBN-2P03B 03018), the National Science Foundation under Grant No. PHY-9870038, the Welch Foundation under Grant No. A-1358, and the Texas Advanced Research Program under Grant No. FY99-010366-0081.


text:
# 1 Contour plot of the time-like component of the net baryon four-current in planes transverse to the longitudinal axis 𝑒⃗₁^∗ of the compactness tensor (Au+Au collision at 𝐸^{𝑘𝑖𝑛}_{𝑙𝑎𝑏}=8 AGeV, 𝑏=3 fm). Left: Relativistic mean-field EoS without phase transition. Right: with phase transition to quark matter.

Observing compact quark matter droplets in relativistic nuclear collisions

K. Paecha, J. Brachmanna, M. A. Lisab, A. Dumitruc, H. Stöckera, W. Greinera

a Institut für Theoretische Physik, Universität Frankfurt a.M., Postfach 111932,
D-60054 Frankfurt am Main, Germany
b Department of Physics, The Ohio State University, 174 W. 18th Avenue,
Columbus, Ohio 43210, USA
c Department of Physics, Columbia University,
538 W. 120th Street, New York, NY 10027, USA

Oct. 2000

## Abstract

Compactness is introduced as a new method to search for the onset of the quark matter transition in relativistic heavy ion collisions. That transition supposedly leads to stronger compression and higher compactness of the source in coordinate space. That effect could be observed via pion interferometry. We propose to measure the compactness of the source in the appropriate principal axis frame of the compactness tensor in coordinate space.

A variety of signatures for the observation of the quark matter transition in relativistic heavy ion collisions have been proposed, see e.g.  and references therein. Nowadays it is widely believed that irregularities due to such a new state of matter should be seen most clearly in excitation functions, i.e. in the bombarding energy dependence, or possibly also in the impact parameter or mass-number dependence of various observables.

However, in spite of circumstantial evidence for such a novel state of matter , undisputable “irregularities” in such experimental excitation functions have not been reported to date. Most probably a convincing argument that quark matter has indeed been found requires several independent measurements of distinct observables pointing to the same conclusion - namely the onset of deconfinement at a particular bombarding energy, impact parameter and system size.

In the present paper we introduce a new observable for the onset of the phase transition. It relies on the measurement of the compactness of the source, which is related to the pressure and density of the system in the compression and expansion stage of the nucleus-nucleus collision. Observables related to the pressure do not only reflect the transient high density but rather relate to the order parameter (field), specifying in which thermodynamical state  the excited matter was created. The compactness can be identified via interferometry: The illumination of the baryon source by the pion radiation is subject to experimental scrutinity via pion interferometry measurements .

To illustrate the basic idea, let us first discuss compression shocks in heavy-ion collisions employing the Rankine Hugoniot Taub shock adiabat (RHTA) solution  of relativistic hydrodynamics,

$$W^2W_{0}^{}{}_{}{}^{2}+P\left(\frac{W}{\rho }\frac{W_0}{\rho _0}\right)=0.$$
(1)
It relates the energy per baryon in the local rest frame of the compressed matter, $`W`$, to its baryon density $`\rho `$ and pressure $`P`$. $`W_0`$ and $`\rho _0`$ are the energy (=923 MeV) and density (=0.16 fm-3) in the ground state, and the compression factor is $`\rho /\rho _0`$.

Solving for $`P`$ we obtain

$$P=W\rho _0\frac{\gamma _{\mathrm{cm}}1/\gamma _{\mathrm{cm}}}{1\gamma _{\mathrm{cm}}\rho _0/\rho }.$$
(2)
Simply speaking, the onset of the transition to quark matter at a given incident energy $`E_{\mathrm{lab}}^{\mathrm{kin}}=2(\gamma _{\mathrm{cm}}^{}{}_{}{}^{2}1)W_0`$ lets the pressure increase less rapidly with $`W`$ and/or $`\rho `$, and consequently higher compression $`\rho /\rho _0`$ can be achieved than for the case without transition. Now, as $`\rho V\pi R_\mathrm{A}^2L\rho `$ must equal $`2A`$ by virtue of baryon number conservation, the longitudinal thickness $`L`$ of the compressed matter is proportional to $`1/\rho `$. Thus, a transition to quark matter leads to a more compact system, just as compact stars with a quark matter core are more compact than pure neutron stars . Of course, in heavy-ion collisions that expectation is based on the behavior of relativistic compression shocks rather than hydrostatic and gravitational equilibrium.

To investigate quantitatively the experimental observable, we perform 3-dimensional (1-fluid) calculations of relativistic hydrodynamics for the compression as well as for the subsequent expansion stage of the reaction. That is, we solve numerically the continuity equations for the energy-momentum tensor, $`_\mu T^{\mu \nu }=0`$, and the net baryon current, $`_\mu N_B^\mu =0`$. Detailed discussions of (3+1)-d numerical solutions for hydrodynamical compression and expansion can be found e.g. in . We shall employ two different equations of state (EoS) $`P(W,\rho )`$: i) a relativistic mean field (RMF) hadron fluid  corresponding to baryons and antibaryons interacting via exchange of massive scalar and vector bosons, plus free thermal pions; the parameters of the lagrangian are fitted to the ground state of infinite nuclear matter, in particular the nuclear saturation density, the energy per nucleon, and the incompressibility. ii) the same EoS as in i) for the low-density phase, supplemented by the Bag Model equation of state with bag constant $`B^{1/4}=235`$ MeV for the quark-gluon (QG) phase. The phase coexistence region corresponding to this first-order transition is constructed employing Gibbs’ condition of phase equilibrium, $`p_{\mathrm{RMF}}(T,\mu _\mathrm{B})=p_{\mathrm{QG}}(T,\mu _\mathrm{B})`$, where $`T`$ and $`\mu _B`$ denote the temperature and the baryon-chemical potential, respectively. For example, for $`\rho =0`$ we find $`T_C170`$ MeV, while at $`T=0`$ phase coexistence sets in at $`\rho 4.6\rho _0`$. A more detailed discussion of these EoS can be found in .

As an example, we study the compactness in the reactions Au+Au at impact parameter $`b=3`$ fm and energy $`E_{\mathrm{lab}}^{\mathrm{kin}}=8A`$ GeV, where the transition to quark matter does occur (within the present model). We define the configuration space sphericity tensor, also called “compactness tensor” below, as the second moment of the net baryon current. On fixed time hypersurfaces we have

$$F_{ij}(t)=\mathrm{d}^3xx_ix_jN_B^0(t,\stackrel{}{x})\mathrm{\Theta }\left(\rho (t,\stackrel{}{x})\rho _{cut}\right).$$
(3)
We apply an additional density cut $`\rho >\rho _{cut}`$ in the integral to discard spectator matter. In the future the cuts and the hypersurface will have to be adapted to the experimental conditions. However, that is not crucial for understanding the effect.

The three eigenvalues $`f_n`$ are the solutions of the cubic equation $`det(F_{ij}f\delta _{ij})=0`$, and the eigenvectors $`\stackrel{}{e}_n`$ follow from solving the linear system of equations $`(F_{ij}f_n\delta _{ij})e_n^j=0`$. Let $`\stackrel{}{e}_1`$, $`\stackrel{}{e}_2`$ be the eigenvectors defining the reaction plane, $`\stackrel{}{e}_1`$ corresponding to the bigger of the eigenvalues $`f_1`$ and $`f_2`$. To simplify the discussion we shall assume that the matter distribution is symmetric with respect to the reaction plane.

We can now rotate the coordinate frame around $`\stackrel{}{e}_3`$ by an angle $`\mathrm{cos}\mathrm{\Theta }=\stackrel{}{e}_1\stackrel{}{e}_z`$, where $`\stackrel{}{e}_z`$ defines the longitudinal (beam-) direction in the lab frame. The rotated compactness tensor $`F^{}`$ can be written in terms of the eigenvalues $`f_n`$ and orthogonal eigenvectors $`\stackrel{}{e}_n^{}`$ as

$$F^{}=f_1\stackrel{}{e}_1^{}\stackrel{}{e}_1^{}+f_2\stackrel{}{e}_2^{}\stackrel{}{e}_2^{}+f_3\stackrel{}{e}_3^{}\stackrel{}{e}_3^{}.$$
(4)
In diagonal form, $`F^{}`$ specifies an ellipsoid in configuration space with principal axis along $`\stackrel{}{e}_n^{}`$ and radii $`\sqrt{f_n}`$. For example, cigar patterns, oriented along the z-axis, would correspond to $`f_1>f_2=f_3`$, $`\stackrel{}{e}_1^{}=\stackrel{}{e}_z`$, $`\stackrel{}{e}_2^{}=\stackrel{}{e}_x`$, $`\stackrel{}{e}_3^{}=\stackrel{}{e}_y`$. The “compactness” can now be defined as the ratio of the in-plane radii, $`\sqrt{f_2/f_1}`$.

Fig. 1 shows typical density distributions in the transverse planes of the coordinate frame where $`F^{}`$ is diagonal, eq. (4). The polar angle relative to the lab frame as well as the eccentricity

$$ϵ=\frac{f_3f_2}{f_3+f_2}$$
(5)
are shown as well. As already indicated in the introduction, we find very different eigenvalues for the two equations of state, the one without and the other with a first order phase transition included. The calculation with transition to quark matter corresponds to higher compactness of the baryon distribution. Because the incompressibility $`p/\rho `$ or $`p/(W\rho )`$ is smaller in the presence of the coexistence phase, the compactness tensor is much flatter (nearly a factor of two !) in model ii) than in model i). Moreover, our (3+1)-dimensional expansion solutions show that the ratio of the in-plane radii $`\sqrt{f_2/f_1}`$ remains much smaller in the case with phase transition, cf. Fig. 2.

This allows to measure directly the density increase in the high density stage of the reaction, if a phase transition occurs. Care must be taken that the investigated range of impact parameters constitutes a moderately small bin of centrality values. One should keep in mind that the compression factor is affected by the incompressibility evaluated on the RHTA (1), but not along a path of fixed specific entropy, because a large amount of entropy is being produced in the compression process.

The predicted change in compactness due to the reduced incompressibility associated with a first-order phase transition can be observed experimentally via the novel interferometry analysis developed recently by Lisa et al. The proposed method is quite robust and incorporates other interesting information as the configuration space tilt angles , which complement the momentum-space flow angles. It avoids cuts in tilted ellipsoids, which are not analysed in the appropriate rotated frame, and were the eccentricity and the RMS-radii are much less distinct for models i) and ii).

Acknowledgements: This work is supported by BMBF, GSI, DFG, and Graduiertenkolleg “Theoretische und Experimentelle Schwerionenphysik”. A.D. acknowledges support from the DOE Research Grant, Contract No. De-FG-02-93ER-40764. We would like to thank Stefan Scherer, Sven Soff and Steffen Bass for stimulating discussions.