text:
# Effect of transients in nuclear fission on multiplicity of prescission neutrons

## Abstract

Transients in the fission of highly excited nuclei are studied in the framework of the Langevin equation. Time-dependent fission widths are calculated which show that a steady flow towards the scission point is established, after the initial transients, not only for nuclei which have fission barriers but also for nuclei which have no fission barrier. It is shown from a comparison of the transient time and the fission life time that fission changes from a diffusive to a transient dominated process over a certain transition region as a function of the spin of the fissioning nucleus. Multiplicities of prescission neutrons are calculated in a statistical model with as well as without a single swoop description of fission and they are found to differ in the transition region. We however find that the difference is marginal and hence a single swoop picture of fission though not strictly valid in the transition region can still be used in the statistical model calculations.

1 Introduction

The fission dynamics of hot compound nuclei continues to be a subject of considerable interest essentially due to the fact that the fission life time determined from the phase space argument of Bohr and Wheeler turns out to be too small to allow for the rather large number of experimentally observed light particles evaporated prior to fission . It was therefore felt necessary to look beyond the statistical model and this gave rise to a revival of interest in the original work of Kramers who considered fission as a diffusive probability flow over the fission barrier . Dissipative dynamical models of fission were subsequently developed to include particle emission and were found to reproduce a host of experimental data .

In a dissipative dynamical model, the dynamics associated with the fission degree of freedom is usually considered to be similar to that of a Brownian particle floating in a viscous heat bath. The heat bath in this picture represents the rest of all the other nuclear degrees of freedom which are assumed to be in thermal equilibrium. The interaction between these large number of intrinsic degrees of freedom and the fission degree of freedom gives rise to a random force and consequently a dissipative drag on the dynamics of fission . The fission trajectory can thus be followed by solving the appropriate Langevin equation . Alternatively, one can study the time evolution of the phase space density of the fission coordinates using the Fokker-Planck equation .

Several time scales are of physical significance as a compound nucleus makes its journey from its formation to its scission. Considering an ensemble of compound nuclei which are all formed at a given instant, it requires a certain interval of time in order to develop a steady probability flow at the saddle point across the fission barrier. During this time interval, also referred to as the transient time ($`\tau `$) , the probability flow at the saddle point increases from zero to its statioanry value. This stationary probability flow also defines the stationary fission width ($`\mathrm{\Gamma }_0`$) and the associated fission life time ($`\tau _f=\mathrm{}/\mathrm{\Gamma }_0`$). However, the entire fission process itself becomes a transient when the compound nucleus has no fission barrier. In such cases, it is usually argued that the full distribution reaches the scission point in a single swoop . The compound nucleus then survives exactly the swooping down time ($`\tau _s`$) in this picture. These transients in nuclear fission were considered earlier in great details in a series of publications . Using a number of simplifying assumptions, analytical expressions were obtained for the time dependent fission widths for both the weakly damped and strongly damped motions . The stationary value of the fission width was obtained by Kramers earlier . The transient time for cases where there is no barrier was also obtained .

An highly excited nucleus can emit a number of light particles (including photons) before it undergoes fission. In order to trace the time evolution of such a hot nucleus, it would therefore be most desirable to couple a dynamic theory of fission with the statistical emission of light particles. Extensive work has been done in this direction using either the Fokker-Planck or the Langevin equations and the importance of the dynamic nature of the fission process at high excitations has been established . Such calculations however are quite involved and require a large computer time. An alternative and easier approach would be to perform a statistical model calculation by modifying a cascade code in which fission is treated as one of the decay channels and the time dependence of the fission width is explicitly taken into account . In such calculations, the input fission widths and the transient times are usually taken from the analytical expressions .

In the present work, we would examine certain issues related to the time dependence of fission widths and its effect on the multiplicity of the prescission neutrons. First, we would study the effect of lowering the fission barrier on the time dependence of the rate of fission. We would use the Langevin equation to model the fission dynamics for our purpose. The motivation for this study is to find the transition from a diffusive process in the presence of a fission barrier to a transient dominated picture when there is no fission barrier. We would indeed find that the diffusive nature of fission continues to some extent even for cases which have no fission barrier. The underlying physical picture that would emerge for fission in the absence of a fission barrier would be as follows. Consider an ensemble of fission trajectories which have started together sliding down the potential (with no fission barrier) towards the scission point. However, the random force acting on the trajectories will introduce a dispersion in their arrival time at the scission point. In other words, the trajectories will cross the scission point at different instants and a flow will thus be established at the scission point. However, the effect of this dispersion will be reduced when the conservative force becomes much stronger than the random force. This would happen at very large angular momentum of the fissioning nucleus due to the strong centrifugal force. Therefore, the single swoop picture for fission becomes more appropriate at very large values of spin of the nucleus. We would establish the above scenario in the first part of our work.

Next, we would perform statistical model calculations for prescission neutron multiplicity using the time dependent fission widths as well as the single swoop description of fission. This would be done with the aim of finding how well the statistical model calculations with and without the single swoop assumption agree with each other. We would subsequently calculate the prescission neutron multiplicity in a dynamical model of fission and would compare the results with those obtained from the statistical calculations. Though one would expect the results from the statistical and the dynamical calculations to be the same, there could be some differences due to reasons such as the following one. The neutron width at any instant of time evolution of a hot nucleus depends upon its temperature. However, a part of the total excitation energy is used to build up the collective kinetic energy in the fission channel in the dynamical model. This reduces the available intrinsic excitation energy, and hence the temperature, in the dynamical model compared to the statistical model calculation where the total excitation energy is taken as the intrinsic excitation energy of the nucleus. This could reduce the neutron multiplicity in a dynamical model than that obtained from a statistical calculation. We would ascertain the magnitude of such differences from our calculation.

The time dependent fission widths will be evaluated by solving the Langevin equation numerically in the present paper. We have to resort to numerical methods rather than using analytical solutions due to the fact that analytical solutions assume a constant friction whereas we would use a strongly shape dependent friction force in our calculation. This shape dependent friction is a modified version of the one-body wall friction, the so-called ”chaos-weighted wall formula”, in which the effect of chaos in the single-particle motion is incorporated . The prescission neutron multiplicity and fission probability calculated from Langevin dynamics using the chaos-weighted wall friction were found to agree fairly well with the experimental data for a number of heavy compound nuclei ($`A200`$) over a wide range of excitation energies . Since the chaos-weighted wall friction contains no adjustable parameter to fit the experimental data, we consider this to be the most appropriate dissipative force to describe the dynamics of fission and accordingly we shall use it here. Further, we have already reported a systematic study of fission widths using this friction . This study was confined to cases with fission barriers whereas we would concentrate upon fission in the absence of a barrier in the present work.

We shall give a brief description of the dynamical model for fission in the next section. Section 3 will contain the numerical results. A summary of the results along with the conclusions will be presented in the last section.
2 Langevin dynamics of fission
2.1 Nuclear shape, potential and inertia
We shall choose the shape parameters $`c,h`$ and $`\alpha `$ as suggested by Brack et al.  as the collective coordinates for the fission degree of freedom. However, we will simplify the calculation by considering only symmetric fission ($`\alpha =0`$). We shall further assume in the present work that fission would proceed along the valley of the potential landscape in $`(c,h)`$ coordinates though we shall consider the Langevin equation in elongation $`(c)`$ coordinate alone in order to simplify the computation. Consequently, the one-dimensional potential in the Langevin equation will be defined as $`V(c)=V(c,h)`$ at valley and the coupled Langevin equations in one dimension will be given  as

$`{\displaystyle \frac{dp}{dt}}`$ $`=`$ $`{\displaystyle \frac{p^2}{2}}{\displaystyle \frac{}{c}}\left({\displaystyle \frac{1}{m}}\right){\displaystyle \frac{F}{c}}\eta \dot{c}+R(t),`$ (1)
$`{\displaystyle \frac{dc}{dt}}`$ $`=`$ $`{\displaystyle \frac{p}{m}}.`$ (2)
The shape-dependent collective inertia and the friction coefficient in the above equations are denoted by $`m`$ and $`\eta `$ respectively. The free energy of the system is denoted by $`F`$ while $`R(t)`$ represents the random part of the interaction between the fission degree of freedom and the rest of the nuclear degrees of freedom considered collectively as a thermal bath in the present picture. The collective inertia, $`m`$, will be obtained by assuming an incompressible irrotational flow and making the Werner-Wheeler approximation . The driving force in a thermodynamic system should be derived from its free energy which we will calculate considering the nucleus as a noninteracting Fermi gas . The instantaneous random force $`R(t)`$ is modeled after that of a typical Brownian motion and is assumed to have a stochastic nature with a Gaussian distribution whose average is zero . The strength of the random force will be determined by the friction coefficient through the fluctuation-dissipation theorem. A detailed description of the various input quantities may be found in a recent publication .
2.2 Nuclear friction
In order to define the friction force in the Langevin equation (eq.(1)), we will consider the one-body wall and window dissipation . We will use the chaos-weighted wall friction as the one-body wall dissipation. We shall give here a brief description of the chaos-weighted wall friction, the details of which is given elsewhere . The chaos-weighted friction coefficient $`\eta _{cwwf}`$ is given as

$`\eta _{cwwf}=\mu \eta _{wf},`$ (3)
where $`\eta _{wf}`$ is the friction coefficient as was given by the original wall formula  and $`\mu `$ is a measure of chaos (chaoticity) in the single-particle motion of the nucleons within the nuclear volume and depends on the instantaneous shape of the nucleus . In the present classical picture, this will be given as the average fraction of the nucleon trajectories that are chaotic when the sampling is done uniformly over the nuclear surface. The value of chaoticity is evaluated for a given shape of the nucleus by sampling over a large number of classical nucleon trajectories while each trajectory is identified either as a regular or as a chaotic one by considering the magnitude of its Lyapunov exponent and the nature of its variation with time . It is observed that the value of chaoticity $`\mu `$ changes from 0 to 1 as the nucleus evolves from a spherical shape to a highly deformed one. This implies that the wall friction is very small for near spherical shapes of a nucleus. The physical picture behind this is as follows. A particle moving in a spherical mean field represents a typical integrable system and its dynamics is completely regular. When the boundary of the mean field is set into motion (as in fission), the energy gained by the particle at one instant as a result of a collision with the moving boundary is eventually fed back to the boundary motion in the course of later collisions . An integrable system thus becomes completely nondissipative resulting in a vanishing friction coefficient. On the other hand, the nucleon motion becomes chaotic in a highly deformed nucleus. Consequently, the energy tranfer from boundary motion to particle motion becomes irreversible giving rise to a large friction force acting on the boundary motion. These aspects were established from various numerical investigations  and was incorporated into the wall formula to give the chaos-weighted wall friction .

In the wall and window model of one-body dissipation, the window friction is expected to be effective after a neck is formed in the nuclear system . Further, the radius of the neck connecting the two future fragments should be sufficiently narrow in order to enable a particle that has crossed the window from one side to the other to remain within the other fragment for a sufficiently long time. This is necessary to allow the particle to undergo a sufficient number of collisions within the other side and make the energy transfer irreversible. It therefore appears that the window friction should be very nominal when neck formation just begins. Its strength should increase as the neck becomes narrower reaching its classical value when the neck radius becomes much smaller than the typical radii of the fragments. We however know very little regarding the detailed nature of such a transition. We shall therefore refrain from making any further assumption regarding the onset of window friction. Instead, we shall define a transition point in the elongation coordinate $`c_{win}`$ beyond which the window friction will be switched on. We shall also assume that the compound nucleus evolves into a binary system beyond $`c_{win}`$ and accordingly correction terms for the motions of the centers of mass of the two halves will be applied to the wall formula for $`c>c_{win}`$ . However, it may be noted that while the window friction makes a positive contribution to the total friction, the centre of mass motion correction reduces the wall friction. The two contributions therefore cancel each other to some extent. Consequently, the resulting wall-and-window friction is not very sensitive to the choice of the transition point. We shall choose a value for $`c_{win}`$ at which the nucleus has a binary shape and the neck radius is half of the radius of either of the would-be fragments.

The chaos weighted wall and window friction will thus be given as

$`\eta (c<c_{win})=\eta _{cwwf}(c<c_{win}),`$ (4)
and

$`\eta (cc_{win})=\eta _{cwwf}(cc_{win})+\eta _{win}(cc_{win}).`$ (5)
The detailed expressions for the wall and window frictions can be found in ref. .
3 Results
3.1 Fission widths from Langevin equation
With all the necessary input defined as above, the Langevin equation (eq.(1)) is numerically integrated following the procedure outlined in ref.. Starting with a given total excitation energy ($`E^{}`$) and angular momentum ($`l`$) of the compound nucleus, the energy conservation in the following form,

$`E^{}=E_{int}+V(c)+p^2/2m`$ (6)
gives the intrinsic excitation energy $`E_{int}`$ and the corresponding nuclear temperature $`T=(E_{int}/a)^{1/2}`$ at each integration step. A Langevin trajectory will be considered as having undergone fission if it reaches the scission point ($`c_{sci}`$) in the course of its time evolution. The calculations are repeated for a large number (typically 100 000 or more) of trajectories and the number of fission events is recorded as a function of time. Subsequently the time dependent fission rates can be easily evaluated.

We have chosen the 200Pb compound nucleus for our study which has been experimentally formed at differnt excitation energies in a number of heavy ion induced fusion reactions . Figure 1 shows the calculated time dependent fission widths at different spins of the compound nucleus for a given temperature. A number of interesting observations can be made from this figure. The time dependence of the fission width of the compound nucleus with a spin of 40$`\mathrm{}`$ (and with a fission barrier) is typical of a diffusive flow across the fission barrier which has been studied extensively on an earlier occasion . The fission width is found to remain practically zero till a certain interval of time ($`t_0`$) which essentially corresponds to the interval after which the fission trajectories start arriving at the scission point. The fission width subsequently increases with time till it reaches its stationary value ($`\mathrm{\Gamma }_0`$). We will use the following parametric form for the time dependent fission width in order to enable us to use it in our later calculations,

$`\mathrm{\Gamma }(t)=\mathrm{\Gamma }_0[1exp((tt_0)/\tau )]\mathrm{\Theta }(tt_0)`$ (7)
where $`\tau `$ is a measure of the transient time after which the stationary flow is established and $`\mathrm{\Theta }(t)`$ is the step function. The intervals $`t_0`$ and $`\tau `$ are obtained by fitting the calculated fission widths with the above expression.

We next note in fig. 1 that the nature of the time dependence of the fission width remains almost same even though the fission barrier decreases and subsequently vanishes with increasing spin. At very large values of spin, however, fluctuations appear at the later stages of time evolution. These fluctuations are statistical in nature because the number of nuclei which have not yet undergone fission decreases very fast with increasing time for higher values of spin and therefore introduces large statistical errors in the measured numbers. The magnitude of the fluctuations can thus be reduced by considering a larger number of fission trajectories. In our calculation, we have taken particular care by using larger ensembles at higher values of nuclear spin in order to enable us to check whether a stationary value of the fission width is attained at all.

The above observation is of particular interest since it shows that the diffusive nature of fission persists even for cases which have no fission barrier. This diffusive nature is a consequence of the random force acting on the fission trajectories as we have discussed earlier. As a compound nucleus is formed having no potential pocket in the fission channel, it starts rolling down the potential towards the scission point. However, the random force acting on these fission trajectories introduces a spread in their arrival time at the scission point. The spread in the arrival time of the fission trajectories gives rise to a finite fission width as we find in fig. 1.

In order to further investigate the above diffusive nature of fission, the fraction of the number of compound nuclei which have survived fission is shown as a function of time in fig. 2. We have considered the same compound nuclei as in fig. 1 for this figure. Here we find a gradual shift in the decay rate with increasing spin of the compound nucleus. Specifically, the exponential decay of the number of compound nuclei having a fission barrier (with spins 40 and 56$`\mathrm{}`$) is found to continue for those without fission barriers (with spins 66, 70 and 90$`\mathrm{}`$). Subsequently we have calculated the fraction of the surviving compound nuclei from the Langevin dynamics by switching off the random force. Figure 3 shows this decay in which all the nuclei have the same life time which is simply the swooping down time ($`\tau _s`$) from the initial to the scission configuration. The spread in the life time of the trajectories around this value when the random force is switched on can also be seen in this figure. It may also be noted that for very large values of the compound nuclear spin, the decay is very fast and consequently, the above spread is very small. For such cases, fission is dominated by the transients and can be approximated by a single swoop process.

We shall now investigate the relevance of the different time scales in order to distinguish between the roles of stationary flow and transients in fission. When the fission life time ($`\tau _f=\mathrm{}/\mathrm{\Gamma }_0`$) is much longer than the transient time $`\tau `$, most of the fission events take place after the establishment of a stationary flow. Evidently, this holds for nuclei with a barrier in the fission channel. However, it is also possible to have $`\tau _f>\tau `$ for cases which have no fission barrier. This is illustrated in fig. 4 where the ratio $`\tau /\tau _f`$ is plotted as a function of the spin of the nucleus. Beyond the critical angular momentum $`(l_c)`$ at which the fission barrier vanishes, we find a window of angular momentum where $`\tau _f`$ is indeed greater than $`\tau `$. This window represents the transition region over which the fission dynamics changes from a steady flow to transients. Fission becomes transient dominated for spin values at which $`\tau >\tau _f`$. A single swoop description of fission can be applied for such cases. However, a single swoop picture would be rather inaccurate in the transition region where a steady flow still persists. In the next subsection, we shall explore the consequences of using the single swoop description of fission in statistical model calculations in terms of the multiplicities of prescission neutrons. It would be of interest for our later discussions to locate the transition region with reference to the spin distributions of the compound nuclei formed in heavy ion induced fusion reactions. We have therefore plotted the spin distribution of the compound nucleus 200Pb obtained in the fusion of $`{}_{}{}^{19}F`$+181Ta at two excitation energies. It is observed that the transition region lies beyond the range of the spin distribution when the compound nucleus is excited to 78 MeV, whereas it is well within the range of the spin values populated at an excitation of 132 MeV. One would thus expect that the number of prescisssion neutrons would be affected more at higher excitation energies when the single swoop picture is used in the transition region.
3.2 Prescission neutrons from dynamical and statistical model calculation
We shall now consider the emission of prescission neutrons from the Langevin dynamics of fission as well as from a statistical model calculation where time-dependent fission widths will be used. We shall consider neutron and giant dipole $`\gamma `$ evaporation in the Langevin dynamical calculation following a random sampling procedure . A Langevin trajectory will be considered as undergone fission if it reaches the scission point in course of its time evolution. Alternately it will be counted as an evaporation residue event if the intrinsic excitation energy becomes smaller than either the fission barrier or the binding energy of a neutron. The calculation proceeds until the compound nucleus undergoes fission or ends up as an evaporation residue. The number of emitted neutrons and photons is recorded for each fission event. This calculation is repeated for a large number of Langevin trajectories and the average number of neutrons emitted in the fission events will give the required prescission neutron multiplicity.

The statistical model calculation of prescission neutron emission proceeds in a similar manner where a time-dependent fission width is used to decide whether the compound nucleus undergoes fission in each interval of time evolution. The intrinsic excitation energy at each step is given by the total excitation energy minus the rotational energy since no kinetic energy is associated with the fission degree of freedom in the statistical model and the compound nucleus is assumed to be in its ground state configuration (zero potenial energy). We shall use two prescriptions for the time-dependent fission widths in our calculation. In the first one, we shall use the parametric form of the width given by eq.(6) for all spin values including those for which there is no fission barrier. The parameters $`\mathrm{\Gamma }_0`$, $`t_0`$ and $`\tau `$ are obtained by fitting the numerically calculated time-dependent widths. In the other statistical model calculation, we shall use the above parametric form only for those spin values which have fission barriers. For higher spin values for which there is no fission barrier including those in the transition region, we shall use the swooping down picture. For these cases, we shall numerically obtain the swooping down time $`\tau _s`$ as explained earlier. In this statistical model calculation, neutron and $`\gamma `$ evaporation can take place during this period $`\tau _s`$ while the nucleus will be considered as undergone fission at the end of this interval.

Figure 5 shows the calculated prescission neutron multiplicity at different excitation energies of the compound nucleus 200Pb formed in the 19F + 181Ta reaction. Results shown in this figure are obtained from the dynamical and statistical model calculations which are continued for a period of 300$`\mathrm{}/MeV`$. This time period is not sufficient for all the nuclei in the ensemble either to reach the fission fate or to become evaporation residues. Pushing the Langevin calculation much beyond the above time period becomes prohibitive in terms of computer time. The above time duration is however much longer than the transient times and hence are adequate for our purpose of comparing the dynamical and statistical results.

We first note in fig.5 that the neutron multiplicity calculated from the statistical model using the time-dependent fission widths with and without swooping down assumption are alomost same at lower excitation energies though they differ marginally at higher excitation energies. Such a difference was anticipated in the earlier subsection since the swooping down assumption is invoked more frequently for compound nuclei at high excitation energies which are mostly formed with large values of spin and consequenlty with no fission barrier. In order to explore this point further, the differential neutron multiplicities are obtained from the statistical model calculations with as well as without the single swoop description and are shown in fig. 6. The two calculated distributions at an excitation energy of 132 MeV are found to be different beyond $`l_c`$ though they merge again at the higher end of the transition region. This difference essentially reflects the approximate nature of the single swoop description in the transition region. However, the magnitude of this difference is found to be rather small ($``$ a few $`\%`$). At a lower excitation of 78 MeV, the two distributions are almost identical as one would expect since they have very little overlap with the transition region. The significance of the above observations is of interest since it shows that for compound nuclei without a fission barrier, considering a sharp valued life time (the swooping down time $`\tau _s`$) instead of a life time with a dispersion does not make any appreciable effect in the number of emitted neutrons before fission. It is next observed in fig.5 that the neutron multiplicity from the statistical (both calculations) and dynamical models are also very close to each other though the statistical models marginally overestimate the neutron multiplicity compared to the dynamical model. A possible explanation for this observation would be the fact that the compound nuclear temperature in the statistical model is higher than that in the dynamical model since a part of the total excitation energy is locked up as kinetic energy of the fission mode in the dynamical model. This reduces the intrinsic excitation energy and hence the temperature in the dynamical model resulting in a smaller number of evaporated neutrons.

We have already mentioned that a full dynamical calculation can take an extremely long computer time particularly for those compound nuclei whose fission probability is small. We shall therefore follow a combined dynamical and statistical model, first proposed by Mavlitov et al. , in order to perform a full calculation. In this model, we shall first follow the time evolution of a compound nucleus according to the Langevin equations for a sufficiently long period during which a steady flow across the fission barrier is established. We shall then switch over to a statistical model description after the fission process reaches the stationary regime. It is possible to continue this calculation for a sufficiently long time such that every compound nucleus can be accounted for either as an evaporation residue or having undergone fission.

The prescission neutron multiplicity calculated with the above combined dynamical and statistical model is shown in fig. 7 along with the full statistical model calculations. The statistical model calculations are made with as well as without the swooping down assumption in the time-dependence of the fission widths. The experimental values are also shown in this figure. The observations made in this figure are similar to those in fig.5, namely, the statistical calculations slightly overestimate the neutron multiplicity compared to the dynamical (plus statistical) calculation. However, the statistical and dynamical results are quite close to each other and are also close to the experimental values. This result therefore shows that the statistical calculation with time-dependent fission width can represent the dynamical calculation with reasonable accuracy. We have also shown the results of a statistical calculation in this figure where the fission widths are assumed to be independent of time and are given by their stationary values. This calculation substantially underestimates the neutron multiplicity and illustrates the importance of transients at higher excitation energies.
4 Summary and conclusions
We have presented in the above a numerical study of the transients in the fission of highly excited nuclei and their effect on the number of neutrons emitted prior to fission. To this end, we first investigated the time-dependence of fission widths using the Langevin dynamics of fission. We have shown that the fission width reaches a stationary value after a transient period even for those nuclei which have no fission barrier. We have discussed the role of the random force acting on the fission trajectories in introducing a dispersion in their arrival time at the scission point and thereby giving rise to a finite rate of fission for such cases. We have also shown that this stationary fission rate for very large values of spin of the nucleus loses significance since the stationary fission life time itself becomes much smaller than the transient time for such cases. Therefore, fission of nuclei rotating with a large angular momentum can be considered to proceed in a single swoop. Our study demonstrates a gradual transition from a diffusive to a single swoop picture of fission with increasing spin of the compound nucleus.

We have subsequently examined the effect of the transients on the multiplicity of the prescission neutrons emitted in heavy ion induced fusion-fission reactions. We used both the diffusive description and the swooping down picture separately in statistical model calculations and found close agreement between the two calculated neutron numbers at low excitation energies whereas they differed marginally at higher excitations. It was also shown that the differential neutron multiplicities calculated with and without the single swoop assumption differ only in the transition region though the magnitude of the difference is small. We therefore conclude that the single swoop description of fission can be used in statistical model calculations without making any significant error in the final observables.

We finally compared the number of neutrons calculated from a dynamical model with that obtained from a statistical model in which time-dependent fission widths are used. We found that the statistical model marginally overestimates the neutron numbers than those from the dynamical calculation. We explained this difference in terms of the temperature which is lower in the dynamical model than the statistical calculation. The temperature turns out to be smaller in the dynamical model because the excitation energy is shared between the collective fission mode and the thermal mode in the dynamical calculation in contrast to the statistical calculation where the full excitation energy is assumed to be available in the thermal mode. However, in most of the fission events in the dynamical calculation, the kinetic energy builds up to values which are a little above the fission barrier before it proceeds to fission. Since the values of the fission barrier (typically a few MeV or less) are much smaller than the excitation energies (a few tens of MeV or more) considered here, the temperature differences between the statistical and dynamical calculations remain small for most of the cases. Consequently the difference between the prescission neutron multiplicities calculated from the dynamical and statistical models become small, as we have observed in our calculation.


format: text

text:
# ISO observations of Hickson Compact Group 31 with the central Wolf-Rayet galaxy NGC 1741

## 1 Introduction

Compact groups are small systems of galaxies with each member apparently in close proximity to the other members. They are presumed to be physically associated and often show morphological peculiarities indicative of gravitational interaction. Hickson (1982), by examining Palomar Observatory Sky Survey plates, identified 100 compact groups of galaxies with membership defined by the following criteria: (i) There are more than 4 galaxies in the group within 3 magnitudes of the brightest member. (ii) The group must be isolated so that no external galaxies fall within a radius of 3 times the angular extent of the group. (iii) The total magnitude per arcsec2 of the group is less than 26. Considering their high densities and low velocity distributions, compact groups provide an opportunity to study galaxy interactions and their effects, especially considering their short dynamical lifetimes of $``$ 109 yr (Rubin et al., 1990).

One facet of such interactions is the presence of intense star formation, and in particular the presence of Wolf-Rayet features within the spectra of these objects. Wolf-Rayet galaxies are a subset of emission line and HII galaxies whose integrated spectra have a broad emission feature around 4650Å which has been attributed to Wolf-Rayet stars. The emission feature usually consists of a blend of lines, namely HeI $`\lambda `$4686, CIII/CIV $`\lambda `$4650 and NIII $`\lambda `$4640. The CIV $`\lambda `$5808 line can also be an important signature of Wolf-Rayet activity. Wolf-Rayet galaxies are important in understanding massive star formation and starburst evolution (Schaerer and Vacca, 1998). The Wolf-Rayet phase in massive stars is short lived and hence gives the possibility of studying an approximately coeval sample of starburst galaxies. The first catalog of Wolf-Rayet galaxies was compiled by Conti (1991) and contains 37 galaxies. A large number of additional Wolf-Rayet galaxies have been identified and a recent catalog containing 139 galaxies has been compiled by Schaerer et al. (1999). Wolf-Rayet galaxies are found among a wide variety of morphological types, from blue compact dwarfs (Metcalfe et al., 1996) and irregular galaxies, to massive spirals and luminous merging IRAS galaxies (Rigopoulou et al., 1999), with a range of metallicities and ages (Cervino and Mas-Hesse, 1994; Mas-Hesse and Kunth, 1999).

Hickson Compact Group (HCG) 31 consists of 8 galaxies and was observed using the Infrared Space Observatory111Based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with the participation of ISAS and NASA. (Fig. 1). Hickson (1982) identified galaxies A to D from the Palomar Sky Survey in order of luminosity, with A being the most luminous because it included the central starburst. However follow-up surveys (Rubin et al., 1990; Conti et al., 1996; Iglesias-Paramo and Vilchez, 1997; Johnson et al., 1999, 2000) included the starburst in galaxy C rather than A thus making C the most luminous member, and we follow this convention. Galaxies A and C are also designated together as NGC 1741, one of the most luminous Wolf-Rayet galaxies in the catalog of Conti (1991). Rubin et al. (1990) identified several further members of the HCG 31 group, galaxies E, F, G and a possible member further north, Q. The galaxies C and G were included in the Markarian survey and were designated as Mrk 1089 and Mrk 1090 respectively. HCG 31 is situated at a distance of 54.8 Mpc, assuming a Hubble parameter of 75 km sec-1 Mpc-1. The redshift of galaxy D (Hickson et al., 1992) places it at a distance of 359 Mpc, indicating it is not a bona fide member of the HCG 31 group. We present ISOCAM and ISOPHOT results on HCG 31 as part of an observing programme addressing several galaxies exhibiting Wolf-Rayet signatures. The observations and data reduction are presented in Sect. 2. The results are contained in Sect. 3 and discussed in Sect. 4. The conclusions are summarised in Sect. 5.

## 2 Observations and Data Reduction

The ISO (Kessler et al., 1996) observations were obtained using the mid-infrared camera ISOCAM (Cesarsky et al., 1996) and the spectrometric mode of the ISO photopolarimeter ISOPHOT (Lemke et al., 1996). The astronomical observing template (AOT) used was CAM01 (Blommaert et al., 2001), for the imaging observations and PHT40, for spectrometry. The observing log of the CAM01 observations at 7.7 $`\mu `$m and 14.3 $`\mu `$m using the LW6 and LW3 filters and the PHT40 observations is presented in Table 1.

### 2.1 ISOCAM

Each ISOCAM observation had the following configuration: 3′′ per pixel field of view (PFOV), integration time of 5.04 s per readout, a raster with $`4\times 4`$ spacecraft pointings with a different stepsize on each axis: 24′′ stepsize between raster points and a 6′′ stepsize between raster lines. The ISOCAM detector was used in association with one of a selection of four imaging lenses. These allowed plate-scales of 1.5′′, 3′′, 6′′ and 12′′ per pixel to be chosen, resulting in turn in overall detector fields-of-view of 48′′, 96′′, and 192′′ square. The 12′′ per pixel field-of-view (12′′ PFOV) did not access much more sky area than the 6′′ PFOV due to vignetting of the instrument beam outside a nominal maximum unvignetted field having 3 diameter. The observations discussed here employed the 3′′ per pixel field-of-view, giving an overall detector field of view of 96′′ square. Apart from these details of the detector plate-scale and FOV, the spatial resolution was determined by the point spread function (PSF), which had a full-width in arcseconds at the first Airy minimum of 0.84x$`\lambda `$, with $`\lambda `$ expressed in $`\mu `$m. So for the 14.3 $`\mu `$m filter, the full width of the first Airy minimum was 12′′ resulting in a Rayleigh criterion of 6′′. The corresponding number for the 7.7 $`\mu `$m filter was 3.2′′. All CAM data processing was performed with the CAM Interactive Analysis (CIA) system (Ott et al., 1997; Delaney, 2000), and the following method was applied during data processing: (i) Dark subtraction was performed using a dark model with correction for slow drift of the dark current throughout the mission. (ii) Glitch effects due to cosmic rays were removed following the method of Aussel et al. (1999). (iii) Transient correction for flux attenuation due to the lag in the detector response was performed by the method described by Abergel et al. (1996). (iv) The raster map was flat-fielded using a library flat-field. (v) Pixels affected by glitch residuals and other persistent effects were manually suppressed. (vi) The raster mosaic images were deconvolved with a multi-resolution transform method (Starck et al., 1998). Deconvolution or image restoration is a difficult step in image processing especially in dealing with images with many instrumental and noise defects such as those found in CAM data. Consider a signal $`I`$, the result of the observation of an object $`O(x,y)`$ through an optical system. This system may be characterised by the point spread function (PSF) $`P(x,y)`$. The observed signal results from the degradation, or convolution, of $`P(x,y)`$ with $`O(x,y)`$ and the inclusion of noise $`N(x,y)`$:

$`I(x,y)=(OP)(x,y)+N(x,y)`$

To correct for the degradation of the optical system, $`O(x,y)`$ must be determined knowing $`I(x,y)`$ and $`P(x,y)`$. The presence of the additive noise means that there is an infinite number of solutions to this equation. One of the most widely used techniques for determining a solution is the multi-resolution transform method (Starck et al., 1998). It is based on the idea that a given realization of a random variable carries an amount of information that is inversely proportional to its entropy. Among an infinite number of solutions, the one that maximises its entropy will minimise the amount of information required to make it compatible with the observed data. The deconvolved map yielded resolution enhancements of a factor of 1.2 to 1.6 with the higher values associated with the 7.7 micron map. Spatial resolution for the 7.7 micron map was 2.2′′, and 5′′ for the 14.3 micron map. Similar resolution enhancements by factors of 1.2 to 2.0 were obtained using the multi-resolution transform method on HST and also ISO data (Wu, 1994; Pantin and Starck, 1996). There was a loss of flux in this process for the weaker, point-like sources (D, E). On the other hand, more extended sources such as the galaxies A and C have their fluxes conserved through the deconvolution process. The deconvolved maps are presented here, with the source photometry performed on the non-deconvolved maps. Colour ratios taken on the deconvolved maps are valid because any attenuation of source flux occurs similarly for a given source in each filter. Photometry was performed by integrating the pixels containing source flux exceeding the background by 3$`\sigma `$. Fluxes measured in broadband filters are essentially valid either for the nominally assumed $`\lambda `$-1 source spectrum, or even if the measured and modelled galaxy spectrum is taken into account, since colour correction terms are negligible in this case, compared with the measurement precision. The lowest contour level is approximately twice the standard deviation of the background noise obtained using pixels from around the border of the image. The ISOCAM fluxes presented in Table 2 have a photometric accuracy of about 15%. A comparison of photometry for each source before and after deconvolution is given in Table 3. IRAS fluxes for NGC 1741 are also given for comparison in Table 4 (Moshir et al., 1990). The difference between the 14.3 $`\mu `$m and the IRAS 12 $`\mu `$m fluxes is probably due to the strong 8.5 and 11.3 $`\mu `$m PAH features in the IRAS band. Comparison of images from different filters was obtained by matching the resolution of the two deconvolved images before creating ratio maps. The durations of the observations in Table 1 give the length of the actual observation including instrumental, but not spacecraft, overheads.

### 2.2 PHT-S

PHT-S consists of a dual grating spectrometer with a resolving power of 90 (Laureijs et al., 1998). Bands SS and SL cover the range of 2.5 - 4.8 $`\mu `$m and 5.8 - 11.6 $`\mu `$m respectively. The PHT-S spectrum of the nucleus of NGC 1741 was obtained by operating the $`24^{\prime \prime }\times 24^{\prime \prime }`$ aperture of PHT-S in rectangular chopping mode using the ISOPHOT focal plane chopper. In this mode, the satellite pointed to the position between the source and an off-source position, and the chopper moved alternately between these two positions. The source was always in the positive beam in the spacecraft Y-direction. For the PHT-S spectrum of the lobe of galaxy A, PHT-S was operated in a similar manner. In Figure 1 the nominal PHT-S aperture (24”x24” square) is shown overlaid on the 7.7 $`\mu `$m map, for the case of the off-nucleus measurement. However, the beam profiles differ in detail per PHT-S detector and are in general narrower than the nominal aperture (Laureijs et al., 2001). The calibration of the spectrum was performed by using a spectral response function derived from several calibration stars of different brightness observed in chopper mode (Acosta-Pulido et al., 1999). The relative spectrometric uncertainty of the PHT-S spectrum is about 20% when comparing different parts of the spectrum that are more than a few microns apart. The absolute photometric uncertainty is about 30% for bright calibration sources. All PHT data processing was performed using the ISOPHOT Interactive Analysis (PIA) system, version 8.1 (Gabriel, 1998). Data reduction consisted primarily of the removal of instrumental effects such as radiation events which result in a jump of a minimum of two consecutive read-out values, i.e. in a step within the integration ramp. The disturbance is usually very short and the slope of the ramp after the glitch is similar to the slope before it. Once the instrumental effects had been removed, background subtraction was performed and flux densities were obtained. These fluxes were plotted to obtain the spectra for the central starburst region of NGC 1741 and the nucleus of HCG 31A.

## 3 Results

A map of the HCG 31 group consisting of a deconvolved 7.7 $`\mu `$m LW6 map overlaid on a red Deep Sky Survey (DSS) image of the group is presented in Fig. 1. Deconvolved 7.7 and 14.3 $`\mu `$m maps of the central galaxies overlaid on DSS and H$`\alpha `$ images are presented in Figs. 2 and 3.

The bulk of the 7.7 $`\mu `$m and 14.3 $`\mu `$m emission comes from the interacting galaxies A and C, which are together designated as NGC 1741. Peak emission, which is denoted as N on each map, coincides with the central region of galaxy C. A secondary peak is associated with the nucleus of galaxy A, and forms an extension to the emission from galaxy C. The emission at 14.3 $`\mu `$m is primarily confined to NGC 1741, though galaxy E and a star forming knot within galaxy B are also weakly detected. Since emission at 14.3 $`\mu `$m is a primary tracer of warm dust consisting of very small dust grains stochastically heated in regions of low excitation and/or heated dust close to luminous stars, this indicates that most of the dust in the system is within the galaxy NGC 1741. On the 7.7 $`\mu `$m map, emission is much more widespread, with galaxies B and E well detected along with two sources forming an extension of emission from galaxy A. Since the 7.7 $`\mu `$m PAH feature is indicative of star formation (Charmandaris et al., 1999), this map provides a clear representation of the star formation within the system, a fact outlined by the good agreement between the 7.7 $`\mu `$m contours and the underlying H$`\alpha `$ map (Iglesias-Paramo and Vilchez, 1997). However, a surprising finding was the non-detection of galaxy F on the 7.7 $`\mu `$m map, given its easy detection in H$`\alpha `$. Two other sources were also detected on the 7.7 $`\mu `$m map and are denoted by HGC 31:IR1 and HGC 31:IR2. These sources do not have any known optical or H$`\alpha `$ counterparts, and a careful inspection of the data cube before coadding of all map pointings suggests that they are the residue of cosmic ray glitches that were not completely eliminated, and we are not treating them as confirmed sources. As regards the remaining members of the group, G and Q were outside the ISOCAM map and thus not observed (Fig. 1). The background galaxy D was detected at 7.7 $`\mu `$m, but not in H$`\alpha `$. The Tycho catalogue star TYC 04758-466-1, with mV = 11.3 and spectral type F6 was detected at 7.7$`\mu `$m and 14.3 $`\mu `$m, with a flux ratio that is consistent with the Rayleigh-Jeans tail of a blackbody, which serves as a valuable cross-check on the ISOCAM photometric calibration at these faint fluxes. The ISOCAM fluxes for each source from the non-deconvolved maps are given in Table 2.

A PHT-S spectrum of the central starburst region of NGC 1741 is presented in Fig. 4. Detections were made of unidentified infrared bands (UIBs) at 6.2, 7.7, 8.6 and 11.3 $`\mu `$m that are usually attributed to PAHs. The 8.6 $`\mu `$m feature has a peak flux roughly equal to the 7.7 $`\mu `$m feature, which is unusual in starburst spectra. Such features may be indicative of newly created PAHs (Verstraete et al., 1996). The feature at 10.65 $`\mu `$m may be a blend of features including \[S IV\]. A PHT-S spectrum of the nucleus of galaxy A is presented in Fig. 5. PAHs at 6.2, 7.7, 8.6, and 11.3 $`\mu `$m and a feature at 9.7$`\mu `$m, possibly PAH, were detected but at a level much weaker than in the starburst region. For both spectra, line fluxes for each feature were obtained by interpolating the continuum between the ends of the wavelength range. The line fluxes for the central starburst region are given in Table 5, while the fluxes for the nucleus of galaxy A are given in Table 6.

The spectral energy distribution of NGC 1741 using ISOCAM, PHT-SL from the central starburst and IRAS fluxes is presented in Fig. 6. The dust model, denoted by the solid curve (Krugel and Siebenmorgen, 1994; Siebenmorgen et al., 1998) was fitted to the data and contains three separate dust populations: the PAH bands (Boulanger et al., 1998), a warm dust component at 110 K and a cooler dust component at 40 K, and each component is indicated in Fig. 6.

## 4 Discussion

### 4.1 Morphology and star formation in HCG 31

The interaction between galaxies A and C has been the primary cause of the severely disturbed morphology and high rates of star formation within the group. The most obvious morphological features created by the interaction are the tidal tails; the northern tidal tail extending from A and containing two ISO sources, and a southern tidal tail containing dwarf galaxies E and F. The stripping of H I from galaxies A and C allowed the dwarf galaxies to condense and form within the southern tidal tail. Galaxies B and G have also been affected by the central interaction of A and C. Galaxy B appears quite warped in optical and continuum maps (Rubin et al., 1990; Iglesias-Paramo and Vilchez, 1997), with star formation located in three H$`\alpha `$ knots. The galaxy is also linked to galaxy C by a bridge detectable in H$`\alpha `$. G exhibits star formation weighted towards the edge of the galaxy facing the southern tail, indicating that star formation within the galaxy has been triggered by interaction with the tail. The peculiar system morphology and resultant star formation led to investigations of the star formation history of the group (Iglesias-Paramo and Vilchez, 1997) with three distinct modes of star formation suggested. The strongest mode of star formation, involving galaxies A and C, was a product of the collision between the two, with star formation being concentrated in two strong bursts coinciding with the nuclei of the galaxies. The secondary mode, directly linked to this interaction, involves the formation of dwarf galaxies from H I stripped from both galaxies. Galaxies E and F seem to have been formed in this manner. Finally, a third mode shows star formation displaced along the edges of the affected galaxies, concentrated in knots. This mode, responsible for star formation in galaxies B and G, may be due to a stripping mechanism as the galaxies pass through a denser medium such as H I stripped from galaxies A and C. Given an interaction timescale of about 400 Myr and an age of the central burst of about 4 to 5 Myr (Johnson et al., 1999), it is clear that these modes of star formation are very recent, ongoing, events in the dynamical history of HCG 31. However, the occurance of previous periods of high star formation may be noted by the presence of stellar super-clusters in the central burst with ages over 10 Myr (Johnson et al., 1999, 2000).

NGC 1741 was previously surveyed in X-ray (Stevens and Strickland, 1998), UV (Conti et al., 1996), optical/H$`\alpha `$ (Rubin et al., 1990; Iglesias-Paramo and Vilchez, 1997), radio (Williams et al., 1991) and CO (Yun et al., 1997). It has been estimated that $``$ 700 Wolf Rayet and $``$ 104 O-type stars are located within the central burst (Conti et al., 1996). With some knowledge of the content of the burst, it is instructive to determine the current state of star formation in NGC 1741 using a diagnostic such as the 14.3/6.75 $`\mu `$m flux ratio for the system. The 14.3/6.75 $`\mu `$m ratio generally decreases as interactions develop and starbursts age (Vigroux et al., 1999; Cesarsky and Sauvage, 1999; Charmandaris et al., 1999; Helou, 1999). To determine the 14.3/6.75 $`\mu `$m ratio for NGC 1741, the 6.75 $`\mu `$m (ISOCAM LW2 filter) flux was synthesized from the PHT-S spectrum in the wavelength range of that ISOCAM filter (5 to 8.5 $`\mu `$m). In order to check the accuracy of the synthesised flux, 6.75 $`\mu `$m and 7.7 $`\mu `$m fluxes were synthesised for several sources, including Mkn 297 and NGC 7714 (O’Halloran et al., 2000), from their PHT-S spectra, and the derived values were found to agree with the directly measured ISOCAM fluxes for the sources within the photometric uncertainties. For NGC 1741, the 14.3/6.75 $`\mu `$m ratio was found to be 1.6$`\pm `$0.4. This value is quite low. A ratio above 3 is indicative of high star formation leading to the heating of dust in the nuclear region by hot, young ionising stars (Vigroux et al., 1999; Laureijs et al., 2000).

Since the 14.3 $`\mu `$m flux is dominated primarily by emission from very small dust grains and the 6.75 and 7.7 $`\mu `$m fluxes are dominated by PAHs, the 14.3/7.7 $`\mu `$m ratio provides a diagnostic similar to the 14.3/6.75 $`\mu `$m ratio, allowing a determination of the current state of the starburst. To determine the 14.3/7.7 $`\mu `$m ratio, fluxes were obtained in strips on each of the raster maps with each strip bisecting NGC 1741 through the central starburst region. Each strip was equally separated by a position angle (PA) of 45, measured from north, where PA = 0. The strips were divided to obtain the 14.3/7.7 $`\mu `$m ratio. One such ratio-strip is presented in Fig. 7, taken from northwest to southeast through the galaxy. The slope of the 14.3/7.7 $`\mu `$m ratio changes over the length of the strip and attains a peak value of 1.9 at the center of the starburst. To the northeast of the burst the ratio drops to 1.5, while to the southwest the slope is steeper, reaching a minimum of 1.3. The higher value of the ratio in the central starburst region may be due to destruction of PAHs by UV ionisation from O and B type stars, with an added component due to heating of very small dust grains. In comparison with other galaxies with very active starbursts where the 14.3/6.75 $`\mu `$m ratio can exceed 3, the peak 14.3/7.7 $`\mu `$m ratio value of 1.9 for NGC 1741 is quite low (Charmandaris et al., 1999). The low ratio values and the age of the burst (Johnson et al., 2000) suggest that, while the small central region of NGC 1741 is still in the throes of a very recent starburst, the central region may be moving into a post-starburst stage with cessation of the formation of the most massive stars (Charmandaris et al., 1999). This is in contrast to previous higher resolution surveys which suggest strong ongoing star formation, especially in compact regions (Rubin et al., 1990; Conti et al., 1996; Iglesias-Paramo and Vilchez, 1997; Johnson et al., 2000). However, the lower ratio values for NGC 1741 may be due to ISO being sensitive to a more extended region with less activity. The aging of the burst may be due to supernovae and stellar winds (Taniguichi and Kawara, 1988) disrupting the interstellar medium, preventing further star formation as soon as the first generation of O stars have formed and evolved (Cesarsky and Genzel, 2000).

### 4.2 Does HCG 31 harbour an AGN?

Several surveys of HCGs have led to a debate about the nature of the activity within compact groups. Sulentic and Rabaca (1993) and Venugopal (1995) contest the claim by Hickson et al. (1989) that far infrared emission is enhanced in compact groups, arguing that while tidal interactions and mergers did occur in the past in compact groups (Rubin et al., 1991; Menon, 1995), most compact group galaxies seem to exhibit the behaviour of normal galaxies. However, a significant proportion of the brightest HCGs do display some sort of activity, be it starburst or AGN (Coziol et al, 1998a).

CO surveys of HCGs (Yun et al., 1997; Verdes-Montenegro et al., 1998) have shown that CO deficient groups are also the most dynamically evolved and show a history of strong star formation. An extreme example is HCG 16, which includes one Seyfert 2 galaxy, two luminous LINER galaxies and 3 starburst galaxies (Ribeiro et al., 1996). Those groups with AGN components tend to have their AGN concentrated within the group core. This is likely since AGN tend to be found within luminous and early type galaxies.

Several diagnostic tools are available to probe the nature of the activity within the central region. The ratio of the integrated PAH luminosity and the 40 to 120 $`\mu `$m IR luminosity (Lu et al., 1999) provides a tool to discriminate between starbursts, AGN and normal galaxies because the lower the ratio, the more active the galaxy. For NGC 1741, the ratio is 0.08 and is consistent with the values found for other starbursts (Vigroux et al., 1999). Similarly, the ratio of the 7.7 $`\mu `$m PAH flux to the continuum level at this wavelength can provide a measure of the level of activity within the nucleus (Genzel et al., 1998; Laureijs et al., 2000; Clavel et al., 2000). The ratio for NGC 1741 is 3.5 and is indicative of an active starburst. The results from these diagnostic tools indicate that NGC 1741 is home to a compact burst of star formation and that an AGN is not required within the central region of the merging galaxies HCG 31 A and C.

### 4.3 How will HCG 31 evolve?

Coziol et al. (1998b) discuss a possible galaxy evolution scenario for the HCGs. The first phase involves triggering of starbursts due to mergers, as in HCG 31. When the galaxy’s mass is high enough and a large enough gas reservoir is available, the gas falls inwards towards the centre to form and nourish an AGN. As the starburst fades, the AGN can remain as a Seyfert or LINER depending on the amount of gas available and the level of accretion. When the gas reservoir has finally dissipated, a low luminosity AGN or a normal galaxy may be the final result. Most HCGs are dominated by such low luminosity AGNs or non-active galaxies in their cores indicating that the activity in the system has passed and they are now in a quieter phase.

As regards the eventual fate of HCG 31, models by Barnes (1989) suggest that after continuing mergers, the final result will be an elliptical galaxy. After star formation activity has ceased and the remaining stars have evolved, the luminosity will fade (Rubin et al., 1990). The HCG 31 group sits within a large common H I envelope totaling $``$ 1010 M of gas, which may end up within the center of the object where it can fuel a starburst or an AGN. Given such a large gas reservoir, and the lack of an AGN within the central burst of HCG 31, it would seem that HCG 31 is in the early stages of such a scenario.

## 5 Conclusions

The 14.3/6.75 $`\mu `$m and 14.3/7.7 $`\mu `$m flux ratios (where the 6.75 $`\mu `$m flux was synthesized from PHT-S measurements) suggest that the central burst within NGC 1741 may be moving towards the post-starburst phase. The ratio of the integrated PAH luminosity to the 40 to 120 $`\mu `$m infrared luminosity and the far-infrared colours reveal that despite the high surface brightness of the nucleus, the properties of HCG 31 can be explained in terms of a starburst and do not require the presence an AGN. The star TYC 04758-466-1, with mV = 11.3 and spectral type F6 was detected at 7.7$`\mu `$m and 14.3 $`\mu `$m.

## 6 Acknowledgements

We thank J. Iglesias-Paramo for kindly allowing the use of the H$`\alpha `$ map of HCG 31. The ISOCAM data presented in this paper were analysed using CIA, a joint development by the ESA Astrophysics Division and the ISOCAM Consortium. The ISOCAM Consortium is led by the ISOCAM PI, C. Cesarsky. The ISOPHOT data presented in this paper were reduced using PIA, which is a joint development by the ESA Astrophysics Division and the ISOPHOT consortium with the collaboration of the Infrared Processing and Analysis Center (IPAC). Contributing ISOPHOT Consortium members are DIAS, RAL, AIP, MPIK, and MPIA.


text:
# Future Evolution of Nearby Large-Scale Structure in a Universe Dominated by a Cosmological Constant

## 1 Introduction

Recent data on the temperature anisotropies of the cosmic microwave background (e.g. Hanany et al., 2000; de Bernardis et al., 2000), the luminosity distance to Type Ia supernovae (Perlmutter et al., 1998; Riess et al., 1998; Garnavich et al, 1998), and the large-scale distribution of galaxies (Peacock et al., 2001; Verde et al., 2002) favors a flat universe with present density parameters of $`\mathrm{\Omega }_{M,0}=0.30\pm 0.15`$ in matter and $`\mathrm{\Omega }_{\mathrm{\Lambda },0}=1\mathrm{\Omega }_{M,0}`$ in a cosmological constant. The energy density of the vacuum (the so-called cosmological constant or “dark energy”) appears to be currently dominating the expansion rate of the Universe.

Given these specific values for the cosmological density parameters, it has now become possible to predict quantitatively the future evolution of the visible Universe. The existence of a cosmological constant has profound consequences in this context (see, e.g. Starobinsky, 2000; Gudmundsson & Björnsson, 2002; Loeb, 2002); in particular: (i) when the Universe will age by a factor of a few, the event horizon will stall at a fixed proper distance of $`3.6h_0^1`$ Gpc around us (where $`h_0`$ is the present Hubble constant in units of $`100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$); (ii) all sources with present redshifts larger than 1.8 have already crossed our event horizon and are therefore out of causal contact today; and (iii) even if we continue to monitor the sources with current redshifts of $`z=510`$ into the infinite future, we will only be able to see these sources acquire intrinsic ages of $`46`$ billion years in their rest frame.

In this paper we use N-body simulations to calculate the future gravitational growth of the observed large-scale structure in the local universe. For simplicity, we assume that the vacuum energy density is constant in time (as for a classical cosmological constant) and adopt the values $`\mathrm{\Omega }_{M,0}=0.3`$, $`\mathrm{\Omega }_{\mathrm{\Lambda },0}=0.7`$ and $`h_0=0.7`$. We simulate the evolution of density inhomogeneities within a sphere of radius $`100h_0^1\mathrm{Mpc}`$ around the Milky Way galaxy. Our goal is to determine the time when the accelerated expansion of the Universe will freeze the large-scale structures in distant future, as well as to find the minimum overdensity of matter interior to a shell surrounding an object today that will allow the shell to remain bound to the object at late times despite the repulsive gravitational force of the vacuum. Another issue of particular interest is the question whether the Local Group of galaxies is bound to the nearby Virgo cluster. It has been well-established that the Local Group has a peculiar infall velocity towards the direction of the Virgo cluster, the so-called ‘Virgocentric Infall’ (e.g. Aaronson, et al., 1982; Davis & Peebles, 1983; Lynden-Bell, et al., 1988), but it has not been investigated yet by a direct numerical simulation whether the Local Group will eventually fall into the Virgo cluster in the future given the observed galaxy distribution in the local universe. The possibility remains that the Virgocentric Infall will freeze in comoving coordinates due to the exponential expansion of the universe, and the Local Group will not ultimately merge with the Virgo cluster. We investigate this issue for the first time by running an N-body simulation, starting from initial conditions at z=0 that match the observed galaxy distribution in the nearby universe and ending it in a distant future.

The outline of this paper is as follows. In § 2 we describe the characteristics of our N-body simulations for the evolution of the global structure of the local universe. The results from the simulations are discussed in § 3. We then derive in § 4 the present-day overdensity threshold for a spherical region to collapse in the future, and apply this condition to the Local Group in § 5. Finally, we summarize our main conclusions in § 6.

## 2 Simulations

The simulations were carried out with the parallel tree N-body/SPH code GADGET111http://www.MPA-Garching.MPG.DE/gadget/ (Springel, Yoshida, & White, 2001). The SPH component of the code was turned off as we are only interested in the dark matter particles which dominate the fluctuating mass density. For the initial conditions, we used the simulation output of Mathis et al. (2002)222http://www.MPA-Garching.MPG.DE/NumCos/CR/ at $`z=0`$. This simulation was designed so that at $`z=0`$ the simulated mass density field would best match the observed galaxy overdensity distribution in the IRAS 1.2 Jy survey of Fisher et al. (1994, 1995). Mathis et al. (2002) constructed mock catalogues of the Mark III Catalogue of Peculiar Velocities (Willick et al., 1995, 1996, 1997) using their simulation result, and have shown that the radial peculiar velocities in mock and real catalogues agree very well. The original data set of Mathis et al. (2002) contains 50.7 million dark matter particles of mass $`3.6\times 10^9h_0^1M_{}`$ within the comoving radius of $`80h_0^1\mathrm{Mpc}`$ around the supergalactic center (inner high-resolution region), and 20.5 million particles within the comoving radius of $`200h_0^1\mathrm{Mpc}`$ (outer low-resolution region). The simulation of this spherical region was run with a vacuum boundary condition.

Because we are only interested in the evolution of the structure with mass-scale larger than that of the Local Group, we reduce the original data set by a factor of 100 through random sampling of particles. This reduction of data allows us to run the simulation with a relatively low computational cost while preserving the mass and spatial resolution we need to make predictions for the questions we are interested in, given the observed galaxy distribution in the nearby universe. The resulting mass of each dark matter particle in the high resolution region is $`3.6\times 10^{11}h_0^1M_{}`$ (comparable to the mass of an $`L_{}`$ galaxy), and the number of high resolution particles is half a million. We have set the gravitational softening length to the physical scale of $`100h_0^1\mathrm{kpc}`$, since structures below this scale are not of interest. This mass and spatial resolution is sufficient for following the overall structure of the local universe, but not sufficient for resolving the inner details of galaxy groups such as the Local Group, of which the Milky Way galaxy is the dominant member.

Using the reduced data set, we have simulated the evolution of the local universe from the present time (corresponding to a scale factor value of $`a=1`$ at $`z=0`$) to the 6 Hubble times into the future ($`a=166`$; $`84`$ billion years from today). Because the present Hubble time $`t_H1/H_0=14`$ Gyr (for $`h_0=0.7`$) is very close to the present age of the Universe $`t_0=13.5`$ Gyr in the adopted cosmology, the epoch of $`a=166`$ roughly corresponds to $`t=t_0+6t_H7t_H`$ from the Big Bang ($`t=0`$). Note that we denote all present-day values with a subscript zero.

In addition to the above constrained realization simulation of the local universe, we have also performed a simulation with a periodic box of $`\mathrm{L}=100h_0^1\mathrm{Mpc}`$ containing $`\mathrm{N}=64^3`$ dark matter particles, in order to make sure that the statistical results presented in this paper are not affected by cosmic variance or by the choice of the vacuum boundary condition for the constrained realization run. Statistically, we find that the results from the two runs agree very well with one another, and so we only show the results from the constrained realization run in § 3.

The simulations were performed on the local Beowulf PC cluster located at Harvard-Smithsonian Center for Astrophysics. It is interesting to note that the numerical code speeds up its calculation at late cosmic times. The accelerated expansion of the Universe freezes the large-scale structures in comoving coordinates, and so at late times the variable time step of the code increases and the tree update for the force computation becomes less frequent.

## 3 Future Evolution of Large-Scale Structure in the Local Universe

Figure 1 shows the distribution of particles in a slab covering a width (centered at zero) along the supergalactic Z axis of $`15<\mathrm{SGZ}<15h_0^1\mathrm{Mpc}`$, and projected onto the supergalactic XY plane in comoving coordinates. From top left to bottom right the four panels show snapshots at times $`t=t_0`$, $`t_0+t_H`$, $`t_0+2t_H`$, and $`t_0+6t_H`$, corresponding to $`a=1.0`$, 2.5, 5.8, and 166, respectively. The solid circle in each panel denotes the physical radius of $`100h_0^1\mathrm{Mpc}`$ at each epoch. At $`t=t_0+6t_H`$, the Universe has expanded so much that this circle is no longer visible. Three major density peaks can be easily identified in all panels: the Coma cluster at the top, the Great Attractor (Centaurus) to the left of the center, and the Perseus-Pisces to the right of the center. The Virgo cluster is also apparent just above the center. Close examination of the structures reveals that most of the evolution takes place between $`t_0`$ and $`t_0+2t_H`$, and little evolution is seen after $`t=t_0+2t_H`$. This follows from the fact that after a few Hubble times, the exponential growth of the scale factor, $`a(t)\mathrm{exp}(\sqrt{\mathrm{\Omega }_{\mathrm{\Lambda },0}}t/t_H)`$, damps any peculiar velocity which is gravitationally induced by mild density inhomogeneities. The bottom right panel shows our event horizon by the thick dashed circle around the supergalactic center, which is located at a physical radius of $`3.6h_0^1`$ Gpc. At $`t=t_0+6t_H`$, all clusters except Virgo are already outside our event horizon. The Virgo cluster will also exit from the event horizon at $`tt_0+6.3t_H`$ ($`a=210`$; $`88`$ billion years from today) and its image will fade quickly while remaining frozen on the sky at the time of its exit. Although our calculation does not take general relativistic effects into account, we expect our near-horizon results to be valid, since structures on this scale should at any event be frozen.

The lack of late evolution is even more apparent in the mass distribution of bound objects. Figure 2 shows the cumulative number of dark matter halos with mass above $`M`$ per comoving volume, N($`>`$$`M`$). The dark matter halos are identified through the HOP grouping algorithm (Eisenstein & Hut, 1998) with the outer overdensity threshold of $`\delta _{\mathrm{out}}=80`$ and the peak overdensity threshold of $`\delta _{\mathrm{peak}}=240`$ in comoving coordinates. From left to right, the lines refer to times $`t=t_0`$ (solid), $`t_0+t_H`$ (short-dashed), $`t_0+2t_H`$ (dotted), $`t_0+3t_H`$ (long-dashed), $`t_0+4t_H`$ (dot-short dashed), $`t_0+5t_H`$ (dot-long dashed), and $`t_0+6t_H`$ (short dash - long dash), which correspond to scale factor values of $`a=1.0`$, 2.5, 5.8, 13.5, 31.2, 72.0, and 166, respectively. The last four lines are almost identical within sampling errors, showing that the mass function freezes after two Hubble times from today and that mergers are very rare subsequently.

## 4 Gravitationally Bound Objects

Next we consider the interior overdensity threshold above which a shell of matter around a spherically symmetric object will remain bound to it at late times. In the absence of a cosmological constant and with no shell crossing, the threshold simply amounts to the kinetic$`+`$potential energy of the shell being negative (Gunn & Gott, 1972). The existence of a cosmological constant makes the interior gravitating mass of the vacuum grow larger as the shell radius increases; consequently a shell with a negative energy today may still be pulled away from the object by the cosmic acceleration at late times. We would like to find the minimum mean interior overdensity today, $`\delta _c`$, above which the shell will remain bound to the object at arbitrarily late times. Here $`\delta =(\rho \overline{\rho })/\overline{\rho }`$, where $`\overline{\rho }=\mathrm{\Omega }_{M,0}(3H_0^2/8\pi G)`$ is the mean density of matter in the Universe.

The critical mean overdensity interior to the shell radius, $`\delta _c`$, was derived by Lokas & Hoffman (2002) in their equation (27). Taking the initial time as $`z=0`$, this equation yields $`\delta _c=17.6`$ for the set of cosmological parameters adopted in this paper. We note that Lokas & Hoffman (2002) have ignored the possibility that the mass shell may have a nonzero initial peculiar velocity. For the growing mode of a density perturbation, the inward peculiar velocity will induce collapse at a lower initial overdensity than that found by Lokas & Hoffman (2002).

We have attempted to test the validity of the critical overdensity threshold $`\delta _c=17.6`$ with our N-body simulation. Figure 3 shows the physical radial velocity of particles in and around the two most massive objects in the simulation: Perseus and Centaurus clusters. For the HOP grouping parameter of $`\delta _{\mathrm{out}}=80`$ (see § 3), the total grouped masses of Perseus and Centaurus are $`1.6\times 10^{15}h_0^1M_{}`$ and $`9.0\times 10^{14}h_0^1M_{}`$ at $`t=t_0`$, and $`5.5\times 10^{15}h_0^1M_{}`$ and $`3.0\times 10^{15}h_0^1M_{}`$ at $`t=t_0+6t_H`$, respectively (see § 4 of Mathis et al. (2002) for a detailed comparison between the simulated and observed cluster properties). On the left column of Figure 3 we show the mean physical radial velocity of concentric shells around the density maximum as a function of their physical radius at $`t=t_0+6t_H`$. The short-dashed line indicates the Hubble law at $`t=t_0+6t_H`$, and the horizontal dotted line indicates zero velocity. Inside the radius of $`2h_0^1\mathrm{Mpc}`$ from the cluster center the velocity fluctuates around zero, implying that the particles have virialized. At larger radii, the velocity field approaches the Hubble flow (denoted by the dashed line). Because the physical radius of $`10h_0^1\mathrm{Mpc}`$ at $`t=t_0+6t_H`$ corresponds to comoving $`60h_0^1\mathrm{kpc}`$, dark matter particles at that radius are still infalling towards the density peak and so the physical radial velocity does not perfectly match the Hubble flow.

On the right column of Figure 3 we have selected out all particles that are within a physical radius of $`2h_0^1\mathrm{Mpc}`$ around the density maximum at $`t=t_0+6t_H`$; these particles can be regarded as bound and virialized based on the left column plots. We then trace back the overdensity that these particles had at $`t=t_0`$, and plot it against the radial particle velocity at $`t=t_0+6t_H`$. Zero velocity is indicated by the vertical long-dashed line, and the critical overdensity $`\delta _c=17.6`$ is indicated by the solid horizontal line. For a small number of particles, the value of $`\delta _c=17.6`$ is an overestimate, probably due to them having inward initial peculiar velocities. On average, the analytic estimate for the critical overdensity based on the spherical tophat collapse model appears to provide a good approximation to the actual threshold.

## 5 Future Evolution of the Local Group

### 5.1 Will the Local Group Fall into Virgo?

Although our simulation does not resolve the inner structure of the Local Group, we can get an idea about the future trajectory that its center-of-mass will follow by tracing the particles which reside near the supergalactic center at the present time. In Figure 4, we show the projection of a region with $`15<\mathrm{SGZ}<15h_0^1\mathrm{Mpc}`$ onto the supergalactic XY plane in a similar fashion to Figure 1. The density peak near the top of the panel is the Virgo cluster, which is separated by a comoving distance of $`14h_0^1\mathrm{Mpc}`$ from the supergalactic center. The observational estimates of the distance to the Virgo cluster from us by the methods of surface-brightness fluctuations, planetary nebula luminosity function, Tully-Fisher relation, and Cepheid variable stars, all agree with $`16\pm 2\mathrm{Mpc}`$ (e.g. Jacoby et al., 1992; Graham et al., 1999), in reasonable agreement with the simulated distance within $`2\sigma `$.

We identified all particles enclosed within a sphere of radius $`2h_0^1\mathrm{Mpc}`$ around the supergalactic center at the present time as shown in the left panel of Figure 4). We have found 11 such particles (shown as the open circles), corresponding to a total mass of $`3.9\times 10^{12}h_0^1M_{}`$. For comparison, the Local Group has a zero-velocity surface (separating it from the Hubble expansion) of $`1.18\pm 0.15\mathrm{Mpc}`$ and a total mass of $`M_{\mathrm{LG}}=(2.3\pm 0.6)\times 10^{12}M_{}`$ (Courteau & van den Bergh, 1999; van den Bergh, 1999) based on the radial velocity dispersion of its members. Initially, the sampled 11 particles have a mean peculiar velocity of $`670\mathrm{km}\mathrm{s}^1`$ at $`t=t_0`$ towards the direction of $`(250,502,368)\mathrm{km}\mathrm{s}^1`$ with respect to the cosmic microwave background (CMB) frame. For comparison, the observed motion of the Local Group relative to the CMB frame is $`627\pm 22\mathrm{km}\mathrm{s}^1`$ in the direction of $`(406,352,324)\mathrm{km}\mathrm{s}^1`$ (Kogut et al., 1993). The agreement between the simulated and the observed velocity amplitude at $`t=t_0`$ is impressive given the crudeness of the initial conditions of the simulation; however, the direction of the velocity vector deviates somewhat from the observed orientation.

More recently, Tonry et al. (2000) used the results from the Surface Brightness Fluctuation Survey to model the local and large-scale flows. Their best-fit model includes two attractors, one of which having a best-fit location coincident with the Virgo cluster and the other having a fit location slightly beyond the Centaurus cluster (which is commonly referred to as the Great Attractor). Assuming an attenuated power-law mass distributions for the two attractors, their best-fit model has enclosed mass (in excess of background density) of $`7\times 10^{14}M_{}`$ for Virgo and $`9\times 10^{15}M_{}`$ for the Great Attractor within spheres centered on the attractors with radii reaching the Local Group. For comparison, the corresponding masses in our simulation are $`7\times 10^{14}h_0^1M_{}`$ for the Virgo and $`8\times 10^{15}h_0^1M_{}`$ for the Centaurus, in good agreement with the observational estimate by Tonry et al.

We follow the trajectories of the above 11 particles up to the time $`t=t_0+6t_H`$. The right panel of Figure 4 indicates that the particles have traversed a comoving distance of $`5h_0^1\mathrm{Mpc}`$ towards the direction of the Virgo cluster but have not fallen into it. The physical distance to Virgo (which is proportional to $`a(t)`$) would increase exponentially at late times. Neighboring particles show a similar behavior, so this conclusion is not likely to be affected by the random sampling of particles from the original data set described in § 2. By $`t=t_0+6t_H`$, the peculiar velocity field at moderate overdensities is heavily damped (see § 3). Most of the displacement towards Virgo takes place by $`t=t_0+2t_H`$. Note that the Virgo cluster itself has also moved by a few comoving $`h_0^1\mathrm{Mpc}`$ in response to the pull of the Great Attractor. We conclude that the Local Group will not be gravitationally bound to the Virgo cluster, given the initial conditions of the simulation that are matched to the observed galaxy distribution.

### 5.2 The Milky Way and the Andromeda Galaxies

The Local Group is a centrally concentrated system in which half of all members lie within $`450\mathrm{kpc}`$ of its center-of-mass. Most of the mass is concentrated in the Andromeda and the Milky Way subgroups. Unfortunately, the mass estimate of each subgroup suffers from large uncertainties. For example, the mass of the Andromeda subgroup ranges from a value of $`M_\mathrm{A}=(13.3\pm 1.8)\times 10^{11}M_{}`$ (Courteau & van den Bergh, 1999) to a lower estimate of $`7.0_{3.5}^{+10.5}\times 10^{11}M_{}`$ (Evans et al., 2000). The mass of the Milky Way ranges from $`M_{\mathrm{MW}}=(8.6\pm 4.0)\times 10^{11}M_{}`$ (Zaritsky, 1999) to a larger value of $`19_{17}^{+36}\times 10^{11}M_{}`$ (Wilkinson & Evans, 1999). All of the existing estimates are consistent with each other to within the quoted errors. Under most estimates, the total mass of the Local Group is close to $`M_{\mathrm{LG}}=(23\pm 6)\times 10^{11}M_{}`$ (van den Bergh, 1999).

We may evaluate the present-day overdensity of the Local Group, $`\delta _{\mathrm{LG}}`$, using its total mass and the distance of $`\mathrm{}=740\pm 40\mathrm{kpc}`$ (Binney & Merrifield, 1998) between the Milky Way and Andromeda, yielding $`1+\delta _{\mathrm{LG}}=3M_{\mathrm{LG}}/[4\pi (\mathrm{}/2)^3\overline{\rho }]=130h_0^2`$. The present overdensity of the Local Group clearly exceeds the critical collapse threshold of $`\delta _c=17.6`$.

Adopting the measured radial component of the relative velocity between the Andromeda and the Milky Way galaxies $`v_r=120\mathrm{km}\mathrm{s}^1`$ (Binney & Tremaine, 1987) and ignoring its tangential component, it can easily be shown that the total energy of the Local Group, $`\frac{1}{2}\mu _{\mathrm{LG}}v_r^2GM_\mathrm{A}M_{\mathrm{MW}}/\mathrm{}`$, is negative, where $`\mu _{\mathrm{LG}}=M_\mathrm{A}M_{\mathrm{MW}}/(M_\mathrm{A}+M_{\mathrm{MW}})`$ is the reduced mass of the Andromeda and the Milky Way subgroups. This result is obtained for either of the latter or the former pair of mass values quoted above for $`M_\mathrm{A}`$ and $`M_{\mathrm{MW}}`$. The transverse component of the relative velocity is commonly assumed to be much smaller than the radial velocity (Einasto & Lynden-Bell, 1982). Based on the radial separation and the relative velocity between the two galaxies, one finds that they are likely to merge within a Hubble time. A pedagogical discussion on the dynamical evolution of the Local Group based on the least action principle can be found in Peebles (1993). Dynamical evaporation of stars out of the merger product can be neglected on the time-scale of interest here (Binney & Tremaine, 1987).

## 6 Conclusions & Discussions

We have simulated the future evolution of our cosmic neighborhood in a universe dominated by a cosmological constant using an N-body code and the initial conditions that were reconstructed from the observed galaxy distribution of the IRAS 1.2 Jy survey. We find that the large-scale structure and the mass distribution of bound objects will freeze in $`30`$ billion years from today due to the accelerated expansion of the Universe. The Local Group of galaxies will get somewhat closer to the Virgo cluster of galaxies in comoving coordinates, but will rapidly recede from Virgo in physical coordinates. However, the overdensity inside the Local Group is well above the required threshold ($`\delta _c=17.6`$) for it to resist the repulsive gravitational force of the cosmological constant. Therefore the Milky Way and the Andromeda galaxies are likely to merge within a Hubble time.

If the Universe is dominated by a cosmological constant, then the Local Group is detached from the rest of the Universe and the physical distance from us to all other systems that are not bound to the Local Group will increase exponentially with time in the distant future. Combining our simulations with the results of Loeb (2002), we predict that when the age of the Universe will be $`100`$ billion years, there will only be one massive galaxy for us to observe, namely the merger product of the Andromeda and the Milky Way galaxies. All other systems that are not gravitationally bound to the Local Group will exit through our event horizon and their images will fade rapidly while remaining frozen on the time of their exit.

The precision of our quantitative results is limited by the accuracy of the initial conditions of the simulation at $`z=0`$ compared to the true mass density field in the nearby universe. The locations of the simulated clusters of galaxies are shifted by a few megaparsec relative to the recent observational estimates, and in principle the match between the two can be improved in the future by refining the reconstruction method of the initial conditions of the simulation from the observational data. While we believe that the general predictions derived in this paper are independent of the limited mass and spatial resolution of our simulations, it would be appropriate to perform in the future higher-resolution simulations with improved initial conditions that match the observed galaxy distribution. Since the observational uncertainties in the distance measurements are still substantial, the precision of the observed galaxy distribution will also get better with future observations. The combination of both of the above improvements will allow us to make more accurate predictions for the future evolution of the nearby large-scale structure.

Acknowledgments
We are grateful to the authors of Mathis et al. (2002) paper and the GIF (German-Israel Foundation) collaboration for allowing us to use their $`z=0`$ simulation output as the initial condition of our simulation. In particular, we thank Volker Springel for providing us with the data in a convenient format for the GADGET code, as well as for clearing out some of our technical problems and questions at the initial stage of this work. This work was supported in part by the grants 5-7768 from NASA and AST-9900877, AST-0071019 from NSF for AL.