text:
# Continuous stellar mass-loss in N-body models of galaxies

## 1 Introduction

During the last decade, many N-body simulations (e.g. Friedli & Benz 1993; Junqueira & Combes 1996) as well as theoretical works on two-fluid gravitational instabilities (Jog 1992, 1996) have underlined the importance of the gas mass fraction for the dynamics and evolution of galaxies.

Gas crucially conditions star formation and evolution of large-scale instabilities (e.g. bars, spiral arms). In turn, star formation and large-scale flows induced by these instabilities influence mass and chemical profiles of galaxies (e.g. Martinet & Friedli 1997; Portinari & Chiosi 2000).

The gas content in a given area is determined by the competing processes of star formation and stellar mass-loss and by spatial flows of gas. A coherent picture of galactic evolution must therefore consistently couple stellar and gaseous dynamics with star formation and stellar evolution.

Star formation is, in large-scale models of galaxies, typically implemented through simplistic parametrizations reflecting mainly the energy supply from young stars and the gas mass fraction locked in newly born stars. On the other hand, stellar mass-loss has received much less attention from galactic N-body modellers since it was for a long time considered as a secondary issue. But nowadays, both observations and stellar evolutionary models (see Sect. 2 for references) indicate that the gas mass fraction restituted by stars is huge and may reach, when integrated over the stellar mass spectrum, some 45% over the Hubble time. Stars thus represent not only a place for permanent gas blocking but also an important temporary reservoir of gas that will be gradually reinjected into the interstellar medium.

Recent computer models of galactic dynamics usually approximate the stellar mass-loss as instantaneous, i.e. happening at the moment of stellar birth. Nonetheless, stars lose matter during all their lives. Stellar lifetimes span a very large range and in the case of low-mass stars they compare with or overpass the Hubble time. The use of the instantaneous recycling approximation, which was proposed by Tinsley (1980) for high mass-stars, is therefore not satisfactory for the whole stellar mass spectrum.

From the above emerges an obvious motivation for building a computer code able to follow the galactic structure and dynamics together with a non-instantaneous gas recycling. One can expect that the connection of such a recycling with the dynamics will affect large-scale gravitational instabilities, spatial flows of matter, star formation and gas consumption rates, etc., and thus the long-term evolution of galaxies.

The development and presentation of such a code is the central aim of this paper. Our model is innovative especially in introducing a continuous (i.e. time-dependent) gas recycling scheme, grafted on an underlying N-body code using standard techniques for computing gravitational field (particle-mesh method) and gas dynamics (sticky-particles).

The text is organized as follows. In Sect. 2, we describe data and steps in deriving the time-dependence of mass-loss of a computer stellar particle. This curve, simply parametrized, is the key input information for our N-body code, however it can be useful also for other models (not necessarily N-body) of galaxies or star clusters. Section 3 presents our computer code, the emphasis being put on the implementation of the mass-loss. An instructive simulation of a barred galaxy is shown in Sect. 4. Section 5 summarizes the outcome and future prospects.

## 2 Stellar mass-loss and long-term galactic evolution

### 2.1 Mass-loss rate of a coeval stellar population

Stars of all masses lose matter via stellar winds, high-mass stars, in addition, in supernova explosions. Despite the fact that stellar winds accompany all stages of stellar evolution (see de Jager et al. 1988, for observations of wind rates across the Hertzsprung-Russell diagram), they vary in intensity and duration: gas is not released uniformly over a stellar lifetime but typically in a few relatively short episodes.

The time dependence of mass-loss rate of a coeval stellar population (hereafter CSP), $`\dot{M}(t)`$, could in principle be obtained if we knew the stellar initial mass function (IMF), $`\psi (m_\mathrm{i})`$, and the time dependence of the mass-loss rate of individual stars, $`\dot{m}`$:

$$\dot{M}\left(t\right)=_{m_{\mathrm{min}}}^{m_{\mathrm{max}}}\dot{m}(t,m_\mathrm{i})\psi \left(m_\mathrm{i}\right)dm_\mathrm{i}$$
$`\left(1\right)`$
Since neither the former nor the latter are known with much precision, we shall limit ourselves to construct a model curve based on several simplifying assumptions adequate to the purpose of the work which is to simulate the effects of stellar mass-loss on the long-term evolution of galaxies; we do not pretend to model changes occurring on a time-scale of less than $`10`$ Myr that is short compared to the dynamical time at most galactocentric distances.

Mass-loss for stars of different masses end evolutionary stages is reviewed in more detail in Jungwiert (1998) where references can be found for both observations and models. Here, only a brief summary of the most important mass-loss mechanisms is given: a) in the case of low-mass stars (LMS, hereafter defined as stars with initial masses $`m_\mathrm{i}<2\mathrm{M}_{}`$), winds on the red giant branch (RGB) and asymptotic giant branch (AGB) are the most efficient; b) for intermediate-mass stars (IMS, $`2\mathrm{M}_{}<m_\mathrm{i}<8\mathrm{M}_{}`$), the AGB winds dominate; c) in the case of high-mass stars (HMS, $`m_\mathrm{i}>8\mathrm{M}_{}`$), supernova explosions, Wolf-Rayet winds and main-sequence winds, contribute to the overall mass-loss in proportions strongly depending on $`m_\mathrm{i}`$: mass-loss is dominated by supernovae for $`m_\mathrm{i}<25\mathrm{M}_{}`$, while for more massive stars the winds are more efficient than supernovae.

With regard to the above, we shall assume, for simplicity, that: a) for LMS all the mass-loss takes place in two distinct delta function events corresponding to the tips of RGB and AGB (since the respective winds strongly peak near these tips); b) for IMS and HMS, we treat the mass-loss as one delta function at the end of their lives: this reflects the AGB wind peak for intermediate-mass stars and SN explosion + winds for high-mass stars (as the most massive HMS, for which the winds dominate over supernovae, have very short lifetimes – less than 5 Myr for $`m_\mathrm{i}>40\mathrm{M}_{}`$ – the delta function representation is adequate).

Necessary ingredients to construct an approximate mass-loss curve of a CSP thus reduce to the knowledge of IMF, stellar initial-final mass relation and stellar lifetime-mass relation (in the case of LMS, two lifetimes are necessary, for the ends of the RGB and AGB phases). As for the IMF, we use the power-law with three slopes (0.2, 1.7 and 1.3 for initial mass ranges of 0.1-1 $`\mathrm{M}_{}`$, 1-10 $`\mathrm{M}_{}`$ and 10-100 $`\mathrm{M}_{}`$, respectively) inferred by Scalo (1998) as the “average” of recently published IMFs. For the initial-final mass relation and lifetime-mass relation, we rely on the Padua stellar evolutionary models (Bressan et al. 1993; Marigo et al. 1996) computed for the initial metallicity $`Z=0.02`$.

The mass-loss rate function, $`\dot{M}_\mathrm{n}(t)`$, constructed in this way (the subscript “n” denotes the normalization to the initial mass, $`M_\mathrm{i}`$, of a CSP) is shown in Fig. 1 (left). For the purpose of its efficient implementation into the N-body code, we carry out a least-squares fit. The power-law fit, $`\dot{M}_\mathrm{n}(t)t^y`$ suggests a very simple form since it gives $`y=1.04`$, very close to unity. Rather than using 1.04, we fix $`y=1`$, and perform a new least-squares fit with a hyperbola having a shift $`T_0`$ from the moment $`t_{\mathrm{birth}}`$, at which a CSP was born, to avoid the unphysical singularity of the $`t^1`$ function at $`t_{\mathrm{birth}}`$:

$$\dot{M}_\mathrm{n}\left(t\right)\frac{\dot{M}\left(t\right)}{M_\mathrm{i}}=\frac{c_0}{tt_{\mathrm{birth}}+T_0}.$$
$`\left(2\right)`$
We restrict the least-squares fits to combinations ($`c_0`$, $`T_0`$) giving the correct total gas return (cumulative mass-loss) $`R_\tau `$ over time $`\tau `$ (see Sect. 2.2. and Table 2), corresponding to the time interval between the birth of a given CSP and the envisaged end of an N-body simulation. Results of the fit for $`\tau =`$ 5, 10, 15 and 20 Gyr are summarized in Table 1.

The fit for the simulation length of 20 Gyr is shown in Fig. 1 (left) by a dashed line. It is close at all times and over four orders of magnitude on the $`\dot{M}_\mathrm{n}(t)`$ axis to the original dependence. It is also obvious, from Table 1, that the fits for different $`\tau `$ are very similar. We will take advantage of this when implementing the mass-loss into the N-body code: only one fit will be used for all the particles (see Sects. 3.4. and 3.5.).

### 2.2 Cumulative mass-loss of a coeval stellar population

While the time-dependence of mass-loss rate of a CSP is a vital input for our models, the cumulative mass-loss, i.e. the mass fraction lost by stars until time $`t`$ after the birth of the CSP,

$$R_t=_0^t\dot{M}_\mathrm{n}(t^{})dt^{},$$
$`(3)`$
gives a useful estimate of the importance of gas recycling on different time-scales.

Fig. 1 (right) shows, in the lower half of the plot, $`R_t`$ of a CSP computed for the Scalo’s (1998) IMF and the Padua $`Z=0.02`$ stellar models. The curve in the upper half of the plot gives the fractional mass remaining in stars, i.e. $`1R_t`$. Table 2 gives $`R_t`$ for $`t=`$ 0.1, 1, 5, 10, 15 and 20 Gyr.

A discussion of mass-loss for other IMFs and stellar metallicities is beyond the scope of this paper. $`R_t`$ for several other IMFs widely used in the litterature (Salpeter 1955; Miller & Scalo 1979; Kennicutt 1983) as well as for the Padua stellar models corresponding to $`Z`$=0.0004, 0.004, 0.008, 0.02, 0.05 is tabulated in Jungwiert (1998). Here we note only a few facts. For $`Z=0.02`$, the highest gas return over the Hubble time (let’s say 15 Gyr) is obtained for the Scalo’s (1998) and Kennicutt’s (1983) IMFs ($`R_{15\mathrm{Gyr}}=0.43`$ and 0.44, respectively), followed by the Miller & Scalo’s (1979) IMF (0.41), while the Salpeter’s (1955) IMF gives a markedly lower value (0.30). In all the cases, the returned mass fraction is considerable. As for the relative contributions of LMS, IMS and HMS to the total mass-loss, they are comparable, so that no category can be considered as significantly dominant. For example, the Scalo’s (1998) IMF and $`Z=0.02`$ stellar models lead to ratios LMS : IMS : HMS $`1:1.25:1`$ for the contributions to $`R_{15\mathrm{Gyr}}`$. The stellar metallicity does not change the mass-loss of the whole CSP too much despite the fact that individual stars have their lifetimes and final masses influenced by Z quite noticeably (see Figs. 16 and 17b in Jungwiert 1998).

### 2.3 Galactic gas return rate

We define, for the purpose of simulations presented in Sect. 4, gas return rate ($`GRR`$) as the analog of star formation rate ($`SFR`$): $`GRR`$ (units of $`\mathrm{M}_{}\mathrm{yr}^1`$) is the rate at which stellar mass is converted into gaseous mass. In general, $`GRR`$ depends on the whole history of star formation:

$$GRR(t)=_0^tSFR(t^{})\dot{M}_\mathrm{n}(tt^{})dt^{},$$
$`(4)`$
where $`\dot{M}_\mathrm{n}(tt^{})`$ is the (normalized) mass-loss rate, at time $`t`$, of a stellar population born at $`t^{}`$. Note that when $`GRR`$ refers to a specified area, $`SFR`$ in Eq. (4) does not refer to the same area: $`GRR`$ at time $`t`$ is not determined by the star formation history of this area but by the star formation history of stellar matter that is now present in this area but has, in general, various orbital histories.

## 3 N-body code with star formation and time-dependent stellar mass-loss

### 3.1 Gravitational interaction and equations of motion

We use the standard particle-mesh (PM) technique (Hockney 1970; Hohl 1971; for a review, see: Hockney & Eastwood 1981; Selwood 1987) to compute the gravitational potential on a Cartesian grid via Fast Fourier Transform (FFT). Forces at grid points are then evaluated by differencing the potential. Both the assignment of particle masses to grid points and the finding of forces between them are carried out by the bilinear “cloud-in-cell (CIC)” interpolation. The equations of motion are integrated by means of the leap-frog algorithm (Hockney & Eastwood, 1981).

The gravitational softening that reduces interparticle forces at short distances (see Pfenniger & Friedli 1993) is that of a Plummer sphere. The choice of the softening length is done with regard to the criteria laid down by Romeo (1994, 1997) to optimize the stability and relaxation properties of N-body disks.

### 3.2 Particle species

There are four particle species in our code. Apart from relatively massive stellar particles (hereafter referred to as “standard stars”), typically representing millions of stars, and relatively massive gaseous particles (hereafter “clouds”) with characteristic masses of $`10^510^6\mathrm{M}_{}`$, we introduce low-mass stellar and low-mass gaseous particles (hereafter nicknamed “starlets” and “cloudlets”) that mediate the exchange of matter – converted from gaseous to stellar phase and vice versa by star formation and stellar mass-loss – between the standard stars and clouds. The definition of mass boundaries between cloudlets and ordinary clouds and between starlets and ordinary stars, and thereby the mass ranges for cloudlets and starlets, are precised in Section 4.1. (item 13).

The standard clouds are meant to simulate the cold gas component. They are modelled as finite-size particles that undergo physical collisions in which they dissipate a fraction of kinetic energy. The relative velocity of two colliding clouds is reduced by factors $`\beta _\mathrm{r}`$, $`\beta _\mathrm{t}`$; the former multiplies the velocity component along the line joining the two clouds (radial component), the latter perpendicular to it (tangential component) (for pioneering use of this “sticky-particles” technique in simulations of galaxies, see Schwarz 1981; among more recent applications and discussion on values of $`\beta _\mathrm{r}`$, $`\beta _\mathrm{t}`$, see e.g. Jungwiert & Palouš 1996).

Star formation is allowed to proceed in our cold gas component in the way described in Sect. 3.3. It results in the production of the above mentioned starlets. These in turn join the standard stellar component by merging with a standard star provided some is near by (typically 200 pc; to be precised below). On the other hand, standard stars lose continuously mass in the form of the cloudlets – as detailed in Sect. 3.4. – that join the cold gas component by merging with a standard cloud provided some is near by. The condition of particle “nearness” for star-starlet or cloud-cloudlet merging ensures that the mass exchange between the stellar and gaseous components is local.

In usual situations the existence of our starlets and cloudlets will be ephemeral since they will merge with a nearby standard particle immediately after they are born. This is technically equivalent to a direct exchange of a fraction of mass between a standard stellar and a standard gaseous particle. However, there are important and frequent astrophysical situations in which stars and cold gas are segregated spatially. For example, in some areas cold gas can be depleted by intense star formation or swept away by gravity torques. There can also be regions of cold gas without any or with only small amounts of stars, for instance the whole galaxies in early evolutionary stages or, in the case of present-day galaxies, the area beyond the periphery of optical disks.

Our starlets/cloudlets enable to follow such situations. If a starlet/cloudlet is produced in an area (defined below) with no standard star/cloud, it leads an independent and possibly long-term existence until it meets the respective standard particle (meanwhile, the starlet is subject to mass-loss as standard stars and it produces cloudlets). The long-living starlets/cloudlets can also meet another starlet/cloudlet. Such events are treated by starlet-starlet and cloudlet-cloudlet mergings. In the course of time, a starlet/cloudlet grown by successive mergings of starlets/cloudlets might approach the mass typical of standard stars/clouds. If this occurs, the starlet/cloudlet is converted into a standard star/cloud. Conversely, a standard cloud is changed over to a cloudlet if its mass falls, due to the starlet production, below the value typical of standard clouds. In contrast, we do not consider a conversion of standard stars into starlets since the mass of a standard star remains of the same order during a simulation: it cannot decrease below the quantity $`1R_{t_{\mathrm{end}}}`$ that is larger than $`0.55`$ for any simulation length $`t_{\mathrm{end}}`$ shorter than 20 Gyr (see Table 2).

In this respect, the cloudlets can be viewed as a very crude representation of the warm/hot galactic gas: they are products of star formation and do not participate in the star forming process unless they “cool” in areas of high gas density (either by merging with a standard cloud or with many other cloudlets).

In summary, the four particle species have the following properties. Standard stars produce cloudlets and can absorb starlets. Starlets can be absorbed by standard stars or merge with other starlets; they produce cloudlets; if their mass exceeds a given value ($`M_{\mathrm{sS}}`$, $`10^510^6\mathrm{M}_{}`$, see Sect. 4.1., item 13), they are converted into standard stars. Standard clouds undergo mutual, partly inelastic collisions; they produce starlets and can absorb cloudlets; if their mass falls below a given value ($`M_{\mathrm{cC}}`$, $`10^5\mathrm{M}_{}`$, Sect. 4.1., item 13), they are converted into cloudlets. Cloudlets can be absorbed by standard clouds or merge with other cloudlets; if their mass exceeds the above mentioned $`M_{\mathrm{cC}}`$, they are converted into standard clouds.

All the particle productions, absorptions, mergings and collisions happen in the mass and momentum conserving manner.

Perhaps the most questionable interaction in our code is the merging of a starlet with a standard star or another starlet. The real stellar fluid cannot dissipate energy; nevertheless, if new stars are born with a low velocity dispersion (as expected from the kinematics of parent gas clouds), the velocity dispersion of the real stellar fluid as a whole can decrease. We merge the particles for computational reasons, so as to keep their number low (see below); we are unable to follow every new born stellar particle independently. In this context, the merging of standard stars with starlets can be viewed as very crudely mimicking the rejuvenation of the stellar fluid by newly born stars and the related decrese of its velocity dispersion. The merging of one starlet with another can seem more problematic but this process is relatively rare in simulations we present. However, we admit that the merging of stellar particles could in general bias the evolution of the stellar orbital structure, especially if the merged particles had very different velocities. This problem certainly merits further investigation.

As indicated above, the interparticle interactions (collisions, mergings) require the ability to identify nearby particles. The same is true for the star formation process described in the subsequent subsection. To find nearby particles in a CPU-efficient way, we divide a galaxy into a 2D rectangular grid (hereafter “interaction grid”; for its specifications see Sect. 4.1., items 11 and 6); four lists (standard stars list, starlets list, standard clouds list and cloudlets list) are associated with it: within each of them, the particles are continuously sorted in such a way that their indices are stored in the order of the increasing index of the grid boxes. This permits a fast identification of all the particles in a given box at every time step. Particles located in a given box are considered as “nearby” in the above used sense. It is true that two particles from distinct neighbouring boxes may be closer to each other than particles within a same box. However, our concern is basically to have at most one starlet and one cloudlet per box and to achieve this in a computationally simple manner. The fact of ignoring interactions behind the box boundary is dynamically unimportant since it happens on distances below the dynamical resolution of the code.

The star-starlet, starlet-starlet, cloud-cloudlet and cloudlet-cloudlet mergings permit to keep a low number of particles. If starlets/cloudlets were not merged, they would soon proliferate in an inacceptable number since at each time step at which the star formation and mass-loss procedures are applied, roughly $`N_{\mathrm{boxes}}`$ ($`10^410^5`$) starlets/cloudlets are created. The merging rule usually leads to mergings of nearly all starlets/cloudlets with standard stars/clouds within the time step in which these starlets/cloudlets were created (see Sect. 4.3.) except areas or evolutionary stages where or when there is not enough stars or cold gas.

### 3.3 Star formation scheme

We impose a generalization of the Schmidt law for star formation (Schmidt 1959; for a recent review, see Kennicutt 1998), i.e. the surface density of star formation rate, $`\mu _{SFR}`$, is proportional to a power of the gas surface density: $`\mu _{SFR}=A\mu _{\mathrm{gas}}^\alpha `$.

Kennicutt (1998) gives $`\alpha =1.4`$ as the “best-fit” for the observationally determined Schmidt law relating the mean surface density of $`SFR`$ to the mean surface density of cold (molecular + atomic) gas in both galactic disks and circumnuclear regions of disk galaxies. In our model, we keep the freedom of choosing $`\alpha `$ as an input parameter (typically 1 or 2). The Schmidt law is applied with a time step $`\mathrm{\Delta }t_{\mathrm{SF}}`$ (longer than the time step used to integrate the equations of motion; Sect. 4.1., item 12) separately in each of the boxes of the interaction grid.

Inside a given box, star formation proceeds as follows: 1) The standard clouds list is used to compute the cold gas surface density. 2) If this density is non-zero (i.e. there is at least one standard cloud), a starlet is created with the mass computed according to the Schmidt law. 3) The starlet mass is subtracted from the standard clouds in amounts proportional to their masses. 4) The position and velocity of the starlet are found in such a way that the total momentum of the standard clouds and of the starlet is conserved. 5) Standard stars, if any, are found using their list. If there are some, one of them is randomly chosen to merge with the starlet. If there is no standard star, the starlets list is searched through for present starlets, if any, that immigrated from neighbouring boxes. If the search is positive, all the starlets in the box merge into one. Otherwise the new starlet will temporarily lead an independent existence as described in the previous subsection.

### 3.4 Stellar mass-loss scheme

In Sect. 2 we have evaluated the normalized mass-loss rate, $`\dot{M}_\mathrm{n}(t)`$, of a coeval stellar population. Here our goal is to implement the mass-loss into the N-body code.

As described in two previous subsections, new stellar mass forms, in our code, from gas as starlets that can merge with nearby standard stars or other starlets. Our typical stellar particle (either a standard one or a starlet) therefore does not represent a CSP but has a complex star formation history since it consists of mass of starlets born at various times in addition to the mass it was assigned in the beginning of the simulation.

The exact treatment of the mass-loss – for the code purposes discretized in time with a step $`\mathrm{\Delta }t_{\mathrm{ML}}`$ (Sect. 4.1., item 12) – would then require to keep in memory, for each stellar particle, its star formation history, i.e. initial masses of starlets that were incorporated into it and the times at which these starlets were born. The mass lost by stellar particle $`J`$ at time $`t_l`$ due to incorporated starlets would read:

$$\mathrm{\Delta }M_J\left(t_l\right)=_{t_l\mathrm{\Delta }t_{\mathrm{ML}}}^{t_l}\underset{k<l}{}m_J\left(t_k\right)\dot{M}_\mathrm{n}\left(t^{}t_k\right)\mathrm{d}t^{}=$$

$$=\underset{k<l}{}m_J\left(t_k\right)\mathrm{\Delta }M_\mathrm{n}\left(t_lt_k\right),$$
$`\left(5\right)`$
where $`m_J(t_k)`$ is the mass with which a starlet, incorporated into $`J`$ at $`t<t_l`$, was created at time $`t_k`$ (the incorporation could have happened later than at $`t_k`$), and $`\mathrm{\Delta }M_\mathrm{n}(t_lt_k)`$ is the fraction of mass lost over $`t_l\mathrm{\Delta }t_{\mathrm{ML}}`$ by the stellar mass born at $`t_k`$. Computing this sum for every stellar particle ($`N_{\mathrm{stars}+\mathrm{starlets}}10^5`$) in each step ($`N_{\mathrm{steps}}1000`$ if mass-loss is applied every $`\mathrm{\Delta }t_{\mathrm{ML}}10`$ Myr over the Hubble time) and keeping in memory the arrays of $`m_{}`$ and $`t_k`$ would be computationally very expensive.

We therefore introduce a much simpler manner to evaluate the mass-loss of stellar particles getting benefit from the fact that $`\dot{M}_\mathrm{n}(t)`$ was shown to be reasonably well fitted by a simple analytical function with only two parameters (Eq. 2). The trick consist in using the same functional form for all stellar particles disregarding their star formation histories, while individually adjusting the two parameters whenever a new starlet is incorporated.

It is thus necessary to define how to change the relevant parameters when two stellar particles of different ages and masses merge, and to check whether or not such a change plausibly represents the sum of their mass-loss rates. Formulated mathematically, we want to find, for a particle $`J`$ formed at time $`t_J`$ by merging of particles $`R`$ and $`S`$, born at $`t_R`$ and $`t_S`$, parameters $`C_J`$ and $`T_J`$ such that the mass-loss rate of $`J`$ (hereafter called the combined mass-loss rate), expressed as

$$\dot{M}_J\left(t\right)=\frac{C_J}{tt_J+T_J},$$
$`\left(6\right)`$
would have similar time dependence and total gas return as if particles $`R`$ and $`S`$ were followed independently. Eq. (5) is then replaced by the corresponding integral of eq. (6):

$$\mathrm{\Delta }M_J\left(t_l\right)=C_J\mathrm{ln}\left(1\frac{\mathrm{\Delta }t_{\mathrm{ML}}}{t_lt_J+T_J}\right)$$
$`\left(7\right)`$
Since there are two free parameters, $`C_J`$ and $`T_J`$, we can impose two conditions on the combined mass-loss. After trying several possibilities (for a more detailed analysis, see Jungwiert (1998), pages 52-58) we require Eq. (7) to give the correct cumulative mass-loss at two time moments: at the end of the simulation, $`t_{\mathrm{end}}`$, and at an intermediate time, $`t_{\mathrm{im}}`$, between $`t_J`$ and $`t_{\mathrm{end}}`$; the choice of $`t_{\mathrm{im}}=t_J+0.2(t_{\mathrm{end}}t_J)`$ proved to give good results. These two conditions can be expressed as

$$_{t_J}^{t_\mathrm{x}}\dot{M}_J\left(t\right)dt=_{t_J}^{t_\mathrm{x}}\dot{M}_R\left(t\right)dt+_{t_J}^{t_\mathrm{x}}\dot{M}_S\left(t\right)dt,$$
$`\left(8\right)`$
where $`t_\mathrm{x}`$ denotes either $`t_{\mathrm{end}}`$ or $`t_{\mathrm{im}}`$.

A priori, it was not ensured that our combined mass-loss approximation well represented the real sum of particles’ mass-losses. However, the tests presented in Sect. 3.6. verify the adequacy of our approach.

Knowing how to combine the mass-loss rates of individual particles when they merge, it remains to be said what values of $`C_J`$ and $`T_J`$ we assign to particles that do not form by merging, i.e. to starlets at the moment of their birth and to standard stars, if any, present in the beginning of the simulation. Such can be viewed as representing coeval stellar populations. Comparing eqs. (2) and (6), we see that $`T_J`$ is simply $`T_0`$ and $`C_J`$ is $`c_0`$ multiplied by the mass of a particle at the moment of its birth, $`t_{\mathrm{birth}}=t_J`$. Values of $`c_0`$ and $`T_0`$ are chosen from Table 1 (as we have mentioned in Sect. 2.1., the fits presented in Table 1 are very close to each other; we will use that for $`\tau =10`$ Gyr).

For the standard stars present in the beginning of the simulation, there is an additional issue consisting in defining their $`t_{\mathrm{birth}}`$; we postpone it to a special subsection below.

The last point is to explain what happens with masses $`\mathrm{\Delta }M_J`$ lost by individual stellar particles at every time step $`\mathrm{\Delta }t_{\mathrm{ML}}`$. We use a strategy similar to that for star formation (Sect. 3.3.). The procedure runs through the interaction grid box by box: 1) All the standard stars and starlets are found using the respective lists and their $`\mathrm{\Delta }M_J`$ (computed according to Eq. (7)) are summed up. 2) If the sum is non-zero (i.e. there is at least one stellar particle in the box), the mass is attributed to a new cloudlet. 3) It is positioned to the center of mass of contributions $`\mathrm{\Delta }M_J`$ coming from individual stellar particles and its velocity is computed in such a way that the total momentum of stellar particles in the box and of the cloudlet itself is conserved. 4) The standard clouds, if any, are found using the clouds list. If there are some, one of them is randomly chosen to merge with the cloudlet. If there is no standard cloud, the cloudlets list is searched through for cloudlets, if any, that immigrated from neighbouring boxes. If the search is positive, all the cloudlets in the box merge into one. Otherwise the new cloudlet will temporarily lead an independent existence.

### 3.5 Mass-loss of pre-existing stars

The purpose of this subsection is two-fold: first, to drive attention to the fact that when a simulation of a galaxy starts from a state in which a stellar component is already present, as is the case of numerous published works, it is natural and important to take into account its mass-loss; second, to define $`t_{\mathrm{birth}}`$ and $`C_J`$ for standard stellar particles representing this component (hereafter “pre-existing stars” and “pre-existing stellar component”).

To know at what rate the pre-existing stellar component should lose mass during the simulation, an assumption must be made about the star formation history of the model galaxy. Here we restrict our attention to two extreme cases, a constant star formation rate over some epoch prior to the start of the simulation (defined by $`t=0`$) and a starburst at some moment prior to the start of the simulation. To be more specific, we suppose a simulation of length $`t_{\mathrm{end}}=10`$ Gyr and calculate what gas mass, $`\mathrm{\Delta }M_{\mathrm{gas}}`$, is released during the simulation (i.e. between $`t=0`$ and $`t_{\mathrm{end}}`$) by the pre-existing stellar component. Table 3 gives this mass (relative to the mass of the pre-existing stellar population at $`t=0`$) for two constant $`SFR`$ histories (lasting 5 and 10 years) and two starbursts (occuring 5 and 10 years prior to the start of the simulation):

The values of $`\mathrm{\Delta }M_{\mathrm{gas}}`$ range between 5 and 15% of the pre-existing stellar mass. This can represent a huge amount of gaseous matter in absolute terms. For example, starting a simulation with a stellar component of $`410^{10}\mathrm{M}_{}`$ (as will be the case of our illustrative simulation in Sect. 4) implies that the gas mass released from these stars between $`t=0`$ and $`t_{\mathrm{end}}=10`$ Gyr would range between $`210^9\mathrm{M}_{}`$ and $`610^9\mathrm{M}_{}`$ which is typical of gas mass in disks of observed present-day spiral galaxies. The mass-loss of pre-existing stars is thus quite important and should not be neglected in simulations.

Technically, to be able to use Eq. (7), we need to assign to each of the pre-existing star (having mass $`M_0`$ at $`t=0`$) a birth time, $`t_{\mathrm{birth}}=t_J`$, and a birth mass, $`M_{\mathrm{birth}}`$ (the latter will give us $`C_J`$ after multiplying by $`c_0`$). In the starburst case, all such particles have $`t_{\mathrm{birth}}`$ simply equal to the time moment of the starburst; in the case of a constant $`SFR`$ history, we draw $`t_{\mathrm{birth}}`$ for each pre-existing particle randomly within the period of constant star formation epoch. Once we assign $`t_{\mathrm{birth}}`$, $`M_{\mathrm{birth}}`$ is found by integrating Eq. (2) backwards in time (from $`t=0`$ to $`t_{\mathrm{birth}}<0`$).

Pre-existing stars are thus supposed to have already lost mass $`M_{\mathrm{birth}}M_0`$ prior to the start of the simulation, however their mass-loss will continue also during the simulation.

### 3.6 Tests

Several tests of our mass-loss algorithm are described in Jungwiert (1998). Especially, it has been tested that the combined mass-loss (Sect. 3.4.) of a stellar particle built by merging (“composite particle” for the purpose of this subsection) of two or more stellar particles (“merged particles”) well approximates the sum of individual mass-losses for many mass and age ratios of merged particles and for many subsequent mergings. Five examples are shown in Figs. 2-4. In each row, the left pannel shows the mass-loss rate of a composite particle while the right pannel shows its mass and cumulative mass-loss (i.e. the gas mass it has released). Full lines represent quantities computed as if all the merged particles were followed independently (i.e. according to Eq. 5) while the dashed lines correspond to the combined mass-loss approximation (eqs. 6 - 8) used in our code. All the quantities are normalized to the sum of initial masses of merged particles.

Fig. 2 corresponds to a merging of two stellar particles of the same initial mass. The first one is born at $`t=0`$, the second at $`t=5`$ Gyr when the merging takes place. Fig. 3 shows the case of multiple mergings, equidistant in time, with all the merged particles having the same initial mass. The number of particles, $`N_{\mathrm{merge}}`$, that have merged into the composite particle is 10 in the top row and 10 000 in the bottom row. Fig. 4 also corresponds to $`N_{\mathrm{merge}}=10`$ (top row) and 10 000 (bottom row) but now the first particle is ten times, respectively ten thousand times more massive than the others it merges with. In all the examples presented in Figs. 2-4 the combined mass-loss approximation well represents the exact mass-loss. A few additional examples can be found in Jungwiert (1998, pages 52-58).

We have also carried out several runs of the N-body model in which the motion of particles was inhibited while the star formation and mass-loss schemes, producing cloudlets and starlets and combining them with standard clouds and stars, were operational. We have verified that these runs yield temporal evolutions for $`SFR`$ and $`GRR`$ close to numerical solutions of models not using any particles.

## 4 An example: a disk galaxy

As an illustrative application of our code, we present 4 simple 2-D models of a disk galaxy having no gas in the beginning of simulations. The first two models (denoted A and B) are rather artificial but useful reference cases: in model A stellar mass-loss is supressed, the disk remains 100% stellar (classical purely stellar N-body model); in model B mass-loss from initial stellar particles produces gas but star formation is supressed. The standard models (denoted C1 and C2) are those in which both star formation and stellar mass-loss operate: model C1 uses the linear Schmidt law ($`\alpha =1`$), model C2 the quadratic one ($`\alpha =2`$); stellar mass-loss is due both to initial and newly formed stellar particles. Table 4 summarizes basic attributes of the models. The simulation length is 10 Gyr for all of them.

### 4.1 Input data

For the sake of completness and reproducibility, we enumerate, before describing the simulations, the code input parameters and specify their numerical values. The majority of them will be common to a larger set of models that will be published in a separate paper. Items 1-7 refer to parameters characteristic of collisionless N-body models, items 8-13 are specific to our gas dynamics and star formation/mass-loss schemes.

1) Initial stellar surface density (2D Toomre-Kuzmin disk: $`M_{\mathrm{disk}}=410^{10}\mathrm{M}_{}`$, scale-length of 4 kpc, truncation radius of 16 kpc);

2) Initial number of stellar particles, $`N_{\mathrm{stars}}=50000`$; all of them have initially the same mass $`M_0=M_{\mathrm{disk}}/N_{\mathrm{stars}}=810^5\mathrm{M}_{}`$.

3) Initial velocity distribution. We fix the Toomre (1964) stability parameter $`Q_{\mathrm{stars}}=1.5`$ (constant in space), and then obtain velocities of particles by finding the time-independent axisymmetric solution of the Jeans equations; the ratio of azimuthal-to-radial velocity dispersions is computed using the epicyclic approximation (see e.g. Binney & Tremaine 1987);

4) Bulge and halo masses ($`M_{\mathrm{bulge}}=M_{\mathrm{disk}}`$, $`M_{\mathrm{halo}}=510^{11}\mathrm{M}_{}`$) and scale-lengths (Plummer spheres of 1 and 25 kpc, respectively). The bulge is implemented essentially to partly stabilize the inner disk, the halo to keep the rotation curve flat in the outer disk. Both components (introduced by analytical formulae) are treated as unevolving;

5) Gravitational softening length (Sect. 3.1.), $`s0.2`$ kpc. It corresponds to the Romeo’s (1994, 1997) “safe” softening computed for the center of our disk.

6) Number of divisions of the Cartesian grid used to compute the gravitational potential (Sect. 3.1.): $`ND=256`$ (in each dimension). This number relates to the “active” grid (the computation proceeds on a grid doubled in each dimension as required by the FFT technique). The linear size of a grid cell, $`\mathrm{\Delta }_{\mathrm{grid}}`$, is set equal to the softening length. The above choice of $`ND`$ and $`s`$ gives the linear size of the active grid, $`D_{\mathrm{grid}}=50`$ kpc (allowing for maximum disk radius of 25 kpc);

7) Time step for the integration of equations of motion, $`\mathrm{\Delta }t`$. It is constant over the whole disk. Its value varies with time in such a way that it is 1/25 of the minimum orbital period. Initially $`\mathrm{\Delta }t=0.6`$ Myr. It decreases as the central mass concentration increases;

8) Parameters for cloud collisions (Sect. 3.2.): a) $`\beta _\mathrm{t}=1.0`$. This choice assures the conservation of angular momentum in a collision; b) $`\beta _\mathrm{r}=0.6`$. This choice is rather arbitrary. Varying $`\beta _\mathrm{r}`$ falls beyond the scope of this paper. We just note that a positive sign for $`\beta _\mathrm{r}`$ means, in our notation, that the direction of the cloud-to-cloud momentum in a collision is not reversed, unlike in most works implementing this collisional scheme. For more comments on values of $`\beta _\mathrm{r}`$, see Jugwiert & Palouš (1996);

9) Two parameters for the Schmidt star formation law (Sect. 3.3.): we choose $`\alpha =1`$ or $`2`$ and compute the proportionality constant $`A`$ in such a way that the $`SFR`$ value coincides with that of the Kennicutt’s (1998) fit ($`\alpha =1.4`$) at cold gas surface density of $`5\mathrm{M}_{}\mathrm{pc}^2`$ (a value typical of present day-spiral galaxies);

10) Parameters related to the stellar mass-loss rate (Sects. 3.4. and 3.5): to assign $`t_{\mathrm{birth}}`$ to initial stellar particles, we use a simplifying arbitrary assumption that $`SFR`$ was constant over a period $`T_{\mathrm{SF}}=5`$ Gyr before the beginning of the simulation. $`c_0`$ and $`T_0`$ (that determine $`C_J`$ and $`T_J`$ of individual standard stars and starlets) are taken from Table 1 (for $`\tau =t_{\mathrm{end}}=10`$ Gyr);

11) The interaction grid (Sect. 3.2.) used for cloud collision, star formation and stellar mass-loss schemes has the same size and the same number of cells as the grid described in item 6;

12) Time steps $`\mathrm{\Delta }t_{\mathrm{SF}}`$ (Sect. 3.3.) and $`\mathrm{\Delta }t_{\mathrm{ML}}`$ (Sect. 3.4.) for applying star formation and mass-loss schemes. We fix $`\mathrm{\Delta }t_{\mathrm{SF}}=\mathrm{\Delta }t_{\mathrm{ML}}=6`$ Myr (10 times the initial integration step $`\mathrm{\Delta }t`$);

13) Masses $`M_{\mathrm{sS}}`$ and $`M_{\mathrm{cC}}`$ for star-starlet and cloud-cloudlet conversions (Sect. 3.2.). We fix $`M_{\mathrm{sS}}=(1R_{t=10\mathrm{Gyr}})M_03.310^5\mathrm{M}_{}`$, i.e. equal to the minimum mass to which the cumulative mass-loss $`R_t`$ (Table 2) can drive a stellar particle with initial mass $`M_0`$. As for $`M_{\mathrm{cC}}`$, we set $`10^5\mathrm{M}_{}`$; the overwhelming part of cold gas in galaxies seems to dwell in clouds with higher masses than this value.

### 4.2 Results

The disk evolution in the purely stellar model A and in the standard model C1 is shown in Figs. 5 and 6, respectively. Stellar disks (Fig. 5 and upper row of Fig. 6), initially axisymmetric, first develop transient spiral arms and then a long-lasting bar. Meanwhile, in model C1, stars release gas in the form of cloudlets. These gradually merge and in regions of high gas density form standard clouds (see Sect. 3.2.) that kinematically react to the non-axisymmetric potential and accumulate predominantly in a dense nuclear ring (Fig. 6, bottom row) and in an outer ring (Fig. 6, middle row), connected respectively to the inner and outer Lindblad resonances (ILR, OLR). A hint of another, weaker gaseous ring, surrounding the stellar bar and related to the 4/1 resonance, can also be seen.

For models B and C2 (Fig. 7, left and right pannels, respectively), we show only an advanced evolutionary stage ($`t=9`$ Gyr) that can be compared with the last column snapshots of Figs. 5 and 6. In model B, one can note that the bar is shorter and rounder (than in all the other models), the outer ring is less regular (than in C1 and C2 models) while the nuclear ring has collapsed to a smaller size. All these differences are related to a huge gas release and associated bar weakening (for a review, see Combes 2001) that are discussed below. The disk structure and ring shapes in model C2 are rather similar to those of model C1, except that the standard (i.e. relatively massive) clouds in the nuclear ring are considerably depleted by a more intense star formation (resulting from a higher star formation efficiency at high gas densities).

As indicated under item 10 of Sect. 4.1., we have chosen the assumption of constant $`SFR`$ over 5 Gyr preceding the start of the simulation. Given the initial disk stellar mass of $`410^{10}\mathrm{M}_{}`$, Table 3 implies that the gas mass released by the pre-existing stars over the simulation is about $`610^9\mathrm{M}_{}`$. While our choice of the star formation history is arbitrary and favouring a rather huge mass-loss, Table 3 indicates that other choices would lead to a gas restitution of the same order.

We show in Fig. 8, separately for the central disk region (here defined by $`R<2`$ kpc) and for the outer disk ($`R>2`$ kpc), temporal evolutions of $`SFR`$, $`GRR`$, total (i.e. stellar + gaseous mass) and gas mass fraction. The gas restitution by pre-existing stars in models B, C1 and C2 has important effects. First, it opens the way to star formation. Second, the presence of gas affects the radial transport of matter.

As can be seen in Fig. 8, $`SFR`$ is of the order of $`1\mathrm{M}_{}\mathrm{yr}^1`$ in the disk as a whole and, in the outer disk, it is pretty constant over the whole simulation, except early stages. The relative constancy of $`SFR`$ is a result of the time-dependent recycling as was demonstrated by Jungwiert (1998) on a larger set of models, starting from either a two-component disk with $`M_{\mathrm{gas}}/M_{\mathrm{stars}}=0.1`$ or from a purely gaseous disk. As for $`GRR`$, we note that, except early stages, contributions of the pre-existing stellar population and of stars formed during the simulation are comparable (the latter contribution corresponds to the difference between $`GRR`$ curves for models C1 (or C2) and for the non-SF model B).

The transport of matter towards the center is quantified in Table 5 that shows the mean transfer rate of total (stellar + gaseous) mass into the central region ($`R<2`$ kpc) and the relative strength of the bar, $`Q_{\mathrm{tot}}`$ (defined as the ratio of the maximum bar tangential force over azimuthally averaged radial force of the total force field; Combes & Sanders 1981; Buta & Block 2001) at 10 Gyr.

The radial tranport is the least efficient in model A, as expected, since there is no gas. Model B is the most gas rich since mass-loss operates but star formation is inhibited, however it is not the one with the largest inflow. The most efficient are models C1 and C2. This behaviour can be understood with the help of Fig. 9 that shows evolutions of the bar strength $`Q_{\mathrm{tot}}`$. In model A, the bar instability appears after about 1.5 Gyr, the bar strength steeply grows until about 3 Gyr, then saturates at the value of $`0.33`$ and after 5 Gyr enters a phase of only slow secular weakening (at 10 Gyr, $`Q_{\mathrm{tot}}=0.30`$). The bar evolution in models B and C differs from model A in two respects. First, the bar appears earlier (instability starts before $`t=1`$ Gyr and the saturation is reached at $`2`$ Gyr) since the gas production helps to destabilize the disk. Second, after saturation, there is more important weakening of the bar, especially pronounced in the gas rich model B (bar strength decreases from 0.35 to 0.14). This is connected to the known fact that a too important accumulation of matter near the center disturbs the bar orbital structure (Pfenniger & Norman 1990; Hasan & Norman 1990; Friedli & Benz 1993). The high gas fraction in the disk of model B thus first leads to a rapidly growing central concentration of gas (in the nuclear ring) which in turn disturbs the bar thus slowing down the future mass transfer (see Heller & Shlosman 1996).

The recycling models C1 and C2 thus provide, from the point of view of the efficiency in transferring matter inwards, the best combination of gas mass fraction and bar force, both quantities being interdependent.

Table 6 gives two gas mass fractions at the end of simulations: 1) ratio of gas-to-total mass inside the central region (here, “total” means gaseous + stellar but without the bulge contribution); 2) ratio of gas masses inside the central region and in the disk as a whole. In model B, about one quarter of the central region dynamical mass is in the form of gas. Note also that in models B and C1 roughly one half of the gas mass is located in the central region, i.e., more or less in the nuclear ring.

### 4.3 Computational efficiency and statistics of particle species

The presented algorithm for the continuous stellar mass-loss has a low computational cost. In the standard model C1, with 50 000 standard stars and active grid of size 256 x 256, the mass-loss and star formation related subroutines used about 12% and 3%, respectively, of the total CPU, its main part being spent on calculating the gravitational potential and advancing the particles. The simulations were carried out on a relatively slow one-processor PC Pentium II 366 MHz with one time step, $`\mathrm{\Delta }t`$, taking on the average 0.9 sec.

Fig. 10 gives the statistics (masses and numbers) of the four particle species for the model C1. Fig. 11 shows spatial distributions of starlets and cloudlets.

The number of standard stars was constant. In principle, it cannot decrease since these particles have no way to disappear (except escapes that however did not occur) but it could rise due to starlets-star conversion after a starlet reaches the critical mass $`M_{\mathrm{sS}}`$ (Sects. 3.2. and 4.1., item 13); such conversions however did not happen. The mass in standard stars slowly varies due to their mass-loss and absorption of starlets.

The number of standard clouds quickly (on a time-scale of $`500`$ Myr) grows (due to the initially copious stellar mass-loss) from $`0`$ to $`2000`$ and then continues climbing only mildly (to 3300 at $`t=10`$ Gyr). On the other hand, the overall mass in these particles reaches a maximum ($`1.810^9\mathrm{M}_{}`$) at about 4 Gyr and then declines (to $`1.310^9\mathrm{M}_{}`$ at $`t=10`$ Gyr) due to intense star formation and decreasing mass-loss (cf. Fig. 8).

The number of starlets and the total mass in them are very low throughout the simulation ($`200`$ and $`10^6\mathrm{M}_{}`$) and negligible compared to standard stars. The starlets are thus almost virtual particles: despite a big number of them (of the order of 1000; one new starlet per box which hosts at least one standard cloud) is formed every time step $`\mathrm{\Delta }t_{\mathrm{ML}}`$ (Sect. 4.1., item 12), the overwhelming majority is immediately incorporated to nearby standard stars. The only surviving starlets are located in outer parts of the galactic disk (Fig. 11, left) where the number of standard stars is low for the star-starlet mergings to work with the 100% efficiency.

The situation is somewhat different for cloudlets since their number is relatively high, between 10 000 and 13 000. However, their total mass is (except very early evolutionary stages) much lower (it is scaled up by a factor of 10 in Fig. 10 !) than the mass in standard clouds. Initially, the cloudlets happen to gather $`10^8\mathrm{M}_{}`$ due to the mass-loss from initial standard stars, however they efficiently lose it in favour of the standard clouds (initially only by cloudlet-cloud conversions, later on also by cloud-cloudlet mergings). Cloudlets survive only in regions of low total gas density and most of the time harbour only a few $`10^7\mathrm{M}_{}`$. Looking at Fig. 11 (right panel), one sees that the distribution of cloudlets is largely complementary to that of the standard clouds: cloudlets are depleted in high-gas-density regions, i.e. in the nuclear and outer gaseous rings. Note also that the cloudlet distribution extends beyond the outer edge of the standard cloud distribution (cf. Fig. 6).

Overall, our star-formation and mass-loss schemes with four particle species keep the total number of all particles ($`\mathrm{65\hspace{0.17em}000}`$) close to the initial value ($`\mathrm{50\hspace{0.17em}000}`$) and the various types of the particle species interactions do not put a heavy load on the CPU time.

## 5 Concluding remarks

We have presented a gas dynamical N-body model including star formation and time-dependent stellar mass-loss. The implementation of the latter feature is, both physically and technically, the most innovative element of our simulations and, as such, is the major result of this work. The related description should allow other people studying the galactic evolution by means of N-body simulations to transpose our mass-loss scheme into their codes.

Our code is obviously not free of shortcomings, for instance the ignorance of gas cooling and heating, inherent to sticky-particles representation of gas (we however recall – see Sect. 3.2. – that our subdivision of gas particles into standard clouds and cloudlets takes into account, on a very rudimentary level, at least some aspects of the gas cooling and of the gas transition from hot/warm to cold). Questions about the mergings of stellar particles were pointed out in Sect. 3.2. Other potentially important simplification is the neglection of the growth of metallicity on the stellar mass-loss. This will be remedied in the forthcoming research that will include the chemical enrichment.

We have carried out a simple 2D simulation of a disk galaxy starting with no gas. Its purpose was essentially illustrative. It should demonstrate the functionality and potential applicability of the developed code. Our intention was also to show, how the dynamics of classical purely stellar simulations could be altered if stellar mass-loss were taken into account.

It goes without saying that the choice of special initial conditions – on one hand, zero gas mass in the beginning of the simulation and, on the other hand, presumed constant star formation rate prior to the simulation – is not coherent. Nevertheless, such conditions can be approached in some astrophysical situations: during strong interactions between galaxies, gas is stripped from the outer parts of galactic disks (e.g. Moore et al. 1996) while that falling towards the center can be largely consumed by a nuclear starburst (e.g. Mihos & Hernquist 1996). A quietly star forming galaxy can thus be turned into a state without much gas. It can however be replenished by the stellar mass-loss as is the case of our simple model. Moreover, our results as for the rate of mass inflow, enhanced by the time-dependent recycling (relative to pure N-body calculations), seem not to qualitatively depend on our special initial conditions since they point in the same direction as the results we obtained for two-component disk models (Jungwiert 1998).

This work will be followed by a second paper dealing with the application of the presented code to disk galaxies with the emphasis on comparing the instantaneous recycling approximation with our time-dependent recycling, especially as to effects on star formation rate and radial transport of matter.

###### Acknowledgements.

This research was carried out in the framework of the Czech-French project No. 9902 that is co-financed by the Academy of Sciences of the Czech Republic (AV ČR) and by the Centre National de la Recherche Scientifique (CNRS) and in the framework of the Key Project K1048102 of AV ČR. It was also partly supported by the grant No. A3003705 of the Grant Agency of AV ČR. We thank the anonymous referee for careful reading of the manuscript and for his suggestions.


format: text

text:
# Growing random networks under constraints

\[

## Abstract

We study the evolution of a random graph under the constraint that the diameter remain constant as the graph grows. We show that if the graph maintains the form of its link distribution it must be scale-free with exponent between $`2`$ and $`3`$. These uniqueness results may help explain the scale-free nature of graphs, of varying sizes, representing the evolved metabolic pathways in $`43`$ organisms.

\]

In recent years measurements on a wide variety of networks such as the world wide web, the internet backbone , social networks  and metabolic networks have shown that they differ significantly from the classic Erdos-Renyi model of random graphs . While the traditional Erdos-Renyi model has a Poisson link distribution, with most nodes having a characteristic number of links, these networks have scale-free link distributions following a power law $`p(x)x^\gamma `$. To account for these observations, the traditional Erdos-Renyi growth process  has been replaced by newer processes  relying on the intuitively appealing idea of preferential attachment. These processes, have been extensively studied in .

While models of preferential attachment provide an explanation of the prevalence of scale-free networks, the models are largely endogenous and do not take account of global exogenous selection pressures which might shape the form of evolving and growing networks. Such selection pressures would be especially relevant in a biological context. Recent measurements of the topological properties of graphs representing the metabolic networks of $`43`$ organisms have demonstrated their scale-free nature . These metabolic networks are a rare example of different graphs of varying size shaped by similar selection pressures, and allow the testing of explanations for their generic features.

The main selection-based explanation  for the metabolic network topologies relies on the fact that scale-free networks are robust with respect to random malfunction of nodes . Robustness is identified with the diameter of the network, and scale-free networks maintain their diameter when nodes are eliminated at random. However, while scale-free graphs are robust in this sense, it has not been shown that robust graphs must be scale-free. This leaves lingering the question of why metabolic networks are scale-free.

In this paper we consider the evolution of random graphs under the constraint that the diameter remain constant as the graph grows. We show if the graph maintains the form of its link distribution it must be scale-free with exponent between $`2`$ and $`3`$. These uniqueness results may help explain the (apparently universal) scale-free nature of graphs, of varying sizes, representing the evolved metabolic pathways of different organisms. Our assumptions and results are consistent with experimental findings.

We first present a brief introduction to the study of metabolic networks, review the findings of previous investigations, and present the basic definitions necessary for the rest of the paper.

A cell is a complex system composed of numerous organic constituents thickly interwoven in a web of reactions. The processes underlying the life of the cell, which include the generation of mass and energy, and information transfer, are a result of this network of complex interactions.

Much work has been done on understanding the control processes underlying the workings of a cell . However there are many open fundamental questions. While it is of importance to uncover the fundamental design principles underlying the organization of a cell, progress in this direction has been limited because of the immense complexity and the lack of good abstractions which capture certain aspects of the large scale organization.

One possible abstraction of this web of interactions is to represent the gamut of chemical reactions by a graph, where each node represents a chemical constituent of the cell and a directed edge from one chemical constituent A to another constituent B implies that B is a product of a reaction between A and other chemical constituents. Such a graph representation is referred to as a metabolic network.

Large scale sequencing projects have furnished integrated pathway-genome databases from which metabolic networks can be inferred.

Recently, such databases have been used  to analyze the topological properties of the metabolic networks of $`43`$ different organisms including E-coli (bacterium) and Caenorhibditis elegans (eukaryote). They found remarkable similarities in these properties.. In short they found that these networks were uniformly scale free with exponents between $`2`$ and $`3`$.

We present now a brief recap of basic definitions necessary to understand our results.

A Directed Graph $`G(V,E)`$ is a collection of points V connected by edges E such that each edge points from one point to another. We will also use the word node to denote a point in the graph.

The degree of a node is the number of edges attached to it.

The outgoing degree is the number of edges going out and the ingoing degree is the number of edges coming in.

The link distribution $`p\left(k\right)`$ of a graph is the probability that a given node chosen at random has k edges going into it or going out. Note that there are two different link distributions ingoing and outgoing.

The diameter of a graph is the average number of steps in takes to go from a node to any other node.

The nth moment $`M_n`$ of a distribution $`p\left(k\right)`$ is defined as

$$M_n=p\left(k\right)k^n$$
(1)
A metabolic network is a directed graph with nodes representing the various chemical species. There is a directed edge from $`A`$ to $`B`$ if $`A`$ participates in a reaction which leads to $`B`$.

We now proceed to the main result of this paper. We study the evolution of a random graph under the constraint that its diameter is constant as the graph grows.

Let us assume that the outgoing link distribution of the growing graph is $`p_k`$, and it has a variable number of nodes denoted by $`N`$. If these graphs are random, their diameter at a size $`N`$ is approximately given by the formula

$$D=\frac{\mathrm{log}N}{\mathrm{log}E}+1$$
(2)
where the quantity $`E`$ is:

$$E=\frac{M_2}{M_1}$$
(3)
The quantity $`E`$ represents the expected degree following a link. From a node $`A`$ on a graph, the expected degree of a node $`B`$, found by a random edge traversal, will in general be different from the average degree of the graph. This is because high degree nodes have more links. Thus, following a link chosen at random, the probability that the resulting node is a high degree node is higher than if the node was chosen at random .

We now impose the constraint that the diameter is constant as the graph grows. This amounts to demanding that $`D`$ be independent of $`N`$ in equation 2. This would be constant with respect to $`N`$ if the denominator scales as $`\mathrm{log}N`$, i.e.

$$\mathrm{log}E=\alpha \mathrm{log}N$$
(4)
where $`\alpha `$ is a constant. Equation (4) implies

$$E=N^\alpha $$
(5)
Note that $`\alpha 1`$. This is because no node can have degree greater than the total number of nodes, $`N^\alpha N`$ or $`\alpha 1`$.

Now we have because of (3) and (5)

$$\frac{^{k_c}k^2p\left(k\right)}{M_1}=N^\alpha $$
(6)
where $`k_c`$ is the largest degree in the graph.

Since $`N`$ is a variable here, we can differentiate both sides of (6) with respect to $`N`$ in order to say something about the relationship between the various quantities. Note that $`k_c`$ is dependent on $`N`$ in some way. This differentiation gives an equation entirely in terms of $`k_c`$ since the derivative of the integral depends only on the value of the integrand at the endpoints.

Differentiating the left-hand side of (6) gives

$$\frac{}{N}\left(\frac{M_2}{M_1}\right)=\frac{1}{M_1}\left(\frac{M_2}{N}\right)\frac{M_2}{M_1^2}\left(\frac{M_1}{N}\right)$$
(7)
For a normalizable distribution $`M_2M_1^2`$ by definition. In the worst case, when $`M_1`$ is the largest it can be, $`M_2=M_1^2`$. Substituting into (7), the first term is $`2\frac{M_1}{N}`$ and the second term is $`\frac{M_1}{N}`$. Thus, in the worst case, neglecting the second term does not change the scaling of the left-hand side of (7) with respect to $`N`$. We may thus simplify our equation by dropping the derivatives of $`M_1`$ and substituting $`M_1=a`$ where $`a`$ is a constant. (In most real cases, the average degree ($`M_1`$) is only weakly dependant on size, making this realistic.) The resulting equation is:

$$k_c^2p\left(k_c\right)\frac{k_c}{N}=a\alpha N^{\alpha 1}$$
(8)
To derive the nature of $`p_k`$ we need to know something about the dependance of $`k_c`$ on $`N`$. We observe that the expected degree $`E`$ of a node chosen at random will always be smaller than the largest degree $`k_c`$. Also note that the largest degree $`k_c`$ will be smaller than the total number of nodes $`N`$. Thus $`k_c`$ is bounded below and above by power laws

$$EN^\alpha k_cN$$
(9)
which means that $`k_c`$ itself will scale with $`N`$ as a power law according to some intermediate exponent $`\beta `$ as $`k_c=bN^\beta `$ where $`\alpha \beta 1`$

Putting this into (8) we get an equation describing the function $`p()`$ in terms of $`N`$.

$$p(bN^\beta )=\left(\frac{a\alpha }{\beta b^2}\right)\frac{1}{N^{3\beta \alpha }}$$
(10)
Substituting $`bN^\beta =x`$ we get

$$p(x)=\left(\alpha ab^{1+\frac{\alpha }{\beta }}\right)\frac{1}{x^\gamma }$$
(11)
where $`2\gamma =3\frac{\alpha }{\beta }3`$. This shows that under the constraint that a growing graph has a constant diameter the probability distribution assumed constant in functional form $`p(x)`$ must have a power law distribution where the probability of a node having $`x`$ links is inversely proportional to $`x`$ with an exponent $`\gamma `$ between $`2`$ and $`3`$.

We next consider our results in the context of evolved metabolic networks. Jeong, et. al.  have measured link distribution, average degree and diameter for the metabolic networks of $`43`$ different organisms. That they found the link distributions to be uniformly scale-free with exponents between $`2`$ and $`3`$. Furthermore, they found that the diameter was constant with respect to size (see Figure 1).

The fact that the metabolic network diameter is constant across sizes suggests an evolutionary selection pressure on the organism as it evolved. Our results suggest that such a constant diameter constraint (the possible biological reasons for which we discuss later) leads to a scale-free link distribution with exponents between $`2`$ and $`3`$, as has been observed. Thus, our results help to explain why such networks are likely to be scale-free.

We now justify the various formulae and assumptions used in the derivation. Figure 1 shows the metabolic network diameters as a function of the size of the network. As is evident, the diameter is constant across sizes at around $`3.4`$. According to the Newman, et. al. formula (2), the diameter for these graphs should be around $`3.29`$ ***This is because the typical cutoff $`k_c`$ for the metabolic network graphs scales as $`N^{\frac{1.2}{\gamma }}`$ according to the data (see Figure 2). Using this fact(Figure 2) the quantity $`E`$ (see equation (3)) is easily calculated in terms of $`N`$ for the metabolic networks as $`EN^{\frac{3.6}{\gamma }1.2}`$. Substituting this into equation (2) we get
$$D=\frac{1}{\frac{3.6}{\gamma }1.2}+1$$
Substituting $`\gamma 2.2`$, consistent with measurements, yields $`D=3.29`$.. Since these quantities are close it means that the Newman et. al formula (2) applies.

For further validation of our underlying assumptions, Figure 2 demonstrates that for these data, the cutoff $`k_c`$ does scale algebraically with $`N`$ (consistent with our assumption that $`k_cN^\beta `$). In the figure $`k_c`$ has been obtained from the graph data by using the average degree $`d`$ and the exponent of the power law $`\gamma `$, which are all related by the formula $`d=kp_k\frac{1}{\gamma 2}\left(1\frac{1}{k_c^{\gamma 2}}\right)`$.

Figure 3 shows the variation of the average outgoing degree for the metabolic networks on a log-linear plot. Because the data is approximately linear on a log-linear scale, we conclude that the average outgoing degree is proportional to $`\mathrm{log}(N)`$, confirming that the variation of the average degree ($`M_1`$) is weak as we have assumed in the proof.

Figures 1-3 show that the formulae and the assumptions used to derive the scale-free nature hold true for metabolic networks. Our reasoning for the scale-free link distribution observed is consistent with measurements.

Further we provide some biological reasons as to why metabolic networks might be constrained to have small constant diameters. With more steps in the network, pathways are longer. For $`A`$ to be converted to $`B`$ it takes more steps. Because the driving force at each reaction in the cell gets smaller on average, under this scenario (long pathways) the cell would be inefficient at doing its job. It has been shown that certain metabolic paths (ie routes leading from a chemical A to a chemical B) are the shortest possible to carry out a transformation; other pathways can be designed, but involve more steps and intermediates .

According to Fell  minimizing transition times between different states, and reducing the time for perturbations to die out, is a consideration. Watts  has shown that perturbations die out rapidly in networks with small diameters. Thus having small diameters would reduce transition times between different states, suggesting a selective advantage to maintaining the network diameter.

Another reason for maintaining network diameter relates to the minimization of metabolite concentrations . With hundreds of metabolites in cells, their average concentration must be kept low to avoid osmotic and solvation problems. Some metabolites are even toxic if allowed to accumulate. A small diameter leads to quick dissipation.

We have presented a plausible reason for the scale-free distribution observed in metabolic networks, with our assumptions and conclusions being consistent with experiments and with other biological facts. Our argument addresses the issue of why robust networks are likely to be scale-free. Combined with endogenous models of preferential attachment, and the error tolerance of scale-free networks, our results help explain the prevalence of scale-free networks in selective environments.


text:
# PRODUCTION OF LITHIUM, BERYLLIUM, AND BORON BY HYPERNOVAE

## 1 Introduction

An unusual class of very energetic supernovae (“hypernovae,” hereafter HNe) has recently been observed (Iwamoto et al. 1998, 2000). Observationally, these events are identified by their high luminosities and peculiar light curves. Theoretically, these events seem to be the highly energetic core collapse explosion of C-O cores (Iwamoto et al. 1998, 2000; Woosley et al. 1999; Nakamura, Mazzali, Nomoto, & Iwamoto 2001; Tan, Matzner, & McKee 2001). It has been suggested (Iwamoto et al. 1998; Wheeler, Yi, Höflich, & Wang 2000) that these may be associated with at least some gamma-ray bursts (GRBs).

The purpose of this paper is to point out that these very energetic stellar explosions are good sites for the copious production of the light elements lithium, beryllium, and boron (LiBeB). This occurs through the collision of the HN ejecta with the circumstellar medium. Fields et al. (1996) noted that such nucleosynthesis occurs in the explosion of all supernova (hereafter SN) ejecta, when the fastest ejecta collide with the surrounding medium and undergo spallation reactions. Fields et al. found the LiBeB production is particularly large for exploding C-O cores (resulting from, e.g., WR explosions or binary interactions). Even so, for explosions of C-O cores with “normal” energies, the net light element yields are too small to significantly affect the Galactic evolution of LiBeB.

As we will see, for hypernovae the LiBeB production efficiency is much higher, due to: (1) a surface composition of hypernovae which is essentially composed of C and O, ideal parent objects for spallation into lighter isotopes; and (2) a significant fraction of the outer envelope is propelled to high velocities (energies higher than nuclear reaction threshold) due to the very high kinetic energy released in their explosion. For the case of hypernovae, the astrophysical context is reasonably well-defined because there are only two key physical parameters (explosion energy, ejected mass), both of which are constrained by observations of the supernova light curves. Adopting a calculated velocity (energy) spectrum of the ejected C and O it is straightforward to evaluate the absolute yield of light elements by spallation. The only difficulty is that the fast nuclei are slowed down in the course of their propagation, and that the cross sections are energy dependent. The procedure adopted to take into account these effect is explained in Fields et al. (1996).

We combine our theoretical LiBeB yields with a simple model of Be chemical evolution to quantify the hypernova contribution to Galactic Be. By comparing these results with observed Be abundance determinations in very metal poor stars in the halo of our Galaxy, we place an upper limit on the ratio of HNe to Type II SNe. This limit holds assuming a constant HN/SN ratio. In addition, if we assume a correlation between hypernovae and gamma ray bursts, we can constrain the fraction of GRB that can be identified as LiBeB producing hypernovae.

The term “hypernova” has been used in different ways by different authors (e.g., Paczyński 1998; Iwamoto et al. 1998, 2000), so it is important to clarify the meaning used here. We define a hypernova as a core-collapse explosion whose detailed mechanism is unknown, whose kinetic energy is much higher than usual, and whose envelope is dominated by carbon and oxygen (rather than hydrogen and helium). Such events are possibly associated with long-timescale GRBs, and we will discuss this possible association in detail in §4.

## 2 Production of LiBeB by Different Stellar Progenitors

We now compute the production of LiBeB by hypernovae. As seen in Table 1, the bulk properties of HNe are diverse. In particular, the explosion energies and ejected masses apparently span a considerable range. One would expect that the LiBeB yields are very sensitive to both of these parameters. We will show this to be the case, and we will use analytic expressions to derive the scaling of the yields with these parameters.

Fields et al. (1996) noted that the fastest ejecta of a supernova explosion have energies above the thresholds for nuclear spallation reactions. When these fast particles interact with the surrounding medium, they will therefore produce LiBeB. The light element production depends on the velocity spectrum of the explosion, particularly that of the outermost layers. The LiBeB yields also scale as the local ISM (target) density, while the irradiation timescale is that of the ionization energy losses and thus scale inversely with density. These density effects cancel, giving a LiBeB nucleosynthesis which is independent of the local density but which does depend on the fast particle (and ISM) composition. The Fields et al. (1996) study is based on the numerical simulation of the Type Ic event SN 1994I (Nomoto et al. 1994), which is modeled as the explosion of a 2.1 $`M_{}`$ C-O core. In this model the outermost and fasterst layers (the only ones that count in our problem) are a mixture of C and O with no H and a small ($`10\%`$) admixture of He.

The velocity spectrum of the outermost layers can be calculated analytically, as was shown by Imshennik & Nadyozhin (1988, 1989) and reviewed in Nadyozhin (1994); these results have recently been confirmed and extended into the relativistic regime by Matzner & McKee (1999). Fields et al. (1996) give a full derivation of the relevant formulae for our problem; here, we will only summarize key inputs and results. The particle spectrum is determined by the bulk hydrodynamics of the problem, namely the velocity (or kinetic energy) distribution as a function of mass shell $`M(>v)`$. The ejected particles have a energy spectrum $`dN/d\epsilon =m^2v^1dM/dv`$ where $`m`$ is the mean particle mass. Nadyozhin (1994) and Matzner & McKee (1990) find that the velocity spectrum of the fastest ejecta is a power law

$$M(>v)=\zeta ^sM_{\mathrm{ej}}\left(v/v_{}\right)^s$$
(1)
where $`M_{\mathrm{ej}}`$ is the ejected mass and $`v_{}(E/M_{\mathrm{ej}})^{1/2}`$ is a characteristic speed associated with the ejecta of an explosion having energy $`E`$. The constants $`\zeta `$ and $`s`$ take on different values depending on the polytropic index relevant to the problem; for our case of $`n=3`$, we have $`\zeta =1.92`$ and $`s=7.2`$. Eq. (1) holds for the fastest, outermost ejecta, i.e., for $`M(>v)M_{\mathrm{ej}}`$.

The key points here are that (1) the particles follow a steep power law spectrum in kinetic energy per nucleon $`\epsilon `$, with $`dN/d\epsilon \epsilon ^{(s+1)/2}=\epsilon ^{4.1}`$ and (2) the fraction of ejected particles above a particular velocity (or energy) threshold–and thus the fraction available for LiBeB production spallation reactions–scales as the very strong power $`v_{}^s=v_{}^{7.2}`$. Thus, once we adopt the appropriate value of $`s`$, the spectrum of particles is fixed, as are the ratios among the LiBeB isotopes produced for a given projectile and target composition. These results are (almost) independent of the explosion energy or ejected mass, and thus should not vary much from one HN to the next.111In fact, a dependence does remain since the spectrum of eq. (1) is cut off at an energy $`E_{\mathrm{max}}v_{}^2`$. However, this is more difficult to calculate accurately as it depends on the details of shock breakout. Also, for the steeply falling spectra and high energies we consider here, the results are only mildly sensitive to $`E_{\mathrm{max}}`$. By contrast, the total LiBeB yield, e.g., $`m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{HN}}`$, depends very strongly on the explosion energy and ejected mass, with scaling

$$m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{HN}}M_{\mathrm{ej}}^{(s2)/2}E^{s/2}=M_{\mathrm{ej}}^{3.1}E^{3.6}$$
(2)
and thus we can expect strong variations in LiBeB yields among HNe.

So if the hypernova energy is a factor of 10 higher than the usual $`10^{51}`$ erg, then the mass ejected above LiBeB thresholds – and thus the yields – goes up by a factor of $`10^{3.6}`$ = 4000. This suggests that (low-mass) hypernovae can be prolific LiBeB sources. Given the composition and spectrum of the projectiles and the known composition of the target, one calculates the spallation yield in the thick target approximation (excellent in our problem), taking into account the energy dependent spallation cross sections.

A grid of models, comprising various kinds of exploding C-O cores (observed and not) is presented in Table 2. One can verify that the yields obey the scalings given in eq. (2). The most copious LiBeB producers are obviously those with lower mass and higher kinetic energy. Type Ia SNe are a relatively interesting source due to their frequency, but they are less productive than low mass hypernovae, since they eject comparable masses but have 10 times less energy. The high energy explosions of massive hypernovae is overcompensated by their heaviness. Except for the energy, normal Type Ic SNe are events very similar to the SN 1998bw massive hypernova. The extremely large LiBeB production by low mass HNe, if they exist, makes them the most efficient LiBeB-producing events known. As such they could have played a role in the evolution of light elements in the early Galaxy, and possibly the intergalactic medium if there were Pop III HNe. Note that the calculated B/Be ratio (Table 2), around 30, is consistent with the same ratio observed in stars all along the metallicity scale (Duncan et al. 1997, Primas et al. 2000, Cunha et al. 2001). Moreover the isotopic ratios of lithium ($`{}_{}{}^{7}\mathrm{Li}/{}_{}{}^{6}\mathrm{Li}2.1`$) and boron ($`{}_{}{}^{11}\mathrm{B}/{}_{}{}^{10}\mathrm{B}4.1`$) are in good agreement with these observations.

## 3 LiBeB Abundance Constraints on Hypernova Rates

We now turn to the contribution of HNe to the Galactic evolution of LiBeB. From this point of view, it is important to note that HN represent a primary LiBeB production mechanism. That is, due to their self-produced C-O cores, the HNe ejecta are always enriched in C and O, and thus the yields of Be are essentially independent of the ambient interstellar medium metallicity. Consequently, we expect a linear scaling between the HN ejecta of Be and O, $`\mathrm{Be}\mathrm{O}`$, and thus a constant Be/O ratio in the Galaxy. Of course, HNe are not the only primary mechanism; another is LiBeB production via metal-rich particles accelerated in superbubbles (Vangioni-Flam et al. 2000). In addition, standard Galactic cosmic rays, with a composition which reflects the ISM metallicity, give a secondary contribution which does depend on the interstellar metallicity, and so scales as $`\mathrm{Be}\mathrm{O}^2`$.

The relative contribution of primary and secondary processes to Be nucleosynthesis thus depends on the Be-O relation. Unfortunately, O/H is difficult to measure in cool stars, and controversy has arisen as two different O/H (and O/Fe) trends have been claimed. If O/Fe changes in Pop II (e.g., Israelian et al.  igr (1998, 2001); Boesgaard, King, Deliyannis, & Vogt 1999b ; Mishenina, Korotin, Klochkova, & Panchuk mkkp (2000)) Fields et al. (2000) showed that both primary and secondary components are needed, with primary dominating at $`\text{[O/H]}1.5`$, and secondary dominating above. On the other hand, if O/Fe is constant in Pop II (e.g., Carretta, Gratton, & Sneden carretta (2000); Fulbright & Kraft fk (1999)), then a primary source of LiBeB dominates until the roughly solar metallicities (Vangioni-Flam et al. 1998). Thus, regardless of the O/Fe behavior, there is a need for primary Be at some level; the quantitative amount does depend on the details of O data. In what follows we will consider the implications of both possibilities for O/Fe.

One can place LiBeB in full chemical evolution context (e.g., Vangioni-Flam et al. 2000; Fields & Olive 1999) but a simplified approach, appropriate for Pop II, allows one to focus on the physics of the HN contribution to Be. In this approximation, we neglect the (small) astration of Be, and thus the primary production of LiBeB species $`i`$ is described by

$$M_{\mathrm{gas}}\frac{\mathrm{d}}{\mathrm{d}t}X_im_{\mathrm{ej},i}_{\mathrm{HN}}_{\mathrm{HN}}+m_{\mathrm{ej},i}_{\mathrm{SN}}_{\mathrm{SN}}$$
(3)
where $`m_{\mathrm{ej},i}_{\mathrm{HN}}`$ is the mean mass in $`i`$ created by one HN, and $`m_{\mathrm{ej},i}_{\mathrm{SN}}`$ is the same quantity for one (superbubble) SN; the rates of each event are given by $`_{\mathrm{HN}}`$ and $`_{\mathrm{SN}}`$.

Since we are considering primary production, the yields are independent of the initial ISM metallicity, and in fact eq. (3) applies not only to primary LiBeB but also to metals such as O and Fe. Thus we can write

$$\beta \frac{X_{\mathrm{Be}}}{X_\mathrm{O}}=\frac{ϵ_{\mathrm{HN}}m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{HN}}+m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{SN}}}{m_{\mathrm{ej},\mathrm{O}}_{\mathrm{SN}}+ϵ_{\mathrm{HN}}m_{\mathrm{ej},\mathrm{O}}_{\mathrm{HN}}}$$
(4)
where we have assumed a constant ratio

$$ϵ_{\mathrm{HN}}=_{\mathrm{HN}}/_{\mathrm{SN}}$$
(5)
which we will refer to as the “HN rate parameter.”

Given information about spallation and stellar yields, eq. (4) allows us to relate the observed $`X_{\mathrm{Be}}/X_\mathrm{O}16/9\mathrm{Be}/\mathrm{O}`$ to the hypernova rate parameter $`ϵ_{\mathrm{HN}}`$:

$$ϵ_{\mathrm{HN}}=\frac{\beta m_{\mathrm{ej},\mathrm{O}}_{\mathrm{SN}}m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{SN}}}{m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{HN}}\beta m_{\mathrm{ej},\mathrm{O}}_{\mathrm{HN}}}$$
(6)
Unfortunately, eq. (6) as it stands is difficult to evaluate due to the model-dependence of the superbubble $`m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{SN}}`$ and the unknown nature of the HN oxygen yield $`m_{\mathrm{ej},\mathrm{O}}_{\mathrm{HN}}`$. We can still make progress, however, by setting an upper limit to $`ϵ_{\mathrm{HN}}`$, as follows. First, we note that the largest possible oxygen yield is when the ejecta is pure oxygen: $`m_{\mathrm{ej},\mathrm{O}}_{\mathrm{HN}}m_{\mathrm{ej},\mathrm{tot}}_{\mathrm{HN}}`$. We also note that the HN contribution is maximized when we ignore the SN contribution. It thus follows that we may limit the HN rate parameter to be

$`ϵ_{\mathrm{HN}}`$ $``$ $`{\displaystyle \frac{\beta m_{\mathrm{ej},\mathrm{O}}_{\mathrm{SN}}m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{SN}}}{m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{HN}}\beta m_{\mathrm{ej},\mathrm{tot}}_{\mathrm{HN}}}}`$ (7)
$``$ $`{\displaystyle \frac{\beta m_{\mathrm{ej},\mathrm{O}}_{\mathrm{SN}}}{m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{HN}}\beta m_{\mathrm{ej},\mathrm{tot}}_{\mathrm{HN}}}}`$
With eq. (7) in hand, we can now place limits on the HN rate parameter. We adopt the SN oxygen yield $`m_{\mathrm{ej},\mathrm{O}}_{\mathrm{SN}}=2M_{}`$, which is insensitive to the choice of initial mass function. For the total hypernova ejected mass we adopt the large and thus conservative value $`m_{\mathrm{ej},\mathrm{tot}}_{\mathrm{HN}}=10M_{}`$. Finally, we must adopt a Be yield for HNe. As Table 2 illustrates, the wide range of HN masses and energies implies a huge range in Be yields, spanning orders of magnitude. We will adopt $`m_{\mathrm{ej},\mathrm{Be}}_{\mathrm{HN}}10^5M_{}`$, the lower of the two values found for SN 1998bw in Table 2. The energy and ejected mass dependence for this value are as in eq. (2).

With these parameters, we can evaluate eq. (7) once we have made a choice of $`\beta `$. As noted above, this depends on the oxygen data. The weaker limit to $`ϵ_{\mathrm{HN}}`$ comes from the constant O/Fe case, in which Be is primary over all of Pop II. In this case, we have $`\beta 3\times 10^8`$, and thus our fiducial numbers give

$$ϵ_{\mathrm{HN}}6\times 10^3$$
(8)
This evaluation is coherent with the upper limit which can be derived from the beryllium abundance in extremely metal-poor stars, as observed with the VLT by Primas et al. (2000). On the other hand, if O/Fe varies, then the relevant Be/O ratio is that of the primary component, which Fields et al. (2000) showed to be $`\beta 8\times 10^{10}`$. As this is smaller, we get a tighter limit:

$$ϵ_{\mathrm{HN}}1.6\times 10^4$$
(9)
There are various ways one can physically interpret a limit to $`ϵ_{\mathrm{HN}}`$. If one attributes the origin of a HN to a mass effect, then $`ϵ_{\mathrm{HN}}`$ is essentially the fraction (by number) of massive stars which become a HN. If we assume that stars above some lower mass limit $`m>m_{\mathrm{HN}}`$ become a HN, then for a Salpeter mass function (with massive stars in the range $`10M_{}m100M_{}`$) we derive $`m_{\mathrm{HN}}>91M_{}`$ from eq. (8), and $`m_{\mathrm{HN}}>99M_{}`$ from eq. (9). These lie at the upper edge of the allowed range, reflecting the smallness of the HN contribution. The HN origin could also be related to additional physical parameters, such as binary interactions or rotation. In this case, the limit to $`ϵ_{\mathrm{HN}}`$ would reflect not only a mass effect but also the fraction of systems where such conditions are present.

The variation in HN energy and ejected mass will also have an important effect on the limits quoted. Had we adopted weaker explosions, or more massive ejecta, these would lead to lower Be yields and thus weaker limits on $`ϵ_{\mathrm{HN}}`$. One could even imagine turning the problem around: given an independent measurement of $`ϵ_{\mathrm{HN}}`$, one could use these limits to infer the mean $`v_{}`$ for HNe.

## 4 On the Possible Association Between GRBs and Hypernovae

The Burst and Transient Source Experiment (BATSE) on board the Compton Gamma-Ray Observatory (CGRO) detected more than 2500 gamma-ray bursts (hereafter GRBs) from 1991 to 2000. The distribution of these bursts over the sky is highly isotropic but the number of faint bursts is notably smaller than the expected number if the distribution of bursts was homogeneous in a Euclidean universe (Fishman and Meegan 1995 and references therein). These two facts provide strong evidence that GRBs occur at cosmological distance (Paczyński 1991). This cosmological origin is now firmly established for the long bursts ($`>2\mathrm{s}`$) thanks the discovery of their afterglows made possible by the Beppo-SAX satellite. These late and fading counterparts are first detected in the X-ray range, then in the optical and later in the radio range. About twenty optical afterglows have been discovered to date. The redshifts of most of these events have been measured, and range from $`z=0.433`$ (GRB 990712) to $`z4.5`$ (GRB 000131). These values represent either a direct measure of the redshift of the afterglow or in a few cases the redshift of the host galaxy.

The two most popular models for the source of GRBs associate them with the coalescence of two compact objects (NS-NS or NS-BH, Eichler et al. 1989; Paczyński 1991; Narayan et al. 1992; Mochkovitch et al. 1993) or the collapse of a very massive star into a black hole (collapsar, Woosley 1993). Such collapsars could be either the collapse of a single Wolf-Rayet star endowed with rotation, or the merger of the core of a massive star with a black hole or a neutron star, and should lead to hypernovae as defined in this paper. The recent observations of the optical afterglows of long bursts and the observations of their host galaxies provide several pieces of evidence in favor of the association with massive stars: the indication of dust extinction in optical afterglows and gas absorption in X-ray afterglows suggest that GRBs occur near star-forming regions (Paczyński 1998; Bloom et al. 1998). The first direct evidence for the GRB-massive star association comes with the supernova SN 1998bw which is probably associated with GRB 980425 (Galama et al. 1998). GRB 980326 and GRB 970228 might also be associated with supernovae (Bloom et al. 1999; Reichart 1999).

As the sample of long bursts with a determined redshift is still small, the distribution of GRBs as a function of redshift $`z`$ must be estimated indirectly. This has been done by many authors (for a review, see Piran 1999) who fit the observed peak flux distribution assuming a given rate of bursts $`\rho _{\mathrm{GRB}}(z)`$. The other parameters for such a calculation are the luminosity distribution of bursts $`\mathrm{\Phi }(L)`$, which is usually taken to be independent of $`z`$, the assumed spectral shape for the GRBs and the usual cosmological parameters. The burst rate obtained by this method is very uncertain. Whereas the results of these calculations are weakly sensitive to the adopted values of the cosmological parameters (Cohen & Piran 1995), it has been shown that the BATSE sample is not large enough to distinguish between to extreme assumptions : a constant rate

$$\rho _{\mathrm{GRB}}(z)=\rho _0$$
(10)
or a rate proportional to the cosmic star formation rate

$$\rho _{\mathrm{GRB}}(z)\rho _{\mathrm{SFR}}.$$
(11)
This is true in particular when relaxing the assumption that GRBs are standard candles (Krumholz et al. 1998), which was usually made for the first calculations (Sahu et al. 1997; Wijers et al. 1998) but is not supported by the observations.

In this paper we will only consider the case of long GRBs. As we are interested in putting constraints on the association between GRBs and hypernovae, we will assume that the burst rate is indeed proportional to the cosmic star formation rate. We will use the results obtained by Porciani & Madau (2001). They did a recent estimation of the GRB rate under this assumption. They used $`\mathrm{\Omega }_\mathrm{M}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`H_0=65h_{65}\mathrm{km}/\mathrm{s}/\mathrm{Mpc}`$. The spectral shape of the GRBs was given by the so called GRB-function (Band et al. 1993). This is a phenomenological $`4`$-parameter function which is known to fit very well the observed spectra. They used the following parameters : $`\alpha =1.0`$ for the low energy slope, $`\beta =2.25`$ for the high energy slope and $`E_\mathrm{b}=511\mathrm{keV}`$ for the break energy, which corresponds to the typical values obtained by Preece et al. (2000) who did a detailed spectral study of a large sample of long GRBs. We know from the few GRBs with a measured redshift that the luminosity of GRBs is strongly variable but the luminosity distribution is very poorly constrained. Porciani and Madau used a power-law distribution :

$$\mathrm{\Phi }(L)=C\left(\frac{L}{L_0}\right)^\gamma $$
(12)
where $`C`$ is correctly normalized to have $`_0^+\mathrm{}\mathrm{\Phi }(L)𝑑L=1`$. With all these assumptions and using different estimations of the SFR, they found that their best fits give a GRB rate of about 1-2 bursts per million Type II supernovae, i.e.,

$$ϵ_{\mathrm{GRB},4\pi }=12\times 10^6$$
(13)
Clearly, this is much lower than the HN parameter found in the previous section.

However, this does not demand that we reject the possible association between HNe and GRB. If the GRB emission is focussed in an opening angle $`\mathrm{\Omega }`$, the rate parameter in eq. (13) has to be corrected upward by a factor $`(\mathrm{\Omega }/4\pi )^1`$. To reconcile the rates, we require that $`\mathrm{\Omega }/4\pi \mathrm{2\hspace{0.17em}10}^4`$$`10^2`$. This range has a considerable overlap with observed broad distribution of beaming factor which spans $`\mathrm{\Omega }/4\pi 10^310^1`$ (Frail et al. 2001). Moreover, our estimation of the hypernovae rate $`ϵ_{\mathrm{HN}}`$ scales with the HN parameters in the same way as the Be yield, i.e., as $`M_{\mathrm{ej}}(E/M_{\mathrm{ej}})^{s/2}=M_{\mathrm{ej}}(E/M_{\mathrm{ej}})^{3.6}`$. This strong dependence means that a modest change in $`E/M_{\mathrm{ej}}`$ produces a large shift in $`\mathrm{\Omega }/4\pi `$. Thus one might hope to turn the problem around, and use an accurate measure of the mean $`\mathrm{\Omega }/4\pi `$ to infer the mean $`E/M_{\mathrm{ej}}`$.

Despite this encouraging result, important uncertainties remain, especially due to our poor knowledge of the explosion mechanism. If only a subclass of GRB progenitors leads to a hypernova as we have defined, then $`ϵ_{\mathrm{HN}}`$ has to be compared to a fraction only of $`ϵ_{\mathrm{GRB}}`$ and the constrain on the beaming angle or the $`E/M_{\mathrm{ej}}`$ ratio becomes more severe. On the other hand, the opposite situation cannot be excluded: the isotropic envelope expansion that we associate with a HN is always present, but the GRB is produced only when certain unknown conditions allow the acceleration of an ultrarelativistic outflow. In this case $`ϵ_{\mathrm{HN}}`$ has now to be compared with $`ϵ_{\mathrm{GRB}}`$ divided by the fraction of explosions producing a GRB and large beaming angle are allowed, even with the $`E/M_{\mathrm{ej}}`$ ratio that we have adopted here.

## 5 Conclusion

Motivated by recent theoretical and observational interest in hypernovae, we have considered the LiBeB production by these objects. We find that the LiBeB yields are very sensitive to the explosion energy and ejected mass. If these parameters typically take values as found for SN 1998bw, then the Be yields can be very large due to the high explosion energy.

Using the yields found for SN 1998bw, we have calculated the impact of HNe on LiBeB evolution in the galaxy. HNe represent a primary source of Be, and thus are constrainted by the observed primary component of Be vs O. Using the observed Be data at low metallicities, we are thus able to place limits on the HN/SN ratio. If we further associate HNe with GRBs, we can infer a limit on the beaming angle of the GRB emission $`\mathrm{\Omega }/4\pi 10^2`$ which is consistent with independent estimates. This agreement is encouraging, though of course significant uncertainties remain.

Under the simple assumption of a constant HN/SN ratio which has been made to derive these limits, there are potentially important consequences for LiBeB evolution. The Li-Fe relation shows a small slope (Ryan, Norris, & Beers 1999) which is consistent with standard GCR production of Li (Ryan et al. 2000) but within errors also allows room for other primary Li contributions. As the Li-Fe relation is measured more precisely, one may be able to detect or limit the Li production by HNe, in addition to that of standard and superbubble cosmic rays.

The assumption that the HN/SN ratio is constant would be true, if both arise from massive star formation and a constant, universal initial mass function. While this is probably the simplest assumption, other scenarios are possible. For example, if a first generation of HN associated with very massive stars (Population III) has existed, it could have produced proto-galactic beryllium, along with C and O. If these Pop III HNe produce strictly C and O but little or no iron, then this Be component would manifest itself under as a plateau in the Be-Fe correlation at the lowest metallicity (but a linear, primary trend in Be-O). Consequently, if this Be-Fe plateau were observed, it would not necessarily imply that BBN has contributed. In this scenario, the HN rate would not follow the SN rate during the Pop III phase. Thus, this picture could be tested by measuring the cosmic SN and HN rates at high redshifts.

We warmly thank Robert Mochkovitch for illuminating discussions. This work has been supported in part by PICS 1076 from the CNRS.