text:
# N-body simulations of self-gravitating gas in stationary fragmented state

## 1 Introduction

From many studies over the last two decades, it has been established that the interstellar medium (ISM) has a clumpy hierarchical structure, approaching a fractal structure independent of scale over 4 to 6 orders of magnitude in sizes (Larson 1981, Scalo 1985, Falgarone et al 1992, Heithausen et al 1998). The structure extends up to giant molecular clouds (GMC) of 100 pc scale, and possibly down to 10 AU scale, as revealed by HI absorption VLBI (Diamond et al 1989, Faison et al 1998) or extreme scattering events in quasar monitoring (Fiedler et al. 1987, Fiedler et al. 1994). It is not yet clear which mechanism is the main responsible for this structure; it could be driven by turbulence, since the Reynolds number is very high, or self-gravity, since clouds appear to be virialized on most of the scales, with the help of magnetic fields, differential rotation, etc…

One problem of the ISM turbulence is that the relative velocities of clumps are supersonic, leading to high dissipation, and short lifetime of the structures. An energy source should then be provided to maintain the turbulence. It could be provided by star formation (stellar winds, bipolar flows, supernovae, etc… Norman & Silk 1980). However, the power-law relations observed between size and line-width, for example, are the same in regions of star-formation or quiescent regions, either in the galactic disk, or even outside of the optical disk, in the large HI extensions, where very little star formation occurs. In these last regions alternative processes of energy input must be available. A first possibilty is the injection by the galactic shear. Additionally, since there is no heating sources in the gas, the dense cold clumps should be bathing at least in the cosmological background radiation, at a temperature of 3K (Pfenniger & Combes 1994). In any case, in a nearly isothermal regime, the ISM should fragment recursively (e.g. Hoyle 1953), and be Jeans unstable at every scale, down to the smallest fragments, where the cooling time becomes of the same order as the collapse time, i.e. when a quasi-adiabatic regime is reached (Rees 1976).

In this work, we try to investigate the effect of self-gravity through N-body simulations. We are not interested in star formation, but essentially in the fractal structure formation, that could be driven essentially by gravity (e.g. de Vega et al. 1996). Our aim is to reach a quasi-stationary state, where there is statistical equilibrium between coalescence and fragmentation of the clouds. This is possible when the cooling (energy dissipated through cloud collisions, and subsequent radiation) is compensated by an energy flux due to external sources : cosmic rays, star formation, differential shear… Previous simulations of ISM fragmentation have been performed to study the formation of condensed cores, some with isolated boundary conditions (where the cloud globally collapses and forms stars, e.g. Boss 1997, Burkert et al. 1997), or in periodic boundary conditions (Klessen 1997, Klessen et al. 1998). The latter authors assume that the cloud at large scale is stable, supported by turbulence, or other processes. They follow the over-densities in a given range of scales, and schematically stop the condensed cores as sink particles when they should form structures below their resolution. We also adopt periodic boundary conditions, since numerical simulations are very limited in their scale dynamics, and we can only consider scales much smaller than the large-scale cut-off of the fractal structure. We are also limited by our spatial resolution: the smallest scale considered is far larger than the physical small-scale cut-off of the structures. Our dynamics is computed within this range of scales.

Our goal is to achieve long enough integration time for the system to reach a stationary state, stationary in a statistical sense. This state should have its energy confined in a narrow domain, thus presenting the density contrasts of a fragmented medium while avoiding gravitational collapse after a few dynamical times. Only when those conditions are fulfilled, can we try to build a meaningful model of the gas. Let us consider for example the turbulence concepts.

Theoretical attempts have been made to describe the interstellar medium as a turbulent system. While dealing with a system both self-gravitating and compressible, the standard approach has been to adapt Kolmogorov picture of incompressible turbulence. The first classical assumption is that the rate of energy transfer between scales is constant within the so-called inertial range. This inertial range is delimited by a dissipative scale range at small scales and a large scale where energy is fed into the system. In the case of the interstellar medium, the energy source can be the galactic shear, or the galactic magnetic field, or, on smaller scales, stellar winds and such. In this classical picture of turbulence, one derives the relation $`v_LL^{1/3}`$, where $`v_L`$ is the velocity of structures on scale $`L`$. If we consider now that the structures are virialized at all scales, we get the relation $`\rho _LL^{4/3}`$, where $`\rho _L`$ is the density of structures on scale $`L`$. This produces a fractal dimension $`D=1.66`$. A more consistent version is to take compressibility into account in Kolmogorov cascade; then $`\rho _Lv_L^3/Lcst`$. Adding virialization, we get the relations $`v_LL^{3/5}`$ and $`\rho _LL^{4/5}`$. This produces the fractal dimension $`D=2.2`$. It should be emphasized again that all these scenarios assume a quasi-stationary regime.

However, numerical simulations of molecular cloud fragmentation have been, so far, carried out in dissipative schemes which allow for efficient clumping (e.g Klessen et al 1998, Monaghan & Lattanzio 1991). As such, they do not reach stationary states. Our approach is to add an energy input making up for the dissipative loss. After relaxation from initial conditions we attempt to reach a fragmented but non-collapsing state of the medium that does not decay into homogeneous or centrally condensed states. It is then meaningful to compute the velocity and density fields power spectra and test the standard theoretical assumptions. It is also an opportunity to investigate the possible fractal structure of the gas. Indeed, starting from a weakly perturbed homogeneous density field, the formation time of a fractal density field independent of the initial conditions is at the very least of the order of the free-fall time, and more likely many times longer. A long integration time should permit full apparition of a fractal mass distribution, independent of the initial conditions.

This program encounters both standard and specific difficulties. Galaxies clustering as well as ISM clumping require large density contrasts, that is to say high spatial resolution in numerical simulations. This point is tackled, be it uneasily, by adaptative-mesh algorithms or tree algorithms, or multi-scale schemes ($`\text{P}^3\text{M}`$). Another problem is the need for higher time resolution in collapsed region. This is CPU-time consuming or is dealt with by multiple time steps. This point is of primary sensitivity in our simulation since we need to follow accurately the internal dynamics of the clumps and filaments to avoid total energy divergence. Finally we need a long integration time to reach the stationary state. This is directly in competition for CPU-time with spatial resolution for given computing resources.

## 2  Numerical methods

### 2.1  Self-gravity

In our simulation we use N-body dynamics with a hierarchical tree algorithm and multipole expansion to compute the forces, as designed by Barnes & Hut (1986, 1989). The number of particles is between $`\mathrm{10\hspace{0.33em}000}`$ and $`\mathrm{120\hspace{0.33em}000}`$ particles. We are limited by the fact that we need long integration time in our study.

Within the N-body framework the tree algorithm is a $`N\mathrm{ln}(N)`$ search algorithm picking the closest neighbours for exact force computation and sorting farther out particles in hierarchical boxes for multipole expansion. The size of the boxes used in the expansion is set, for the contribution of a given region to the force, by a control parameter $`\theta `$ defining the maximal angular size of the boxes seen from the point where the force is computed. Typical values used for $`\theta `$ are 0.5 to 1.0. Multipole expansion is carried out to quadrupole terms. As usual a short distance softening is used in the contribution of the closest particle to the interaction. The softening length is taken as $`1/10`$ of the inter-particle distance of the homogeneous state.

### 2.2 Boundary conditions

There are two ways to erase the finite size effect and avoid spherical collapse in a N-body simulation. We can choose either quasi-periodic boundary conditions or fully periodic boundary conditions following Ewald method (Hernquist et al 1991). In quasi-periodic conditions the interaction of two particles is computed between the first particle and the closest of all the replicas of the other particle. The replicas are generated in the periodization of the simulation box to the whole space. In the fully periodic conditions, the first particle interacts with all the replicas of the other one. The relevance of each method will be discussed in section 3.1. We will compare them for simple free-falls using the fractal dimension as diagnosis.

### 2.3  Time integration and initial conditions

The integration is carried out through a multiple time steps leap-frog scheme. Making a Keplerian assumption about the orbits of the particles, we have the following relation between time step and the local interparticle distance: $`\delta \tau (\delta l)^{\frac{3}{2}}`$. In a fractal medium, the exponent could be different. We do not take this into account since it would lead us to modify the integration scheme dynamically according to the fractal properties of the system. As the tree sorts particles in cells of size $`2^nd_0`$ according to the local density, we should use $`(2\sqrt{2})^n\delta t_0`$ time steps. For simplicity we are using $`3^n\delta t_0`$. Then, $`n`$ different time steps allow us to follow the dynamics with the same accuracy over regions with density contrasts as high as $`2^{3(n1)}`$.

We will use several types of initial conditions. The usual power law for the density power spectrum will be used for the free falls. It is usually implemented by using Zel’dovitch approximation which is valid before first shell crossing for a pressureless perfect fluid. The validity range of this approximation is difficult to assess in an N-body simulation. We use a similar implementation that does not however extend to non-linear regimes but also gives density power spectra with the desired shapes. A gaussian velocity field is chosen with a specified power spectrum $`P_v(k)k^{\alpha 2}`$, then the particles are positioned at the nodes of a grid and displaced according to the velocity field with one time step. The resulting density fluctuations spectrum is $`P_\rho (k)k^\alpha `$ according to the matter conservation.

### 2.4  Gas physics

There can be several modelizations for the ISM: either it is considered as a continuous and fluid medium and simulated through gas hydrodynamics (with pressure, shocks, etc..) or, given the highly clumped nature of the interstellar medium, and its highly inhomogeneous structure (density variations over 6 orders of magnitude), it can be considered as a collection of dense fragments, that dissipate their relative energy through collisions. Our choice of particle dynamics assumes that each particle stands for a clump and neglects the interclump medium. Another reason is that we expect a fractal distribution of the masses. This means a non-analytical density field which is not to be easily handled by hydrodynamical codes.

The dissipation enters the dynamics through a gridless sticky particles collision scheme. The collision round is periodically computed every 1 to 10 time steps. The frequency of this collision round is one way to adjust the strength of the dissipation. Since, as the precise collision schemes will reveal , we adopt a statistical treatment of the dissipation, it is not necessary to compute the collision round at each time step. Two schemes are investigated. The first uses the tree search to find a candidate to collision with a given particle, in a sphere whose radius $`r_c`$ is a fixed parameter. Inelastic collisions are then computed, dissipating energy but conserving linear momentum. In the second scheme, each particle has a probability to collide proportional to the inverse of the local mean collision time. If it passes the probability test, the collision is computed with the most suitable neighbour. In this second scheme, dissipation may occur even in low density regions; in the first it occurs only above a given density threshold. In practice this implies that, in the first scheme, dissipation happens only at small scales, while in the second, it happens on all scales but is still stronger at small scales.

Special attention must be paid to the existence of a small scale cutoff in the physics of the system. Indeed, at very small scale ( $`20`$ AU ), the gas is quasi-adiabatic. Its cooling time is then much longer than the isothermal free fall time at the same scale. Moreover, the mean collision time between clumps at this scale, in a fractal medium, is much shorter than the cooling time. As a result, collisions completely prevent the already slowed down processes of collapse and fragementation in the quasi-adiabatic gas: the fractal structure is broken at the corresponding scale. To take these phenomena into account in our simulation, we need to introduce a large-density cutoff. This is achieved by computing elastic or super-elastic collisions at scale $`<0.5r_c`$, instead of inelastic collisions. The practical implementation follows the same procedure as for the first inelastic collision scheme. In the case of super-elastic collisions, the overall energy balance is still negative at small scales (due to a volume effect), however we do introduce a mechanism of energy injection at small scale. This choice can be furthermore justified by physical considerations.

When there is no energy provided by star formation (for example in the outer parts of galaxies, where gas extends radially much beyond the stellar disk), the gas is only interacting with the intergalactic radiation field, and the cosmic background radiation. The latter provides a minimum temperature for the gas, and plays the role of a thermostat (at the temperature of 2.76K at zero redshift). The interaction between the background and the gas can only occur at the smallest scale of the fragmentation of the medium, corresponding to our small scale cutoff; the radiative processes involve hydrogen and other more heavy elements (Combes & Pfenniger 1997). To maintain the gas isothermal therefore requires some energy input at the smallest scale.

Several schemes for the energy input have been tried. Reinjection at small scale through added thermal motion has proved unable to sustain the system in any other state than an homogeneous one. Then, turning to the turbulent point of view, we have tried reinjection at large scale. Two main methods have been tested; reinjection through a large scale random force field, and reinjection through the action of the galactic shear. Their effect will be discussed in section 5.

Finally, the statistical properties of the system can be described in many different way. We will use the correlation fractal dimension as diagnosis. Definition and example of application of this tool are given in the appendix.

## 3  Simple 3-D free falls

### 3.1  The choice of the boundary conditions

A fractal structure such as observed in the ISM or in the galaxies distribution obeys a statistical translational invariance. If we run simulations with vacuum boundary condition, this invariance is grossly broken. On the other hand the quasi-periodic or fully-periodic conditions restore this invariance to some extent. They have been used in cosmological simulations. Hernquist et al (1991) have compared the two models to analytical results and found that the fully-periodic model simulates more accurately the self-gravitating gas in an expanding universe. This is in agreement with expectations, since uniform expansion is automatically included in a fully-periodic model, while it is not in a quasi-periodic one. Klessen (1996,1998) applied Ewald method to simulations of the interstellar medium. This is not natural since expansion is necessary to validate Ewald method. Nevertheless one can argue that the dynamics is not altered by the use of Ewald method in regions with high density contrasts.

Hernquist’s study was on the dynamics of the modes. We are more interested in the fractal properties. As far as we are concerned the question is wether the same initial conditions produce the same fractal properties in a cosmological framework (fully-periodic conditions) and in a ISM framework (quasi-periodic conditions). We will try to answer this question in section 3.3.

### 3.2  The choice of the initial conditions

The power spectrum of the density fluctuations is most commonly described as a power law in cosmological simulations. Then, different exponents in the power law produce different fractal dimensions. The comoslogist can hope to find the fractal dimension from observations and go back to the initial power spectrum. In the case of the interstellar medium however, it is unlikely that the fractal dimension is the result of initial conditions since the life-time of the medium is much longer than the dynamical time of the structures. In this regard, all initial conditions should be erased and produce the same fractal dimension. We will check whether this is the case for power law density spectra with different exponents. Simple free falls with such initial conditions will allow the comparison between the two types of periodic boundary conditions and will allow us to compare between 3D and 2D models and different dissipation schemes.

### 3.3  Simulations and results

The four first simulations are free falls from power law initial conditions in either fully or quasi-periodic boundary conditions for two different power law exponents: $`\alpha =1`$ and $`\alpha =2`$. About $`\mathrm{120\hspace{0.17em}000}`$ particles were used in the simulations. The time evolution of the particle distribution for $`\alpha =2`$ and fully periodic condition is plotted in Fig. 1. If we compare the qualitative aspect of the matter distribution with those found in SPH simulations (Klessen 1998), a difference appears. Filaments are less present in this N-body simulation than in SPH simulations. We believe that this is due to the strongly dissipative nature of SPH simulations. Indeed filaments form automaticaly, even without gravity, with $`P(k)k^2`$ initial conditions for the density field. Then, the strong dissipation of SPH codes is necessary to retain them; $`N`$-body codes seem to show that gravity by itself is not enough.

#### 3.3.1 Fractal dimension

For each of the four simulations, the fractal dimension as a function of scale has been computed at different times of the free fall. Results are summarized in Fig. 2.

The fractal dimension of a mathematical fractal would appear either as an horizontal straight line or as a curve oscillating around the fractal dimension if the fractal has no randomization (like a Cantor set). For our system, it appears that at sufficiently early stages of the free fall, and for scales above the dissipative cutoff, the matter distribution is indeed fractal. However the dimension goes down with time to reach values of 2 for $`\alpha =2`$, and between $`2`$ and $`2.5`$ for $`\alpha =1`$. This is the fractal dimension for the total density, not for the fluctuations only. Stabilization of the fractal dimension marginally happens before the dissipation breaks the fractal.

One conclusion is however reasonable: as far as the fractal dimension is concerned, the difference between fully-periodic boundary conditions and quasi-periodic boundary conditions is not important. The deviation is a bit stronger in the $`\alpha =1`$ case, which is understandable since $`\alpha =1`$ produces weaker density contrasts in the initial conditions than $`\alpha =2`$ and so the dynamics is less decoupled from the expansion. The effect of this conclusion is that we will not use the fully-periodic boundary conditions in further simulations since they require more CPU time and have no physical ground in the ISM framework. Clump mass spectra are given in appendix B.

## 4 Simple 2-D free falls

We consider a 2-D system where particles have only two cartesian space coordinates but still obey a $`\frac{1}{r^2}`$ interaction force.

### 4.1 Motivations for the two dimensional study

Turning to 2-D simulations allows access, for given computing resources, to a wider scale-range of the dynamics. Or, for given computing resources and a given scale-range, to larger integration times. The later is what we will need in the simulation with energy injection. In practice we will choose to use $`10^4`$ particles in 2-D simulations to access both a larger integration time and a wider scale range. Indeed, between a 3-D simulation with $`N=10^5`$ particles and a 2-D simulation with $`10^4`$ particles one improves the scale range from $`[1,N^{\frac{1}{3}}=46]`$ to $`[1,N^{\frac{1}{2}}=100]`$. At the same time one can increase the integration-time by a factor of $`10`$ at equal CPU cost.

We can also add a physical reason to these computational justifications. The simulations with an energy injection aim to describe the behaviour of molecular clouds in the thin galatic disk. This medium has a strong anisotropy between the two dimensions in the galactic plane and the third one. As such it is a suitable candidate for a bidimensional modeling.

Moreover, studying a bidimensional system will allow us to check how the fractal dimension, dependent or not of the initial conditions, is affected by dimension of the space.

### 4.2  Simulations with 2 different dissipative schemes

Here again initial conditions with a power law density power spectrum are used. Exponent $`\alpha =2`$ is chosen to allow a comparison with the 3-D simulations. Boundary conditions are quasi-periodic. We have carried out two simulations with the two different dissipative schemes described in section 2.4 . The fractal dimensions at different stages of each of these two simulations are shown in Fig. 3.

In both simulations the fractal dimension of the homogeneous initial density field is 2 at scales above mean inter-particle distance and 0 below. Then the density fluctuations develop and the fractal dimension goes down to about $`1.51.6`$ in both simulations. This value was not reached in 3-D simulations. This shows that the dimension of the space has an influence on the fractal dimension of a self-gravitating system evolved from initial density fields with a power law power spectrum. Indeed spherical configurations are favoured in 3 dimensions and may produce a different type of fractal than cylindrical configurations which are favored in two dimensions.

Another point is that the fractal dimension of the system at small scales is indeed sensitive to the dissipative scheme. The second scheme even alters the fractality of large scales. So even if we believe that the fractality of a self-gravitating system is controlled by gravity, we must keep in mind the fact that the dissipation has an influence on the result. In this regard, discrepancies should appear between $`N`$-body simulations including SPH or other codes.

The clump mass spectra are not conclusive due to insufficient number of clumps formed with 10609 particles.

## 5  2-D simulations with an energy source

As already stated, the ISM is a dissipative medium. Structures radiate their energy in a time-scale much shorter than the GMCs typical life time inferred from observations. Therefore we are led to provide an energy source in the system to obtain a life time longer than the free fall time.

### 5.1  The choice of the energy source

Different types of sources are invoked to provide the necessary energy input: stellar winds, shock waves, heating from a thermostat, galatic shear… We have tested several possibilities.

Our first idea was to model an interaction with a thermostat by adding periodically a small random component to the velocity field. We have tried to add it either to all the particles or only to the particle in the cool regions. If enough energy is put into this random injection, the mono-clump collapse is avoided. However the resulting final state is not fractal, nor does it show density structures on several different scales. It consists in a few condensed points, like ”blackholes”, in a hot homogeneous phase.

The second idea is to model a generic large scale injection by introducing a large scale random force field. We have tried stationary and fluctuating fields. It is possible in this scheme, by tuning the input and the dissipation, to lengthen the medium life time in an ”interesting” inhomogenous state by a factor 2 or 3. However the final state, after a few dynamical times, is the same as for thermal input.

We would like to emphasize that, in all the cases just mentioned, the existence of super-elastic collisions at very short distances, are not at all able to avoid the collapse of ”blackholes”. They are indeed beaten by inelastic collisions which happen more frequently. Thus the density cutoff introduced by superelastic collision, is not able all by itself to sustain or destroy clumps.

Only one of the schemes we have tested avoids the blackhole/homogeneity duality. In this scheme, energy is provided by the galactic shear. It produces an inhomogeneous state whose life time does not appear to be limited by a final collapse in the simulations. We now describe this scheme in detail.

### 5.2  Modelisation of the galactic shear

The idea is to consider the simulation box as a small part of a bigger system in differential rotation. The dynamics of such a sub-system has been studied in the field of planetary systems formation (Wisdom and Tremaine 1988). Toomre (1981) and recently Hubert and Pfenniger (1999) have applied it to galactic dynamics.

The simulation box is a rotating frame with angular speed $`\mathrm{\Omega }_0`$ which is the galactic angular speed at the center of the box. Coriolis force acts on the particles moving in the box. An additional external force arises from the discrepancy between the inertial force and the local mean gravitational attraction of the galaxy. The balance between these two forces is achieved only at points that are equally distant from the galactic center as the center of the box. Let us consider a 2-dimensional case. If $`y`$ is the orthoradial direction, and $`x`$ the radial direction the equations of motion are written:

$`\ddot{y}`$ $`=`$ $`2\mathrm{\Omega }_0\dot{x}+F_y`$
$`\ddot{x}`$ $`=`$ $`2\mathrm{\Omega }_0\dot{y}2\mathrm{\Omega }_0r_0{\displaystyle \frac{d\mathrm{\Omega }}{dr}}|_{r_0}x+F_x`$
The $`F_{x,y}`$ are the projections of the internal gravitational forces. The shear force acting in the radial direction is able to inject into the system energy taken from the rotation of the galaxy as a whole.

Some modifications must also be brought to the boundary condition. It is not consistent with the differential rotation to keep a strict spatial periodicity. To take differential rotation into account, layers of cells at different radii must slide according to the variation of the galactic angular speed between them. One particle still interacts with the closest of the replicas of another particle (including the particle itself). This is sketched on Fig. 4.

### 5.3  Results

We have computed simulations in two dimensions with 10609 particles. The direction orthogonal to the galactic disk is not taken into account. A tuning between dissipation rate and shear strength is necessary to obtain a structured but non collapsing configuration. However it is not a fine tuning at all. We found out that, in a range of values of the shear, the dissipation adapts itself to compensate the energy input so that we didn’t encounter any limitation in the life time of the medium. The simulations have been carried out over more than 20 free-fall times. Snapshots from a simulation are plotted in Fig. 5.

#### 5.3.1 Formation of persistent structures

We can see that tilted stripes appear in the density field. This structure appears also in simulations by Toomre (1981) and Hubert and Pfenniger (1999). The new feature is the intermittent appearance of dense clumps within the stripes. Theses clumps remain for a few free fall times of the total system, then disappear, torn apart by the shear. This phenomenon is unknown in the simple free fall simulations where a clump, once formed, gets denser and more massive with time. It also appeared that the tearing of the clumps by the shear happens only if super-elastic collisions enact an efficient density cut-off. If we use elastic collisions only, below the dissipation scale, thus enacting a less efficient density cut-off, the clumps persist. However they do not collapse to the dense blackhole states encountered with different large scale energy input schemes. So the combination of the galactic shear and super-elastic collisions below the dissipation scale is necessary to obtain destroyable clumps, while the shear alone produces rather stable clumps, and the super-elastic collisions alone cannot avoid complete collapse.

While we cannot consider the density field as a fractal struture, it is indeed the beginning of a hierarchical fragmentation since we have two levels of structures. However the stripes are mainly the result of the shear. On the other hand the clumps are the simple gravitational fragmentation of the stripes. These two structures are not created by the exact same dynamical processes.

It must also be noted that the shear plays a direct role at large scale (of the order of 10 pc), where the galactic tidal forces are comparable to the self-gravity forces of the corresponding structures. But since the structures are fragmenting hierarchically, it is expected that its effect propagates and cascades down to the smallest scales.

#### 5.3.2 Fractal analysis

The fractal dimension of the medium is plotted at different times in Fig. 6. Two characteristic scales are clearly visible. At large scale, the scale of the stripes, the dimension is about 1.9 with some fluctuations in time. At small scale, the dimension fluctuates strongly due to the appearance and disappearance of clumps. This shows through the appearance of a bump at the clump size. This size is about the dissipation scale. It is important to mention that the small scale range, while under the cutoff of the initial homogeneous condition, is not under our dynamical resolution. The mean inter-particle distance decreases from its initial value and is consistently followed by the $`N`$-body dynamics.

#### 5.3.3 Effect of the initial conditions

We have already emphasized that in GMC simulations, the resulting structure should be independent of the initial conditions. We have performed the simulation for various values of the exponant of the power spectrum. In the range $`[0.5,0.5]`$ the long time behaviour is the one described above, thus satisfying the required independence. However, for more steep spectrum value like -2, the medium tend to collapse early into large clumps, thereby preventing the shear to act efficiently.

We see that the behaviour is independent of the initial conditions as long as the continuity of the medium at scale larger than the GMC size (the clump size in the simulation) is not broken so that the shear can act to sustain the GMC against collapse.

## 6 Discussion and Conclusion

We have shown that, in simulations without energy input, the evolution of the fractal dimension depends on the initial conditions and also, to some extent, on the dissipation schemes. Then we have investigated the effects of an energy input on the dynamics of the system. This input is necessary to prevent the total collapse of the system. Moreover, we do not want the input to destroy the inhomogeneous fractal structure of the medium, as done by random reinjections (thermal bath, random force field). Among the solutions we have tested, only the highly regular force field provided by the galactic shear preserves an inhomogeneous state. However this state is quite strongly constrained by the geometry of the shear. Stripes appear as a result. Interestingly, clumps are formed in these stripes and can be destroyed if we have an efficient enough density cutoff .

This state does not have a well-defined fractal dimension. Stripes are the effect of a large scale action on the system and within our dynamic range in scale, we cannot get rid of this boundary effect. If we had more resolution, the clumps could fragment, and we might reach a scale domain free of boundary effects. One of our prospects is to follow the inner dynamic of a clump with simulations in recursive sub-boxes going down in scale. This can produce a great dynamic scale range but only in a small region of space.

## Appendix A  Computing the fractal dimension

By fractal dimension we mean more accurately the correlation fractal dimension. We will use the notation $`D_f`$, defined as follows : if we consider a fractal set of points of equal weight located at $`\text{x}_j`$, the fractal dimension of the ensemble obeys the relation :

$$\underset{r0}{lim}r^{D_f}=\underset{r0}{lim}\underset{j}{}_{|\text{x}\text{x}_i|<r}\delta (\text{x}\text{x}_j)𝑑\text{x}_i.$$

The brackets stand for an ensemble average. In computational applications and in physical systems, the $`r0`$ limit is not reached. However the relation should hold at scales where the boundary effects created by the finite size of the system are negligible.

In practice, for a non-fractal ensemble of points, $`D_f`$ usually depends on $`r`$ (but not always). For a fractal it is independent of $`r`$, or at least it oscillates around a mean value. The latter case can happen for a set with a strong scale periodicity like the Cantor set. Examples are given in Fig. 7. As the above remarks show, $`D_f`$ is not a definitive criterion of fractality.

According to the definition, the computation of the fractal dimension is based on the following method. $`r^{D_f}`$ is the mean mass in a sphere of radius $`r`$ centered on a particle. This mass is actually computed for each particle using the tree search, then it is averaged and the value $`D_f(r)`$ is derived. The average can be made on a subset of points to lower the computation cost, the result is then noisier.

## Appendix B Clump mass spectrum

Although the results are not very conclusive, we give the mass distribution of the clumps in the 3-D simple freefall, to allow comparison with observational data and other numerical simulations. Plotting $`\mathrm{log}_{10}(N)`$ against $`\mathrm{log}_{10}(m)`$, the slope infered from observation data is $`0.5`$. Thus the clump mass spectrum behaves as a power law : $`\frac{dN}{dm}m^{1.5}`$.

To produce the mass spectrum, the first step is to define individual clumps. We have used an algorithm very similar to the one proposed by Williams et al (1994). We compute the density field from an interpolation of the distribution of the particles with a cloud-in-cell scheme. It should be mentione d that this procedure tends to produce an artificially high number of very small clumps associated to isolated particles or pair of particules which should not be considered as clumps at all.

On Fig. 8 the clump mass spectrum is given at different times of the evolution and for the two types of initial conditions ($`\alpha =1,2`$). Quasi-periodic boundary conditions are used. As time increases, as we can see, the clump mass spectrum slope decreases to reach $`0.38\pm 0.03`$ for $`\alpha =1`$ and $`0.18\pm 0.03`$ for $`\alpha =2`$. At later times, the power law is broken by the final dissipative collapse. According to this diagnosis, the case $`\alpha =1`$ is closer to the observed values. However, once again, this simulation provides only a transient state which is unlikely to be a good model of a GMC.


format: text

text:
# References

Percolation and Magnetization for

Generalized Continuous Spin Models

Santo Fortunato and Helmut Satz

Fakultät für Physik, Universität Bielefeld

D-33501 Bielefeld, Germany

Abstract:

For the Ising model, the spin magnetization transition is equivalent to the percolation transition of Fortuin-Kasteleyn clusters; this result remains valid also for the conventional continuous spin Ising model. The investigation of more general continuous spin models may help to obtain a percolation formulation for the critical behaviour in $`SU(2)`$ gauge theory. We therefore study a broad class of theories, introducing spin distribution functions, longer range interactions and self-interaction terms. The thermal behaviour of each model turns out to be in the Ising universality class. The corresponding percolation formulations are then obtained by extending the Fortuin-Kasteleyn cluster definition; in several cases they illustrate recent rigorous results.

Introduction

Percolation theory  has been successfully applied to describe continuous thermal phase transitions as geometrical transitions: at the critical threshold, suitably defined clusters become infinite (see  \- ). It is moreover possible to establish a correspondence between thermal and cluster variables, such that the critical exponents of corresponding variables coincide . This provides an intuitive geometric description of the mechanism of the thermal transition; e.g., in the case of the Ising model, we can identify clusters as magnetic domains, so that magnetization sets in when they span the lattice. In a recent study , we have tried to extend this picture to the confinement-deconfinement transition of $`SU(2)`$ gauge theory. This theory leads to a continuous phase transition which belongs to the universality class of the Ising model, due to the common $`Z(2)`$ global symmetry of the corresponding Hamiltonians . However, our result is valid only in the strong coupling limit of $`SU(2)`$, where the partition function can be well approximated  by that of a continuous spin Ising model. In this model, the correspondence between percolation and thermal variables can be rigorously established . For a general study of $`SU(2)`$ gauge theory, the weak coupling limit must be taken into account, and there the approximations used in  are no longer valid.

To address the deconfinement problem more generally, it may be helpful to look for an effective theory for $`SU(2)`$. An interesting attempt in this direction was proposed some time ago in the search for the fixed point of the theory by means of block-spin transformations . There the original Polyakov loop configurations were projected onto configurations of Ising-like spins according to the sign of the Polyakov loop at each lattice site. This assumes that for the criticality only the global $`Z(2)`$ symmetry of the theory is important. As a result of this projection one obtains an effective Hamiltonian with only short-range couplings, although some of the resulting couplings correspond to interactions beyond the fundamental nearest-neighbour form. Moreover, the assumption that only the $`Z(2)`$ symmetry of the theory is relevant must be tested.

In this work we investigate several spin models in order to get a better understanding of the behaviour of the $`SU(2)`$ effective theory. Our investigations are based on lattice Monte Carlo simulations of the various models.

We focus on continuous spin Ising models, in which the individual spins $`s_i`$ at each lattice site can take on all values in the finite range $`[1,1]`$, with the distribution of the spins governed by a spin distribution function $`f(s_i)`$. For $`f(s_i)=\delta (|s_i|1)i`$ and only nearest-neighbour (NN) interactions, we recover the conventional Ising model. For $`f(s_i)=1i`$ and only NN interactions, we obtain the continuous spin model studied in . Here we will consider four more general models of this type ($`d`$ denotes the space dimension):

A) $`d=2`$, only NN interactions, but a non-uniform spin distribution $`f(s_i)`$;

B) $`d=2`$, $`f(s_i)=1i`$, NN and diagonal next-nearest-neighbour (NTN) interaction (Fig. 1);

C) $`d=3`$, $`f(s_i)=1i`$, NN and two types of NTN interactions (see Fig. 2);

D) as case C), but including an additional self-interaction term proportional to $`s_i^2i`$.

We will first present a detailed numerical investigation of the thermal critical behaviour and show that each model belongs to the same universality class as the Ising model. Subsequently, we provide a suitable cluster definition in order to interpret the thermal phase transition of the models as a geometrical percolation transition. Concerning the second issue, it was recently proven rigorously that the result of  can be extended to the models A, B and C . We believe, however, that it is useful to illustrate it by exploiting the power of the finite-size scaling analysis. On the other hand, model D has not been investigated so far; as self-interactions play an important role in gauge theories, it is interesting to examine what happens in this case.

Model A: Spin Distribution Function

The spin variable $`s`$ can take all values in $`[1,1]`$. It is convenient to rewrite it in the form

$$s=sign(s)\sigma ,$$
(1)
separating the sign from the amplitude of the spin. Since we want to retain the $`Z(2)`$ symmetry of the Hamiltonian, the signs of the spins are equidistributed. In contrast, the amplitudes can in general be weighted in different ways by choosing a distribution function $`f(\sigma )`$. The partition function of such a model then has the form

$$Z(T)=\underset{i}{}_0^1𝑑\sigma _if(\sigma _i)\mathrm{exp}\{\kappa \underset{i,j}{}s_is_j\},$$
(2)
where $`\kappa J/T`$, with $`J`$ denoting the spin coupling energy and $`T`$ the temperature. Since we want to study models with continuous transitions, the choice of the function $`f(\sigma )`$ is not arbitrary. It can be proved that it must obey certain regularity conditions, which are not very restrictive, however . In , the amplitudes were equidistributed, with $`f(\sigma )=1`$. Here we consider the following form for $`f(\sigma )`$:

$$f(\sigma )=\sqrt{1\sigma ^2}$$
(3)
which is the Haar measure of the SU(2) group; it appears in the strong coupling approximation of SU(2) . To define a cluster, we use exactly the same procedure as in , i.e., we join nearest-neighbouring spins of the same sign with the bond probability $`p=1exp(2\kappa \sigma _i\sigma _j)`$, where $`\sigma _i`$ and $`\sigma _j`$ are the amplitudes of the two spins. We use free boundary conditions and define percolation when a cluster spans the lattice in each of the two directions. This is done to exclude the possibility that, due to the finite size of our lattices, one could find more than one percolation cluster, which would make the definition of our variables ambiguous. We shall use this definition for all the models to be considered here111In three dimensions even this definition of a spanning cluster does not exclude the possibility of having more than one such cluster in the same configuration. Nevertheless the occurrence of these cases is so small that we can safely ignore them..

The percolation variables to be studied are: the percolation strength  $`P`$,

$$P=\frac{\mathrm{𝑠𝑖𝑧𝑒}\mathrm{𝑜𝑓}\mathrm{𝑡ℎ𝑒}\mathrm{𝑝𝑒𝑟𝑐𝑜𝑙𝑎𝑡𝑖𝑛𝑔}\mathrm{𝑐𝑙𝑢𝑠𝑡𝑒𝑟}}{\mathrm{𝑣𝑜𝑙𝑢𝑚𝑒}\mathrm{𝑜𝑓}\mathrm{𝑡ℎ𝑒}\mathrm{𝑙𝑎𝑡𝑡𝑖𝑐𝑒}}.$$
(4)
and the average cluster size  $`S`$,

$$S=\frac{\underset{s}{}n_ss^2}{_sn_ss}.$$
(5)
Here $`n_s`$ is the number of clusters of size $`s`$ and the sums exclude the percolating cluster. Besides these two, we evaluate another very useful variable. For a given lattice size and a value of $`\kappa `$, we determine the relative fraction of configurations containing a percolating cluster. This variable, that we shall call percolation cumulant and indicate as $`\gamma _r(\kappa )`$, is a scaling function, analogous to the Binder cumulant $`g_r(\kappa )`$ .

In addition to these percolation variables, we also measure the lattice averages of energy and magnetization, in order to study the thermal transition as well. From the magnetization $`M`$ we calculate in particular the physical susceptibility $`\chi `$,

$$\chi (\kappa )=V(M^2M^2)$$
(6)
and the Binder cumulant $`g_r`$,

$$g_r(\kappa )=\frac{M^4}{M^2^2}.$$
(7)
We have performed simulations on four lattice sizes, from $`64^2`$ to $`240^2`$. Our algorithm consists in heatbath steps for the update of the spin amplitudes followed by Wolff-like cluster updates for the flipping of the signs. That is basically the same method as used in , although the heatbath procedure is slightly modified to take into account the presence of the distribution function $`f(\sigma )`$. The proof that the detailed balance condition is satisfied for such a cluster update is obtained following the derivations in  \- . Flipping the spins reduces drastically the autocorrelation time, which allowed us to collect extensive data and thus reduce the errors.

Fig. 3 shows a comparison between the thermal variable $`g_r(\kappa )`$ and the percolation cumulant $`\gamma _r(\kappa )`$, both as functions of the temperature variable $`\kappa `$, for different lattice sizes. The points where the lines cross are the infinite volume thresholds for the thermal and the geometrical transition, respectively. The agreement between the two points is remarkable.

Since the percolation cumulant is a scaling function, we can already get an indication about the critical exponents of the percolation transition. Given the critical point and the exponent $`\nu `$, a rescaling of the percolation cumulant as a function of $`(\kappa \kappa _{cr})L^{1/\nu }`$ should give us the same function for each lattice size ($`L`$ is the linear dimension of the lattice). Figs. 4 and 5 show the rescaled percolation cumulant, using $`\kappa _{cr}=1.3888`$ with the random percolation exponent $`\nu =4/3`$ and Ising exponent $`\nu =1`$, respectively. The figures clearly show scaling for the Ising exponent and no scaling for the random percolation exponent.

The exponent $`\beta `$ governs percolation strength and magnetization, while $`\gamma `$ governs cluster size and susceptibility. To determine them, we have performed high-statistics simulations around the critical point, with the number of measurements for each value of the coupling varying from 50000 to 100000. We have then used the $`\chi ^2`$ method  to determine the values of the exponents. The results are shown in Table I, where the first row shows the calculated percolation results, the second the corresponding magnetization calculations, and the third the analytically known Ising model values. It is evident that the percolation behaviour coincides fully with the thermal critical behaviour. This conclusion holds in general for all admissable distribution functions, as it was shown in .

Model B: Next-to-Nearest Neighbour Interactions

We now want to study how to correctly define clusters in the presence of more than the standard nearest-neighbour interactions. We have now two terms, with a Hamiltonian of the form

$$=J_{NN}\underset{i,j}{\overset{NN}{}}s_is_jJ_{NTN}\underset{i,j}{\overset{NTN}{}}s_is_j$$
(8)
where the first sum describes nearest-neighbour and the second diagonal next-to-nearest neighbour interactions (Fig. 1). Since longer range interactions are generally weaker, we fix the ratio between the two couplings at $`J_{NN}/J_{NTN}=10`$; however, we do not believe that our results depend on the choice of the couplings, as long as both are ferromagnetic. To define clusters, we now extend the Coniglio-Klein method and define for each two spins $`i,j`$ of the same sign, for NN as well as NTN, a bond probability

$$p_B^x=1\mathrm{exp}\{2\kappa _x\sigma _i\sigma _j\},$$
(9)
where $`x`$ specifies $`\kappa _{NN}=J_{NN}/T`$ and $`\kappa _{NTN}=J_{NTN}/T`$, respectively. We recall that in model B, as well as in C and D, we assume a spin equidistribution, i.e., $`f(\sigma )=1`$.

We have studied model B using two different Monte Carlo algorithms, in order to test if a Wolff-type algorithm can also be applied in the presence of NTN interactions. The first is the standard Metropolis update, while the second alternates heat bath steps and a generalized Wolff flipping, for which the clusters are formed taking into account both interactions. The generalization of the cluster update is trivial. After several runs, some with high statistics, we found excellent agreement with the Metropolis results in all cases. So, the mixed algorithm with heat bath and Wolff flippings appears to remain viable also in the presence of more than the standard NN interaction. Subsequently we have therefore used this mixed algorithm. The lattice sizes for model B ranged from $`100^2`$ to $`400^2`$.

As for model A, we now show for model B the comparison between percolation cumulant and Binder cumulant, in order to test that the critical points coincide (Fig. 6). We then again rescale the percolation cumulant, using the critical $`\kappa `$ determined in Fig. 6 together with the exponent $`\nu `$ from the $`d=2`$ random percolation transition ($`\nu =4/3`$) and the $`d=2`$ Ising model ($`\nu =1`$). Figs. 7 and 8 indicate that again the correct exponent is that of the Ising model.

The determination of the other two exponents $`\beta `$ and $`\gamma `$ by means of the $`\chi ^2`$ method confirms that indeed both the thermal and the geometrical exponents belong to the Ising universality class (Table II).

Model C: Extension to Three Dimensions

We have here repeated the study for a $`d=3`$ model with three different interactions (Fig. 2). To fix the model, we have to specify the ratios of the nearest-neighbour coupling $`\kappa _{NN}`$ to $`\kappa _{NTN}`$ and $`\kappa _{diag}`$. We chose them to be 10:2 and 10:1, respectively. Our $`d=3`$ calculations are performed on lattices ranging from $`12^3`$ to $`40^3`$.

Also here we have first compared the results from a mixed algorithm of the same kind as for the previous case to those from a standard Metropolis algorithm; again, the agreement is very good. Figs. 9, 10 and 11 then show the comparison of the thresholds and the scaling of the percolation cumulant. As before, the correspondence between percolation and thermal variables is evident (Table III).

Model D: Adding Self-Interaction

From what we have seen up to now, it seems to be clear that the correct cluster definition can readily be extended to models including several (ferromagnetic) spin-spin interactions. However, such terms are not the only possible interactions in a model with Z(2) symmetry and a continuous transition. There could be anti-ferromagnetic spin-spin couplings as well as multispin terms, coupling an even number of spins greater than two (four, six, etc.). Moreover, since the spins are continuous, the presence of self-interaction terms is possible, determined by $`s^2`$, $`s^4`$, etc. The treatment for antiferromagnetic and multispin couplings so far remains an open question. In contrast, self-interactions are not expected to play a role in the cluster building, since such terms do not relate different spins. Therefore, we test a cluster definition ignoring any self-interaction term.

We thus consider in Model D the same interactions as in Model C, plus a term proportional to $`J_0_is_i^2`$. We chose a negative value for the self-interaction coupling $`J_0`$; this is the more interesting case since the corresponding interaction tries to resist the approach of the system to the ground state at low temperatures ($`\sigma =1`$ everywhere). The ratios of the NN coupling to the others were chosen as $`J_{NN}:J_{NTN}:J_{diag}:|J_0|=6:2:1:2`$.

Again we first verify the viability of the mixed algorithm and then determine the rescaled percolation cumulant (Fig. 12) and the critical exponents (Table IV). It is evident that the percolation and the thermal transition again fall into the Ising universality class.

Conclusions

We have performed a complete investigation of several continuous spin models, obtained from the simple model studied in  by introducing a distribution function for the spin amplitudes, longer-ranged interactions and self-interactions. The thermal behaviour leads in all cases to the Ising universality class. All systems admit equivalent percolation formulations; for models A, B and C, we illustrate numerically the analytical results found in . All the features of the models we studied are expected to enter in the further extension to SU(2) gauge theories, and the fact that they do not interfere with a percolation formulation of critical behaviour provides support for such efforts.

Acknowledgements

It is a pleasure to thank D. Gandolfo for helpful discussions. We would also like to thank the TMR network ERBFMRX-CT-970122 and the DFG Forschergruppe Ka 1198/4-1 for financial support.

|  | $`\kappa _{cr}`$ | $`\beta /\nu `$ | $`\gamma /\nu `$ | $`\nu `$ |
| --- | --- | --- | --- | --- |
| Percolation | 1.3887$`{}_{0.0001}{}^{}{}_{}{}^{+0.0002}`$ | 0.128$`{}_{0.010}{}^{}{}_{}{}^{+0.007}`$ | 1.754$`{}_{0.008}{}^{}{}_{}{}^{+0.007}`$ | 0.99$`{}_{0.02}{}^{}{}_{}{}^{+0.03}`$ |
| Magnetization | 1.3888$`{}_{0.0003}{}^{}{}_{}{}^{+0.0002}`$ | 0.121$`{}_{0.006}{}^{}{}_{}{}^{+0.008}`$ | 1.745$`{}_{0.007}{}^{}{}_{}{}^{+0.011}`$ | 1.01$`{}_{0.03}{}^{}{}_{}{}^{+0.02}`$ |
| Ising Model |  | 1/8 | 7/4 | 1 |

Table I. Thermal and percolation thresholds and exponents of Model A.

|  | $`\kappa _{cr}`$ | $`\beta /\nu `$ | $`\gamma /\nu `$ | $`\nu `$ |
| --- | --- | --- | --- | --- |
| Percolation | 0.9708$`{}_{0.0002}{}^{}{}_{}{}^{+0.0002}`$ | 0.129$`{}_{0.009}{}^{}{}_{}{}^{+0.008}`$ | 1.752$`{}_{0.011}{}^{}{}_{}{}^{+0.009}`$ | 1.005$`{}_{0.020}{}^{}{}_{}{}^{+0.012}`$ |
| Magnetization | 0.9707$`{}_{0.0002}{}^{}{}_{}{}^{+0.0003}`$ | 0.124$`{}_{0.005}{}^{}{}_{}{}^{+0.007}`$ | 1.747$`{}_{0.007}{}^{}{}_{}{}^{+0.009}`$ | 0.993$`{}_{0.010}{}^{}{}_{}{}^{+0.014}`$ |
| Ising Model |  | 1/8 | 7/4 | 1 |

Table II. Thermal and percolation thresholds and exponents of Model B.

|  | $`\kappa _{cr}`$ | $`\beta /\nu `$ | $`\gamma /\nu `$ | $`\nu `$ |
| --- | --- | --- | --- | --- |
| Percolation | 0.36673$`{}_{0.00010}{}^{}{}_{}{}^{+0.00012}`$ | 0.528$`{}_{0.015}{}^{}{}_{}{}^{+0.012}`$ | 1.9850$`{}_{0.015}{}^{}{}_{}{}^{+0.010}`$ | 0.632$`{}_{0.015}{}^{}{}_{}{}^{+0.01}`$ |
| Magnetization | 0.36677$`{}_{0.00008}{}^{}{}_{}{}^{+0.00010}`$ | 0.530$`{}_{0.018}{}^{}{}_{}{}^{+0.012}`$ | 1.943$`{}_{0.008}{}^{}{}_{}{}^{+0.019}`$ | 0.640$`{}_{0.018}{}^{}{}_{}{}^{+0.012}`$ |
| Ising Model  |  | 0.518$`{}_{0.007}{}^{}{}_{}{}^{+0.007}`$ | 1.970$`{}_{0.011}{}^{}{}_{}{}^{+0.011}`$ | 0.6289$`{}_{0.0008}{}^{}{}_{}{}^{+0.0008}`$ |

Table III. Thermal and percolation thresholds and exponents of Model C.

|  | $`\kappa _{cr}`$ | $`\beta /\nu `$ | $`\gamma /\nu `$ | $`\nu `$ |
| --- | --- | --- | --- | --- |
| Percolation | 0.3005$`{}_{0.0001}{}^{}{}_{}{}^{+0.0001}`$ | 0.524$`{}_{0.011}{}^{}{}_{}{}^{+0.010}`$ | 1.9750$`{}_{0.009}{}^{}{}_{}{}^{+0.008}`$ | 0.636$`{}_{0.017}{}^{}{}_{}{}^{+0.011}`$ |
| Magnetization | 0.3004$`{}_{0.0001}{}^{}{}_{}{}^{+0.0002}`$ | 0.513$`{}_{0.010}{}^{}{}_{}{}^{+0.012}`$ | 1.963$`{}_{0.009}{}^{}{}_{}{}^{+0.014}`$ | 0.626$`{}_{0.010}{}^{}{}_{}{}^{+0.011}`$ |
| Ising Model  |  | 0.518$`{}_{0.007}{}^{}{}_{}{}^{+0.007}`$ | 1.970$`{}_{0.011}{}^{}{}_{}{}^{+0.011}`$ | 0.6289$`{}_{0.0008}{}^{}{}_{}{}^{+0.0008}`$ |

Table IV. Thermal and percolation thresholds and exponents of Model D.


text:
# 1 Magnitudes of the SM cross-sections and SM plus Randall-Sundrum contributions for initial helicities ++ for a lowest resonance of 600GeV.

UM-P 031-2000

RCHEP 007-2000

Signature of Randall-Sundrum Quantum Gravity model in $`\gamma \gamma `$ scattering in the TeV range.

S.R.Choudhury111src@ducos.ernet.in

Department of Physics, Delhi University, Delhi, India,

A. Cornell222a.cornell@tauon.ph.unimelb.edu.au and G.C. Joshi333joshi@physics.unimelb.edu.au

School of Physics, University of Melbourne,
Parkville, Victoria 3108, Australia

$`22^{nd}`$ of July, 2000

## Abstract

We examine the implications of the Randall-Sundrum gravity models on $`\gamma \gamma `$ scattering in the TeV range.

It has recently been proposed that the scale of weak interactions $`M1TeV`$ is the only fundamental scale of energy in the description of particle interactions. The conventional scale of gravitational interactions, the Planck scale $`M_{pl}M`$ is no longer fundamental but arises out of a Kaluza-Klein compactification of a higher dimensional space-time into four dimensions. In the first such proposals of weak scale quantum gravity (WSQG), one starts with $`(4+n)`$-dimensional space-time with the standard model (SM) fields confined to a three-brane that we live in . Gravity, on the other hand, is postulated to propagate in the full space-time including the $`n`$ extra dimensions which are assumed to be compact of size $`r_c`$. For distances $`rr_c`$, gravity is no different from other SM interactions, but for $`rr_c`$ the net effect of propagation of gravitons in the extra compact dimensions is to make strength of the gravitational interactions from $`1/M`$, via the relation:

$`M_{pl}^2M^{n+2}r_c^2.`$ (1)
For $`M110TeV`$, this last equation gives as low as $`1mm`$, which is the current experimental limit of validity of Newtonian gravity law. This scenario has interesting phenomenological implications in TeV-scale physics  as well as in the precise values of certain experimentally measurable low energy parameters like $`(g2)`$ of the muon . The existence of the $`n`$-compact dimensions of course implies that the graviton is accompanied by a tower of massive cousins. These will be individually gravitationally coupled but added up, these effectively give rise to interactions comparable or even greater than in strength of weak interactions at sufficiently high energies in the TeV scale. Amongst the TeV-scale processes, clearly the ones that will be most sensitive to these ‘new’ interactions are the ones which are sufficiently suppressed in the SM. The scattering process $`\gamma \gamma \gamma \gamma `$ is one such example since in the SM, the amplitude for this process is of $`O(\alpha ^2)`$. This is also a process is which there is a possibility of obtaining TeV-scale data in the near future and thus is worth investigating. Such investigations in the context of the ADD-model just outlined have been done in  and .
An alternative scenario of compactification has recently been proposed by Randall and Sundrum  where, the hierarchy of scales $`M_{pl}`$ and $`M`$ are generated in a different way. The space-time is now a 5D-non factorizable geometry, with a compact angular co-ordinate $`\varphi `$ ranging from $`|\varphi |=0`$ to $`\pi `$. Two 3-branes with opposite tensions reside at the points $`\varphi =0`$ and $`\varphi =\pi `$ and the resultant 4-d metric has the form

$`ds^2=e^{2\sigma (\varphi )}\eta _{\mu \nu }dx^\mu dx^\nu +r_c^2d\varphi ^2,`$ (2)
with $`\mu ,\nu =0,1,2,3`$ and $`\sigma (\varphi )=kr_c|\varphi |`$, $`r_c`$ being the compactification scale. Our universe in this scenario is at $`\varphi =\pi `$ and because of the exponential ‘warp’ factor in the metric, physical masses in our 4-d world are related to the fundamental scale $`M_0`$ of the theory by the relation $`M_{phys}=exp(kr_c\pi )M_0`$. This then is the reason for the weak scale arising with a value of $`kr_c10`$. Further, because of this compact dimension, gravitons will occur as a tower of particles. The zero mode is uniform in $`\varphi `$ and behaves like an ordinary graviton with a 4-d $`M_{pl}`$ scale; the excitations are massive and couple with matter with a scale $`\mathrm{\Lambda }=exp(kr_c)M_{pl}`$ which is of the order of $`1TeV`$. The net Lagrangian of the coupling of this tower of gravitons with matter can then be expressed as a Lagrangian :

$`L=\left({\displaystyle \frac{1}{M_{pl}}}\right)T^{\alpha \beta }(x)h_{\alpha \beta }^{(0)}(x)\left({\displaystyle \frac{1}{\mathrm{\Lambda }}}\right)T^{\alpha \beta }(x){\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}h_{\alpha \beta }^{(n)}(x).`$ (3)
The scenario in the Randall-Sundrum (RS) version of weak scale quantum gravity is thus substantially different from the ADD one, where the excited gravitons effectively form a continuum. In contrast the RS-excitations form a discrete set with well defined and calculable spectrum. We of course have no hint where exactly will the RS-excitations lie in the TeV-scale. However, choosing the mass of the lowest excitation as a parameter, it is possible to calculate the spectrum and this is what was done in the first paper on the phenomenology of the RS-model . We will also consider the mass of the first RS-resonance as a parameter and investigate the behaviour of $`\gamma \gamma `$-scattering amplitude in this scenario. This is thus an extension of the investigations in references by Davoudiasl and by us, where changes in $`\gamma \gamma `$ cross sections arising from weak scale quantum gravity were investigated in the ADD scenario . Should there be indications of departure from the SM in experiments in the near future, the difference between the two available pictures would be useful. This present note addresses thus the question of RS-scenario implication for the $`\gamma \gamma `$ scattering.
We shall assume, that as in the first paper on RS-phenomenology, that the first RS-resonance occurs at a value at 600 GeV; changes in that would only shift the graphs but not the behaviour. Further, arguments based on string theory tells us that the parameter $`k/M_{pl}10^2`$. We will therefore investigate the phenomenology in a range of $`k/M_{pl}`$ around that.
We record first the standard model (SM) amplitudes for $`\gamma \gamma `$ scattering. This is dominated by the W-loop and we take a region of phase space such that $`s`$, $`|t|`$ and $`|u|`$ are much greater than $`M_w^2`$ (that would limit the scattering angle in the c.m. to be not too close to either $`0`$ or $`\pi `$, which of course is experimentally very convenient). Under these approximations, the SM expressions for the helicity amplitudes are, in terms of the standard Mandelstam variables $`s`$, $`t`$, $`u`$:

$`M_{++++}(s,t,u)=\left(i16\pi \alpha ^2\right)\left(\mathrm{log}\left|{\displaystyle \frac{u}{M_w^2}}\right|{\displaystyle \frac{s}{u}}+\mathrm{log}\left|{\displaystyle \frac{t}{M_w^2}}\right|{\displaystyle \frac{s}{t}}\right)`$ (4)
$`M_{++}(s,t,u)`$ $`=`$ $`\left(i12\pi \alpha ^2{\displaystyle \frac{st}{u}}\right)+\left(i{\displaystyle \frac{8\pi \alpha ^2}{u^2}}\right)+\left(4u^23s\right)\left(\mathrm{log}\left|{\displaystyle \frac{t}{u}}\right|\right)`$ (5)
$`\left(i16\pi \alpha ^2\right)\left({\displaystyle \frac{u}{s}}\mathrm{log}\left|{\displaystyle \frac{u}{M_w^2}}\right|+{\displaystyle \frac{u^2}{st}}\mathrm{log}\left|{\displaystyle \frac{t}{M_w^2}}\right|\right)`$
$`M_{++}(s,t,u)=M_{++}(s,t,u)`$ (6)
All other amplitudes, not related to these are negligible in the region considered.
The SM-amplitudes have to be added to the contributions arising out of the exchange of RS-resonances in the $`s`$-, $`t`$\- and $`u`$-channels. If we assume that the first of the resonances is in the sub-TeV range (assumed $`600GeV`$ in our calculation), then up to about $`\sqrt{s}=1TeV`$, it is a good approximation to consider the exchange of the lowest RS-resonance only and we will do that. Using the interaction given in equation (3), the RS-resonance contributions to the amplitudes work out to be:

$`M_{++++}^{RS}`$ $`=`$ $`\left({\displaystyle \frac{s^2}{\mathrm{\Lambda }^2}}\right)\left[D(t)+D(u)\right]`$ (7)
$`M_{++}^{RS}`$ $`=`$ $`\left({\displaystyle \frac{u^2}{\mathrm{\Lambda }^2}}\right)\left[D(s)+D(t)\right]`$ (8)
$`M_{++}^{RS}`$ $`=`$ $`\left({\displaystyle \frac{t^2}{\mathrm{\Lambda }^2}}\right)\left[D(s)+D(u)\right]`$ (9)
where $`D(t)=(tm_1^2)^1`$, $`D(u)=(um_1^2)^1`$
and $`D(s)=(sm_1^2+im_1\mathrm{\Gamma }+\frac{1}{4}\mathrm{\Gamma }^2)^1`$.
$`m_1`$ is the mass of the first RS-resonance and $`\mathrm{\Gamma }`$ its width, which can be calculated as:

$`\mathrm{\Gamma }=\left({\displaystyle \frac{3}{10\pi }}\right)\left(m_1x_1^2\right)\left({\displaystyle \frac{k}{M_{pl}}}\right)^2`$ (10)
$`x_1`$ above is the first zero of the Bessel function $`J_1(x)`$.
The total amplitude $`M^{SM}+M^{RS}`$ is therefore determined in terms of the parameter $`(k/M_{pl})`$ for a given value of $`m_1`$.

Figures 1,2 shows the typical way in which the $`\gamma \gamma `$ cross-sections for parallel and antiparallel initial helicity states behave as a function of the c.m. energy. Also shown are the pure SM-predictions for comparison. To see the nature of the difference between the RS- and ADD-scenario, we have also presented in figures 3 and 4, similar cross-sections for the ADD-scenario with a cut-off value of $`M_s=3TeV`$. We have chosen for the calculations values of $`\sqrt{s}`$ up to $`1TeV`$, hopefully because such energies would be accessible at NLC for $`\gamma \gamma `$ processes using scattering of laser photons from $`e^+e^{}`$ beams.
The most striking feature of the RS-scenario, is the huge enhancement of the cross-sections even for low values of the parameter $`k/M_{pl}`$ that we have considered. In comparison for the pure fermionic process $`e^+e^{}\mu ^+\mu ^{}`$ considered by Davoudiasl, Hewett and Rizzo , the cross section rise is not so spectacular even for a higher range of values of the parameter $`k/M_{pl}`$. A second feature worth mentioning is the presence of resonance peaks only in the helicity state $`+`$ which should be a very clear signature should experiments be possible with polarized photons. The energy dependence of the cross-sections, even in the limited range considered is also very distinctive of the RS-scenario. The SM cross-sections decrease with energy whereas the total (SM+RS) cross-sections increase mildly with energy. Of course, these features would appropriately scaled as the mass of the first resonance takes on different values, but the qualitative nature of these features will continue. It should also be noted that these features are quite different from the cross-sections obtained by considering the ADD version WSQG as can be seen by comparison of figure 3 and 4 with their counterparts, figures 1 and 2.
In conclusion, experimental results on $`\gamma \gamma `$ scattering in TeV scale would not only provide clear indications of possible departure from SM results not merely as a correction factor but by a large change in the value of cross-sections. Not only that, the values determined and their behaviour as a function of energy will enable us to distinguish between the two currently available scenarios of weak scale quantum gravity.


text:
# On the energy equation and efficiency parameter of the common envelope evolution

## 1 Introduction

A very important stage of the evolution of close binaries is the formation of a common envelope. Generally, if a star in a binary fills its Roche-lobe (either because of the swelling up of the envelope due to exhaustion of nuclear fuel at the centre – or reduction in the separation due to loss of orbital angular momentum) there will be an attempt to transfer mass to the companion. However, if the thermal re-adjustment timescale of the accreting star is larger than the mass-transfer timescale, the accreted layer piles up above the companion, heats up and expands so that the accreting star will also fill its Roche-lobe. The result is that the transfered matter will form a common envelope embedding both stars (Paczynski 1976; Ostriker 1976). When the mass ratio is large the CE-phase is accompanied by the creation of a drag-force, arising from the motion of the companion star through the envelope of the evolved star, which leads to dissipation of orbital angular momentum (spiral-in process) and deposition of orbital energy in the envelope. Hence, the global outcome of a CE-phase is reduction of the binary separation and often ejection of the envelope. For a general review on common envelopes, see e.g. Iben & Livio (1993).

There is clear evidence of orbital shrinkage (as brought about by frictional torques in a CE-phase) in some observed close binary pulsars and white dwarf binaries (e.g. PSR 1913+16, L 870-2 – see Iben & Livio 1993 and references therein). In these systems it is clear that the precursor of the last-formed degenerate star must have achieved a radius much larger than the current orbital separation. Since the present short separation of degenerate stars in some binaries can not be explained by gravitational wave radiation or a magnetic wind as the main source of the loss of $`J_{\mathrm{orb}}`$, frictional angular momentum loss was probably responsible for the very close orbits observed. An alternative process recently investigated (King & Begelman 1999; Tauris et al. 2000) is a highly Super-Eddington mass transfer on a subthermal timescale. However, as shown by the latter authors this process can not explain the observed properties of all systems. They conclude that CE and spiral-in evolution is still needed to account for the formation of double neutron star systems and mildly recycled pulsars with a heavy white dwarf companion and $`P_{\mathrm{orb}}<3`$ days.

There are many uncertainties involved in calculations of the spiral-in phase during the CE evolution. The evolution is tidally unstable and the angular momentum transfer, dissipation of orbital energy and structural changes of the donor star take place on very short timescales ($`10^310^4`$ yr). A complete study of the problem would also require very detailed multi-dimensional hydrodynamical calculations which is beyond the scope of this paper. Here we aim to improve the simple formalism developed by Webbink (1984) and de Kool (1990) for estimating the orbital evolution of binaries evolving through a CE and spiral-in phase. An important quantity in their formalism is the binding energy between the envelope and the core of the donor star. This quantity is expressed by a parameter $`\lambda `$ which depends on the stellar density distribution, and consequently also on the evolutionary stage of the star. With a detailed evolutionary calculation, Bisscheroux (1998) found that $`\lambda 0.40.6`$ on the red and asymptotic giant branches of a 5 $`\mathrm{M}_{\mathrm{}}`$ star. In the literature this parameter is usually taken to be a constant (e.g. $`\lambda =0.5`$) disregarding the evolutionary stage (structure) of the donor star at the onset of the mass transfer process. This in turn leads to a misleading (over)estimation of the so-called efficiency parameter, $`\eta _{\mathrm{CE}}`$, often discussed by various authors. From observations one is able to put constraints on the product $`\eta _{\mathrm{CE}}\lambda `$ (see eq. 4 below). However, by underestimating $`\lambda `$ the derived values for $`\eta _{\mathrm{CE}}`$ will be overestimated – e.g. van den Heuvel (1994) and Tauris (1996) concluded $`\eta _{\mathrm{CE}}>1`$ by assuming $`\lambda =0.5`$. By calculating $`\lambda `$ as a real function of stellar radius we are able to improve the estimates of post-CE orbital separations and hence put more realistic constraints on $`\eta _{\mathrm{CE}}`$.

Detailed calculations on the binding energy of the envelope to determine the outcome of binaries evolving through a CE-phase have been done by Han et al. (1994, 1995). The method and the stellar evolution code used to calculate the binding energy in their work is similar to the one applied in our work. However, here we discuss the binding energy within the context of the $`\lambda `$-parameter in the Webbink-formalism. We present calculated values of $`\lambda `$ for stars with initial mass in the range of 3 – 10 $`\mathrm{M}_{\mathrm{}}`$ at different evolutionary stages. Such values are usefull for an estimation of the orbital decay during a common envelope and spiral-in phase. We also demonstrate a new approach to find a unique solution for the stellar and binary parameters of systems evolving through a CE and spiral-in phase.

We will in particular consider binaries in which the companion of the evolved star is a neutron star. Such binaries are the progenitors of binary millisecond pulsars (BMSPs) with a massive degenerate companion (e.g. van den Heuvel 1994) and they may appear in nature as short lived intermediate-mass X-ray binaries (IMXBs).

In Sect. 2 we describe the basic concepts of common envelope evolution concerning the energy equation and the problems of defining the binding energy. Sect. 3 introduces briefly the numerical stellar code used to evaluate the parameter $`\lambda `$. The results are presented in Sect. 4 and discussed in Sect. 5. We summarize our conclusions in Sect. 6.

## 2 Common envelope and spiral-in phase

### 2.1 The energy equation and orbital evolution

A simple estimation on the reduction of the orbit can be found by simply equating the difference in orbital energy (before and after the CE-phase) to the binding energy of the envelope of the (sub)giant donor. Following the formalism of Webbink (1984), the original binding energy of the envelope at the onset of mass transfer will be of the order: $`GM_{\mathrm{donor}}M_{\mathrm{env}}/a_\mathrm{i}r_\mathrm{L}`$, where $`M_{\mathrm{donor}}`$ is the mass of the donor star at the beginning of the CE-phase; $`M_{\mathrm{env}}`$ is the mass of the hydrogen-rich envelope of the donor; $`a_\mathrm{i}`$ is the orbital separation at the onset of the CE and $`r_\mathrm{L}=R_\mathrm{L}/a_\mathrm{i}`$ is the dimensionless Roche-lobe radius of the donor star, so that $`a_\mathrm{i}r_\mathrm{L}=R_\mathrm{L}R_{\mathrm{donor}}`$. A parameter $`\lambda `$ was introduced by de Kool (1990) as a numerical factor (of order unity) which depends on the stellar density distribution, such that the real binding energy of the envelope can be expressed as:

$`E_{\mathrm{env}}`$ $`=`$ $`{\displaystyle \frac{GM_{\mathrm{donor}}M_{\mathrm{env}}}{\lambda a_\mathrm{i}r_\mathrm{L}}}`$ (1)
The total change in orbital energy is given by:

$`\mathrm{\Delta }E_{\mathrm{orb}}`$ $`=`$ $`{\displaystyle \frac{GM_{\mathrm{core}}M_2}{2a_\mathrm{f}}}+{\displaystyle \frac{GM_{\mathrm{donor}}M_2}{2a_\mathrm{i}}}`$ (2)
where $`M_{\mathrm{core}}=M_{\mathrm{donor}}M_{\mathrm{env}}`$ is the mass of the helium core of the evolved donor star; $`M_2`$ is the mass of the companion star and $`a_\mathrm{f}`$ is the final orbital separation after the CE-phase.

Let $`\eta _{\mathrm{CE}}`$ describe the efficiency of ejecting the envelope, i.e. of converting orbital energy into the kinetic energy that provides the outward motion of the envelope: $`E_{\mathrm{env}}\eta _{\mathrm{CE}}\mathrm{\Delta }E_{\mathrm{orb}}`$ or, by equating eqs. (1) and (2):

$`{\displaystyle \frac{GM_{\mathrm{donor}}M_{\mathrm{env}}}{\lambda a_\mathrm{i}r_\mathrm{L}}}`$ $`=`$ $`\eta _{\mathrm{CE}}\left[{\displaystyle \frac{GM_{\mathrm{core}}M_2}{2a_\mathrm{f}}}{\displaystyle \frac{GM_{\mathrm{donor}}M_2}{2a_\mathrm{i}}}\right]`$ (3)
which yields:

$`{\displaystyle \frac{a_\mathrm{f}}{a_\mathrm{i}}}`$ $`=`$ $`{\displaystyle \frac{M_{\mathrm{core}}M_2}{M_{\mathrm{donor}}}}{\displaystyle \frac{1}{M_2+2M_{\mathrm{env}}/(\eta _{\mathrm{CE}}\lambda r_\mathrm{L})}}`$ (4)
### 2.2 The binding energy of the envelope

The total binding energy of the envelope to the core is given by:

$`E_{\mathrm{bind}}`$ $`=`$ $`{\displaystyle _{M_{\mathrm{core}}}^{M_{\mathrm{donor}}}}\left({\displaystyle \frac{GM(r)}{r}}+U\right)𝑑m`$ (5)
where the first term is the gravitational binding energy and $`U`$ is the internal thermodynamic energy. The latter involves the basic thermal energy for a simple perfect gas ($`3\mathrm{}T/2\mu `$), the energy of radiation ($`aT^4/3\rho `$), as well as terms due to ionization of atoms and dissociation of molecules and the Fermi energy of a degenerate electron gas (Han et al. 1994, 1995).

The binding energy of the envelope can be expressed by eq. (5), under the assumption that the entire internal energy is used efficiently in the ejection process. However, it is not known how much of the internal energy is involved to unbind the envelope. Han et al. (1995) introduced a parameter $`\alpha _{\mathrm{th}}`$ and expressed the envelope binding energy as:

$`E_{\mathrm{env}}`$ $`=`$ $`{\displaystyle _{M_{\mathrm{core}}}^{M_{\mathrm{donor}}}}{\displaystyle \frac{GM(r)}{r}}𝑑m+\alpha _{\mathrm{th}}{\displaystyle _{M_{\mathrm{core}}}^{M_{\mathrm{donor}}}}U𝑑m`$ (6)
The value of $`\alpha _{\mathrm{th}}`$ depends on the details of the ejection process, which is very uncertain. A value of $`\alpha _{\mathrm{th}}`$ equal to 0 or 1 corresponds to maximum and minimum envelope binding energy, respectively. By simply equating eqs. (1) and (6) we calculated the parameter $`\lambda `$ for different evolutionary stages of a given star. The minimum and maximum derived values of $`\lambda `$ are denoted by $`\lambda _\mathrm{g}`$ ($`\alpha _{\mathrm{th}}`$ = 0) and $`\lambda _\mathrm{b}`$ ($`\alpha _{\mathrm{th}}`$ = 1).

By taking the core mass $`M_{\mathrm{core}}`$ as the lower boundary of the integral, we assume that the structure of the core does not change during the envelope ejection process, and therefore we assume that there is no exchange of energy between the core and the envelope. Here, we define the core mass as the central mass which contains less than 10 % hydrogen – see discussion in Sect. 5.2.

For the case of binary stars we use the radii of the stars to calculate the separations at the onset of the mass transfer assuming a companion star of mass $`M_2`$. An estimation on the Roche-lobe radius is given by Eggleton (1983):

$`R_\mathrm{L}`$ $`=`$ $`{\displaystyle \frac{0.49a_\mathrm{i}}{0.6+q^{2/3}\mathrm{ln}(1+q^{1/3})}}`$ (7)
where $`q=M_{\mathrm{donor}}/M_2`$ is the mass ratio.

## 3 A brief introduction to the numerical computer code

We used an updated version of the numerical stellar evolution code of Eggleton (1971, 1972, 1973). This code uses a self-adaptive, non-Lagrangian mesh-spacing which is a function of local pressure, temperature, Lagrangian mass and radius. It treats both convective and semi-convective mixing as a diffusion process and finds a simultaneous and implicit solution of both the stellar structure equations and the diffusion equations for the chemical composition. New improvements are the inclusion of pressure ionization and Coulomb interactions in the equation-of-state, and the incorporation of recent opacity tables, nuclear reaction rates and neutrino loss rates. The most important recent updates of this code are described in Pols et al. (1995, 1998) and some are explained in Han et al. (1994).

We performed such detailed numerical stellar evolution calculations in our work since they should yield more realistic results compared to models based on complete, composite or condensed polytropes.

## 4 Results

The above-mentioned detailed stellar evolution code has been used to probe the structure of donor stars at the onset of mass transfer. We evolved stars with an initial mass in the interval of 3 – 10 $`\mathrm{M}_{\mathrm{}}`$ assuming a chemical composition of ($`X=0.70,Z=0.02`$) and using a mixing-length parameter of $`\alpha =l/H_p=2.0`$. Convective overshooting is taken into account in the same way as in Pols et al. (1998) using an overshooting constant $`\delta _{\mathrm{ov}}=0.1`$. We used de Jager’s wind mass-loss rate (de Jager et al. 1988; Nieuwenhuijzen & de Jager 1990) to estimate the mass loss prior to the Roche-lobe overflow and CE evolution. The wind mass loss from the donor star (and subsequent widening of the orbit) prior to the Roche-lobe overflow is typically only $`14`$ % on the AGB and hence this effect is not very significant for the stars considered here.

With the Eggleton code, the total binding energy is computed for each mesh-point. Included here are the ionization of $`\mathrm{H}^+`$, H, $`\mathrm{He}^{++}`$, $`\mathrm{He}^+`$, He, dissociation of $`\mathrm{H}_2`$, rotational and vibrational modes of $`\mathrm{H}_2`$, as well as ionization of seven heavier elements (C, N, O, Ne, Mg, Si, Fe) which are assumed to be fully ionized at all temperatures and densities (Pols et al. 1995). By integrating from the core to the photosphere, we found the total binding energy of the envelope. We used this energy to derive the parameter $`\lambda `$. In Fig. 1 we present the relation between $`\lambda `$ and stellar radius for each of the different donor stars considered. Here $`\lambda `$ was determined either from the total binding energy of the envelope (i.e. $`\alpha _{\mathrm{th}}=1`$) or from the gravitational binding energy alone (i.e. $`\alpha _{\mathrm{th}}=0`$).

For a donor star in a binary system the stellar radius at the onset of the CE-phase is roughly equivalent to its Roche-lobe radius, and given a mass of its companion star one finds the corresponding orbital period of the system. In Table 1, assuming that the companion is a 1.3 $`\mathrm{M}_{\mathrm{}}`$ neutron star, we list the relation between the orbital period of the binary and $`\lambda `$ at the onset of the mass transfer (again due to the wind mass loss, this period might be slightly larger than the initial ZAMS period). We distinguish $`\lambda _\mathrm{g}`$ (as derived from gravitational binding energy alone) from $`\lambda _\mathrm{b}`$ (calculated from the total binding energy). The gyration radius of the star, $`k^2=I/MR^2`$, is also presented to provide information about the orbital tidal instability. Mass transfer becomes tidally unstable if the rotational angular momentum of the donor exceeds one third of the orbital angular momentum (Darwin 1908).

## 5 Discussion

### 5.1 On the binding energy of the envelope

In Fig. 2 we have plotted the binding energy of the envelope as well as the density profile of two stars with an initial mass of 4 and 10 $`\mathrm{M}_{\mathrm{}}`$. When the star reaches the asymptotic giant branch, the core is surrounded by a convective envelope, and the density gradient around the core mass is very steep. As the star ascends this giant branch, the convective envelope becomes deeper and the density contrast higher.

Assuming a star to be a non-rotating globe of non-degenerate ionized hydrogen and helium in hydrostatic equilibrium (i.e. nearly an ideal monatomic gas) one would expect from the virial theorem: $`\lambda _\mathrm{b}2\lambda _\mathrm{g}`$. However at late stellar evolutionary stages this picture no longer holds. After the second dredge-up stage of evolution, the gravitational binding energy becomes smaller as the stellar radius expands to a supergiant dimension. Meanwhile the ionization and dissociation energy becomes more significant when the temperature decreases in the envelope (Bisscheroux 1998) – raising the internal energy above 0.5 times the gravitational binding energy in absolute value. Hence, the total binding energy becomes smaller in absolute value relative to the gravitational binding energy. This in turn causes a larger discrepancy between $`\lambda _\mathrm{g}`$ and $`\lambda _\mathrm{b}`$, see upper right panel in Fig. 1.

Especially for low-mass stars (3 – 6 $`\mathrm{M}_{\mathrm{}}`$) we see that the contribution of the internal energy to the total binding energy becomes dominating. At some point on the asymptotic giant branch, the internal energy dominates over the gravitational binding energy, and hence the total binding energy of the envelope reaches a positive value. One might regard this alteration from a negative to a positive value of the total binding energy as the ejection of the giant’s envelope, or as a point of the onset of a pulsationally driven (Mira) superwind (Han et al. 1994 and references therein). However, when the binding energy becomes positive, the derived parameter $`\lambda _\mathrm{b}`$ will be negative and our formalism in Sect. 2.1 breaks down – cf. resulting ”negative” final separation in eq. (4). On the contrary $`\lambda _\mathrm{g}`$ is relatively constant during this late stage of evolution.

The evolution of stars above 7 $`\mathrm{M}_{\mathrm{}}`$ is cut short by carbon burning. Therefore they never move as far up the AGB as lower mass stars, and the internal (recombination) energy never becomes so important as to dominate the total binding energy. This is why the difference between $`\lambda _\mathrm{g}`$ and $`\lambda _\mathrm{b}`$ is less pronounced as we go to the more massive stars.

### 5.2 On the definition of the core mass boundary

Our results are sensitive to the definition of core mass. In our work we have chosen to define the core mass as the central mass which contains less than 10 % hydrogen. Using other boundary conditions such as $`^2\mathrm{log}\rho /m^2=0`$ (e.g. Bisscheroux 1998) or the alternative method by Han et al. (1994) will lead to different values for the binding energy of the envelope and hence different estimates of $`\lambda `$. In general the above mentioned methods will result in a core mass boundary further out in the star compared to our definition (see symbols on the dashed line in Fig. 2). Hence their methods lead to smaller envelope binding energies and thus larger values of $`\lambda `$ (factor $`2`$) on the RGB and EAGB. For a 10 $`\mathrm{M}_{\mathrm{}}`$ star, this difference is much higher than a factor of 2. Further up the AGB, where the density gradient is steeper, the differences become smaller and the result more reliable. Nevertheless, in accordance with our findings using any of the two other boundary conditions mentioned will also result in large values of $`\lambda `$ on the AGB (for low-mass stars: $`\lambda _\mathrm{b}1`$) and this is the important issue.

### 5.3 Dependency of the $`\lambda `$-parameter on the chemical composition and mixing-length parameter

In order to investigate the dependency of the $`\lambda `$-parameter on chemical composition and mixing-length parameter, we have plotted in Fig. 3 the computed $`\lambda _\mathrm{b}(R)`$ for both Pop. I and Pop. II chemical abundancies, as well as for two different values of the mixing-length parameter for a 6 $`\mathrm{M}_{\mathrm{}}`$ donor star.

¿From the ZAMS to the base of the RGB ($`R50`$ $`\mathrm{M}_{\mathrm{}}`$) there are no differences on the calculated values of $`\lambda _\mathrm{b}`$. It is also seen that there is hardly any dependency on the chemical composition for the rest of the evolutionary track far up the AGB. However, from the horizontal branch ($`R150`$ $`\mathrm{M}_{\mathrm{}}`$) and halfway up the AGB ($`R300`$ $`\mathrm{M}_{\mathrm{}}`$) a higher value of $`\lambda _\mathrm{b}`$ is found using a larger value of the mixing-length parameter. This picture is reversed near the top of the AGB ($`R>400`$ $`\mathrm{M}_{\mathrm{}}`$). In general we can conclude that $`\lambda _\mathrm{b}`$ is almost independent of the chemical composition of the star, but it depends on the choice of the mixing-length parameter.

### 5.4 The efficiency parameter of the ejection process

In order to explain the long orbital periods observed among those of the BMSPs which are likely to have evolved through a CE-phase, Tauris (1996) and van den Heuvel (1994) found it necessary to invoke efficiency parameters larger than unity (e.g. $`1<\eta _{\mathrm{CE}}<4`$). However, their simple estimates were based on a constant value of $`\lambda =0.5`$. Here we suggest that a high value of $`\lambda `$ (rather than a high value of $`\eta _{\mathrm{CE}}`$) is the natural explanation for the long orbital periods observed. This in turn indicates that if the internal energy is used efficiently ($`\alpha _{\mathrm{th}}=1`$) to unbind the envelope, there is no need to assume $`\eta _{\mathrm{CE}}>1`$ for these systems. This result was also concluded by Han et al. (1995).

Let us now consider a system with a 9 $`\mathrm{M}_{\mathrm{}}`$ donor star and a 1.3 $`\mathrm{M}_{\mathrm{}}`$ neutron star. Such a binary is the progenitor of a pulsar with a massive degenerate companion. We computed the final orbital separation of the binary after a CE-phase by using $`E_{\mathrm{env}}=\eta _{\mathrm{CE}}\mathrm{\Delta }E_{\mathrm{orb}}`$ combined with eqs. (2) and (6). Here we applied $`\eta _{\mathrm{CE}}=1`$ and $`\alpha _{\mathrm{th}}=1`$ (i.e. including the internal thermodynamic energy, $`\lambda =\lambda _\mathrm{b}`$), and present the result in the upper panel of Fig. 4. As a comparison, we also show the final separation computed by using eq. (4) with a constant efficiency parameter $`\lambda =0.5`$ (commonly used in the literature).

A binary will survive the CE-phase only if the final orbit is large enough so that the core of the donor will not instantly fill its Roche-lobe. In the case of a binary with a 9 $`\mathrm{M}_{\mathrm{}}`$ donor and a 1.3 $`\mathrm{M}_{\mathrm{}}`$ neutron star, this condition is reached if the binary separation at the onset of the CE-phase is larger than 435 $`\mathrm{R}_{\mathrm{}}`$ ($`P_{\mathrm{orb}}>327^\mathrm{d}`$). If we use the constant $`\lambda =0.5`$, this separation is reduced to 175 $`\mathrm{R}_{\mathrm{}}`$ ($`P_{\mathrm{orb}}>84^\mathrm{d}`$). In the latter case the donor star still has a radiative envelope with only a very thin convective shell at the onset of the mass transfer. Therefore the CE is formed only as a result of orbital shrinkage due to mass transfer with a high mass ratio between the donor and the neutron star. In the former case, the donor star had time to develop a deep convective envelope and it will therefore expand in response to mass loss due to its isentropic entropy profile. This expansion, in combination with the orbital shrinkage mentioned above, will immediately lead to a runaway mass transfer, and hence a CE will be formed easily.
At late stages of stellar evolution, the value of $`\lambda _\mathrm{b}`$ has an important effect on the final orbital period. This effect is more significant in binaries with less massive donors as can be seen in the lower panel of Fig. 4 for the case of a 5 $`\mathrm{M}_{\mathrm{}}`$ donor.

### 5.5 A new approach for tracking the evolutionary history of a system evolving through a CE

As we have demonstrated in this paper, $`\lambda `$ is a function of stellar radius for any given mass. Combining this fact with eq. (4), which relates the observed orbital period to the original radius of the Roche-lobe filling donor star as a function of mass and $`\lambda `$ (or equivalently: $`a_\mathrm{f}/a_\mathrm{i}`$), we are able to place severe constraints on the mass and radius (and hence age) of the original donor.

An example of this approach is shown in Fig. 5 where we have used the parameters for the recently discovered binary pulsar PSR J1454–5846 (Camilo et al. 2000). This system has an orbital period of $`P_{\mathrm{orb}}=12\stackrel{d}{.}42`$ and a mass function, $`f=0.130`$ $`\mathrm{M}_{\mathrm{}}`$. The exact mass of the white dwarf in this system is unknown; but as an example we shall assume an orbital inclination angle, $`i=60\mathrm{°}`$ yielding $`M_{\mathrm{WD}}=1.03`$ $`\mathrm{M}_{\mathrm{}}`$(for $`M_{\mathrm{NS}}=1.30`$ $`\mathrm{M}_{\mathrm{}}`$). Hence we must require that the original giant donor star had a core mass of this value by the time its envelope was ejected in the CE-phase. In Fig. 5 the intersection points for a 5 $`\mathrm{M}_{\mathrm{}}`$ and 6 $`\mathrm{M}_{\mathrm{}}`$ donor yield core mass values of 0.93 $`\mathrm{M}_{\mathrm{}}`$ and 1.10 $`\mathrm{M}_{\mathrm{}}`$, respectively. To obtain a core mass of 1.03 $`\mathrm{M}_{\mathrm{}}`$ we must therefore search for a solution in-between. We find $`M_{\mathrm{donor}}=5.6`$ $`\mathrm{M}_{\mathrm{}}`$ and the radius at the onset of the RLO is 340 $`\mathrm{R}_{\mathrm{}}`$. The corresponding age of the donor is 93 Myr and the pre-CE orbital period is 760<sup>d</sup>.

The above solution was found for an inclination angle, $`i=60\mathrm{°}`$. In Fig. 6 we have plotted the solutions for $`M_{\mathrm{donor}}`$ and $`R_\mathrm{L}`$ as a function of inclination angle for two different values of the efficiency parameter, $`\eta _{\mathrm{CE}}`$. It is interesting to notice that $`M_{\mathrm{donor}}`$ is more or less independent of $`\eta _{\mathrm{CE}}`$ whereas $`R_\mathrm{L}`$ gets larger for smaller values of $`\eta _{\mathrm{CE}}`$. This result is expected theoretically since a low value of $`\eta _{\mathrm{CE}}`$ requires either more orbital energy released and/or less binding energy of the envelope, and thus a larger $`R_\mathrm{L}`$, in order to eject the envelope successfully. An upper limit to the mass of the white dwarf is about 1.3 $`\mathrm{M}_{\mathrm{}}`$ which yields $`i>47\mathrm{°}`$. For $`\eta _{\mathrm{CE}}=0.5`$ no solution is found for $`i<59\mathrm{°}`$. Hence it would be interesting if a future detection of a spectroscopic line from the white dwarf would yield an orbital inclination angle less than this value, since small values of $`\eta _{\mathrm{CE}}`$ would be excluded that way.

It should also be noted from Fig. 6 that the binary pulsar must have had an AGB radius at the onset of the RLO. Many of the observed binary pulsars which have survived a common envelope evolution also have long orbital periods. Hence, the density gradient was fairly steep at the onset of the Roche-lobe overflow and therefore the uncertainties regarding the core mass definition as already discussed in Sect. 5.2 are less important in this context.

## 6 Conclusions

We have adapted a numerical computer code to study the stellar structure of stars with an initial mass of 3 – 10 $`\mathrm{M}_{\mathrm{}}`$ in order to evaluate the binding energy of the envelope to the core which determines the parameter $`\lambda `$ in the energy equation of the common envelope evolution. We have presented evidence that the parameter $`\lambda `$ depends strongly on the evolutionary stage. Hence the value of $`\lambda `$, at the onset of the mass transfer and CE-phase, is unique for binaries with different initial orbital periods and donor masses. Taking advantage of this, we have demonstrated a new approach for finding a unique solution to the energy equation and the original mass and radius of the donor star. An application of this approach on PSR J1454–5846 yields a constraint for the donor star and the pre-CE orbital period as a function of the unknown orbital inclination angle.

We have demonstrated that in order to obtain a long final orbital period after a CE-phase (as observed in some BMSPs), we do not require an efficiency parameter larger than unity, if the internal energy can be used efficiently to eject the envelope.

In order to obtain the final orbital separation of a binary system after a CE and spiral-in phase, we advise one to calculate the binding energy from the stellar structure by means of eq. (6) with $`0\alpha _{\mathrm{th}}1`$. In case the detailed structure of the donor is not available (e.g. for a quick back-of-the-envelope-calculation), eq. (4) can be used with $`\eta _{\mathrm{CE}}=1`$ in combination with $`\lambda _\mathrm{g}\lambda \lambda _\mathrm{b}`$. The value of $`\lambda `$ can be found in our Table 1 for a given donor star mass and radius.

###### Acknowledgements.

This work was sponsored by NWO Spinoza Grant 08-0 to E. P. J. van den Heuvel. The authors thank E. P. J. van den Heuvel, Gertjan Savonije, Gijs Nelemans and Winardi Sutantyo for stimulating and fruitful discussions; and Onno Pols for his valuable referee suggestions and comments. T. M. T. acknowledges the receipt of a NORDITA fellowship.