text:
# Measurement of Electron Backscattering in the Energy Range of Neutron 𝛽-Decay

## I Introduction

Backscattering of electrons from the surfaces of bulk materials has been studied at low energies largely in relation to materials science applications (e.g. Auger electron spectroscopy and scanning and transmission electron microscopy). For this reason, most detailed studies have been conducted for incident electron energies less than 40 keV Gérard et al. (1995); Massoumi et al. (1992, 1993). At higher energies ($`E>>`$ 1 MeV) sophisticated Monte Carlo calculations exist Agostinelli et al. (2003) that are constrained to reproduce data obtained for nuclear or particle physics applications. However in the intermediate regime ($`0.04\mathrm{MeV}<E<1\mathrm{MeV}`$) there is little data to constrain the simulations. This intermediate energy regime can be important for various beta spectroscopy applications, in particular neutron beta-decay ($`E_{max}=0.782\mathrm{MeV}`$).

Most of the data that exists in this energy region is based on measurements using electron beams and detecting the backscattering from bulk targets using electrical currents in Faraday cups (for a recent review, see Ref. T. Tabata and Shinoda (1999)). In particular, for carbon and aluminum targets, extensive measurements using this technique exist. We have made measurements for beryllium and silicon targets using this technique, and have investigated the effects of secondary electron emission on these types of measurements. In addition, we have conducted measurements using ion-implanted silicon detectors to detect the energy and angle of backscattered electrons, extending the work of Refs. Gérard et al. (1995); Massoumi et al. (1992, 1993) to higher energies. In this way, we can evaluate the reliability of existing models to reproduce the dependence of the backscattering process on energy and emission angle of the backscattered electrons.

The results of these measurements can have a variety of nuclear and particle physics applications. For example, an upcoming measurement of the electron correlation in polarized neutron decay Bowles and Young co-principal investigators (2000) will require an understanding of electron backscattering at the 10-20% level. In this experiment plastic scintillator and silicon multi-strip detectors will be used, thus requiring detailed information on backscattering for low $`Z`$ targets. As a first step in addressing this problem, we have carried out detailed measurements of the energy and angular distributions of backscattered electrons from light materials (Be and Si) for incident electron energies between 43.5 and 124 keV.

## II Experiment Overview

The experiment consisted of an electron gun and a scattering chamber containing a movable target. Two modes of running were used where backscattering data could be acquired. In one mode the silicon detector was used to detect the energy and angle of backscattered electrons. In a second higher-current mode the electrical current due to the backscattered electrons incident on the chamber walls was detected. These two modes will be referred to as silicon detector mode and current integration mode, respectively. Each component of the experiment will be described in greater detail.

### II.1 The Electron Gun

The Kellogg Electron Gun was constructed in order to perform these measurements. A schematic of the gun is shown in Fig. 1.

Electrons were extracted from a hot filament by a potential of typically 6 kV, in the source head. The electrons were focused by an Einzel lens before passing through an additional acceleration potential of up to 130 kV minus the extractor voltage.

The electron beam current could be varied from a few electrons per second to several $`\mu `$A, by adjusting the current passed through the filament. The resultant beam current remained stable after allowing the filament to warm up for typically half an hour.

After the acceleration column, two magnetic coils allowed steering of the beam in the horizontal and vertical directions.

A Faraday cup could be inserted into the beam to monitor the total current of the beam. A beam position monitor, consisting of a rotating wire which sweeps through the beam twice along two orthogonal axes, could be used to locate the beam within the beam pipe.

The energy of the electron gun was stable and reproducible. The energy was absolutely determined to the 1% level using a novel Helmholtz coil spectrometer. This spectrometer is iron-free and has a momentum resolution of 0.3% Yuan et al. (2001). The energy of the electron gun was found to be monochromatic to better than the resolution of the spectrometer.

### II.2 Chamber, Targets and Detector

A schematic of the chamber and detector arrangement used to perform the backscattering measurements is shown in Fig. 2.

The chamber contained a target and a silicon detector on separately rotatable feedthroughs. The chamber, target, and detector were each separately isolated from ground. The electric currents detected by the chamber and the target were integrated using two Ortec model 439 current digitizers.

The targets The Goodfellow Corporation, 800 Lancaster Avenue, Berwyn, PA 19312-1780  were multiple stopping lengths of beryllium and silicon. Two different beryllium targets of nominal thicknesses 1.0 mm and 2.0 mm were used. Two different silicon targets of thickness 1.0 mm were used. The targets were 25 mm by 25 mm square, and were held at one corner in a slotted rod by a set screw. The beryllium targets were nominally 99.0% pure, and had a metallic appearance and the polycrystalline silicon target was nominally 99.999% pure.

The energy and angular distribution of the backscattered electrons was measured using an Ortec ion-implanted silicon detector (Ortec model BU-13-25-300, 25 mm2 nominal active area and 300 $`\mu `$m nominal depletion depth at 100 V bias). The detector was mounted on a rotatable arm which allowed it to be placed at an arbitrary angle with respect to the target. The target to detector distance was typically 8.5 cm. The detector linearity and resolution was calibrated by placing the detector in the electron beam at very low beam current ($``$ 100 Hz). The detector linearity was confirmed to the sub-kilovolt level for these experiments. The energy resolution of the silicon detector was found to be typically $`\sigma =2.5`$ to 3.5 keV, and the resolution was independent of energy to less than 0.1 keV.

For fixed filament current, the current detected on the target was maximized in order to tune the beam. The beam tunes for each energy were stored on the computer that controlled the power supplies of the accelerator. In addition, tunes were checked using a scintillator target coated with graphite which produced a visible spot ($``$ 3 mm $`\times `$ 3 mm) indicative of the beam spot size when struck with sufficient current (typically 1 nA).

## III Backscattering Measurements

Backscattering measurements were performed for normal incidence upon the target. Measurements were performed for incident electron energies of 43.5, 63.9, 83.8, 104, and 124 keV. For each energy, both silicon detector mode and current integration mode measurements were taken.

### III.1 Silicon Detector Mode

In the silicon detector mode, silicon detector spectra were acquired for backscattered electron angles of 20 to 80 degrees in steps of 10 degrees. The backscattered angle $`\theta `$ was defined with respect to the normal of the target, as shown in Fig. 2.

When the beam was off, the detector rate was typically 100 Hz, due to low energy noise. The detector rate with beam on was typically 5 kHz, and was always kept below 20 kHz, to limit ADC pile-up to below 3% on average across the spectrum. The beam current for each target and detector angle was therefore different, due to the different levels of backscattering. For Si, the current was typically 40 pA, while for Be, it was 300 pA.

Fig. 3 shows normalized spectra taken using the silicon detector for a variety of detector angles $`\theta `$ for 124 keV electron beam energy and the beryllium target.

The data are plotted as a function of the digitized pulse height (ADC channel), which is proportional to the energy $`E`$ deposited in the silicon detector. The data are not corrected for the response of the silicon counter, which had a resolution of about 2.5 keV, and itself suffered from backscattering. Below about 20 keV, contributions from noise in the detector and electronics contributed to the count rate at the few percent level, and so no data is shown beneath that threshold.

A detector rate above electronic noise was also detected for angles $`\theta >90^{}`$, where the detector is shielded from any direct backscattering. It can be seen that the rates at $`100^{}`$ and $`140^{}`$ are significant, especially compared to e.g. the $`80^{}`$ measurements. The flux found at $`\theta >90^{}`$ was determined to be due to rescattering from the walls of the chamber, and was confirmed in Monte Carlo studies. In these studies, backscattering from the steel walls and lid of the chamber were included. The chamber floor, dominantly aluminum, had relatively little effect on the chamber background and was not included. The Monte Carlo calculations showed that electrons at $`\theta >90^{}`$ could be produced, and that the result was numerically about 50% larger than the observed chamber background. The disagreement is attributed to the use of a simplified geometry in the Monte Carlo, and uncalibrated treatment of backscattering from these sources.

Motivated by the simulation, experimental studies were carried out lining the chamber walls with stopping thicknesses of aluminized mylar (which has a much lower backscattering fractions). This was found to reduce the chamber background by roughly 50%, a figure in agreement with expectations based on the Monte Carlo.

Based on the Monte Carlo simulations, the $`100^{}`$ data is expected to represent most closely the chamber background for the $`80^{}`$ data, and to represent an upper bound on the chamber background for the other angles. So the $`100^{}`$ data was subtracted from the $`\theta =20^{}80^{}`$ data to arrive at the background-subtracted spectrum. The systematic error in using this subtraction scheme is 3-5% for 124 keV beam energy, peaked at low energy deposition, and is smaller for lower beam energy (as more of the chamber background is moved into the unmeasured region dominated by detector noise).

From Monte Carlo studies and analytical estimations, it was also determined that contributions from X-rays were negligible in both flux (compared to electron backscattering) and in detection efficiency (for the 300 $`\mu `$m thick detector used). Therefore no corrections were made for X-rays.

### III.2 Systematic uncertainties for silicon detector measurements

Table 1 summarizes the systematic uncertainties for the silicon detector measurements.

Each effect in the table will now be described in more detail.

Measurements of the backscattered electron yield in the silicon detector were found to be reproducible at the 7% level. The detector active area was measured using an alpha source and various collimators, and found to be consistent with the geometry of the detector to about 4%. The beam spot size and detector size were also measured by scanning the beam across the detector. The detector size was again found to be in agreement with the directly measured value. At the same time, it was observed that no significantly different response for electrons could be seen for the detector as one varied the position of the beam on the detector surface. A simple silicon response function including only normal-incidence backscattering was therefore used to convolute Monte Carlo predictions for comparison with the data (see section V). The beam spot radius was found to be about 1.7 mm, and the beam was always well-centered on the target to this level, limiting a solid angle correction from finite spot size to less than $`5\%\times \mathrm{sin}\theta `$. Deadtime corrections varied between $`<10\%`$ and $`50\%`$. However, ADC pile-up was kept below 3% and the deadtime corrections could be performed reliably at that level, using information from fast scalers counting the triggers.

The relative alignment of the beam, target, and detector in the plane of rotation of the detector was measured to the level of 0.5 degrees. This was confirmed by taking measurements of electron backscattering at both positive and negative $`\theta `$. The residual systematic error due to alignment problems was estimated to be about 2%.

The silicon measurements were normalized to the total current detected by the chamber and target. The detection of current was calibrated using a precision picoampere source, and comparison among several well-calibrated picoammeters. The calibration agreed each time at the level 0.3 pA, giving a worst case contribution to the normalization uncertainty of 3%. Integration of low-energy secondary electrons was not a factor for this measurement of current, since the total charge on the chamber and target was collected.

The total normalization systematic uncertainty after background subtraction was therefore 12%, dominated by reproducibility and detector solid angle uncertainties.

### III.3 Current Integration Mode

For each energy, current integration mode measurements were also performed. By separately sensing the currents on the target ($`I_{\mathrm{target}}`$) and the chamber ($`I_{\mathrm{chamber}}`$), the total normal-incidence backscattered fraction $`\eta `$ was determined via:

$$\eta =\frac{I_{\mathrm{target}}}{I_{\mathrm{target}}+I_{\mathrm{chamber}}}.$$
(1)
Low energy secondary electrons may be emitted from the surfaces of materials after higher energy electrons strike the surface. In previous experiments, secondary electrons were typically defined to have $`E<50`$ eV. If the fraction of secondary electrons produced per high-energy interaction is large, the sensed currents will give an erroneous measure of the high-energy backscattered fraction.

To suppress and quantify the effects of secondary electrons, a cylindrical “cage” of wires (referred to as the grid) was inserted into the setup to provide a potential wall between the chamber and the target. This grid was made of 50 $`\mu `$m tungsten wire, wrapped on a cylindrical copper frame of radius 6 cm and height 8 cm. The grid had 22 vertical wires evenly spaced on the sides. The wires met at the bottom, but had an opening at the top so that the target rod (which holds the target) could be inserted. The grid, target, and chamber could each be biased at different voltages up to a difference of 200 V with leakage currents between elements kept to below 1 pA.

### III.4 Systematic Uncertainties for Current Integration Measurements

Table 2 summarizes the systematic uncertainties encountered for the current integration measurements.

Each of these uncertainties will now be described in more detail.

It was found that biasing the grid at roughly -50 V relative to the chamber and target caused secondaries created on the chamber and target to be recollected by the chamber and target. When the target was biased to +50 V, and the chamber held at ground, the effect of changing the grid voltage from zero to -50 V resulted in changes of 10-30%, depending on the target material and the beam energy. This agreed well with estimates of the effect based on the Penelope Monte Carlo code, and based on measurements of secondary electron emission summarized in Refs. Joy (1995); Bronshtein and Denisov (1965b, a).

Due to a small piece of the conducting target rod (held at the same potential as the chamber) penetrating into the top of the grid, there was a residual correction still to be made for secondaries. In the data, this showed up as a residual dependence of $`\eta `$ on the relative target/chamber voltage, even when the grid was set to very large voltages. It was found that a correction could be made using the solid angle subtended by that piece of the target rod, and the value of $`\eta `$ determined when the grid was not used. This correction amounted to a 7% contribution to the systematic uncertainty in the determination of $`\eta `$.

Secondaries due to high-energy electrons striking the grid could be accounted for at the 1% level. Measurements of $`\eta `$ under widely varying beam conditions showed that it was reproducible at the 5% level.

The current dependence of $`\eta `$ was also studied from pA to $`\mu `$A in incident beam current, and found to vary at the 3% level. This indicated that electrostatic charging of various non-conducting components in the setup contributed at a small level.

The total fractional systematic uncertainty in $`\eta `$ determined from these current-mode measurements was therefore 9%.

## IV Results

### IV.1 Silicon Detector Mode

The normalized, background-subtracted spectra accumulated for various detector angles for 124 keV electrons normally incident on silicon and beryllium bulk targets are shown in Fig. 4.

The Monte Carlo curves will be described in Section V. The data are plotted as a function of the dimensionless energy $`q=E/E_{\mathrm{beam}}`$, where $`E`$ is the energy detected by the silicon detector and $`E_{\mathrm{beam}}`$ is the energy of the incident electrons, in this case 124 keV. On the vertical axis, $`\frac{1}{N_e}\frac{dN}{dqd\mathrm{\Omega }}`$, the number of counts per incident electron, per unit $`q`$, per unit solid angle is plotted. In the absence of the effects of detector response (resolution and backscattering), this would be the normal-incidence backscattered fraction per unit $`q`$, per unit solid angle.

In Fig. 4a., for silicon, for small backscattered angles ($`\theta =20^{}30^{}`$), a peak is found near $`q=0.65`$, and a shoulder found near $`q=0.95`$. As the backscattered angle increases, the peak at $`q=0.65`$ tends to disappear and shift slightly to higher $`q`$, while the shoulder at $`q=0.95`$ tends to become more pronounced as events at lower $`q`$ disappear.

The same trends can be seen in the spectra for a beryllium target, shown in Fig. 4b. However, in this case, the low-energy peak appears closer to $`q=0.35`$. This can be explained by the beryllium having smaller Rutherford scattering cross section. Electrons therefore penetrate more deeply into the material before scattering.

The silicon results compare well qualitatively with the results of Refs. Gérard et al. (1995); Massoumi et al. (1992, 1993), which were acquired at lower energies on silicon and aluminum targets. It is difficult to compare directly with these measurements. For Ref. Massoumi et al. (1992, 1993), no data is published above 35 keV and more data are displayed for oblique incidence. Ref. Gérard et al. (1995), while providing data at 40 keV, gives no absolute normalization.

The dependence on beam energy was also investigated. Fig. 5 shows the dependence of the backscattered fraction on both the beam energy and the energy of the backscattered electron, integrated over all possible backscattered angles.

The integration over angle was performed using a finite sum with $`10^{}`$ bins centered on each measured angle, and the appropriate solid angle weighting. A small correction due to the unmeasured regions at small and large angles was included. The systematic uncertainty in approximating the integral by the sum was typically 4%, from comparison with analytical forms. Angular bin centering corrections were found to be negligible, due to the smoothness of the angular behavior. The overall systematic uncertainty was consequently increased to 13%.

As seen in Fig. 5, when plotted in terms of the dimensionless variable $`q`$, the curves nearly overlap. The same qualitative behavior can be seen as a function of beam energy for individual backscattered angles. The data are therefore observed to follow a near-scaling behavior. The overall normalization follows the energy dependence of the total normal-incidence backscattered fraction.

The angular dependence of the missed backscattered fraction can be determined by integrating over energy the data of Fig. 4. The result is shown in Fig. 6.

A linear fit based on the first 20 keV of data above the analysis cut was used to extrapolate to 50 eV (the defined threshold for secondary electrons), so that these integrals and subsequent integrals could be compared with the current integration measurements. In order to estimate the reliability of the fit, an additional systematic uncertainty was assigned to the extrapolation, based on comparison of this fit to a fit constrained to pass through zero at zero energy. For 124 keV beam energy, this extra systematic uncertainty was of order a few percent.

In order to better estimate the systematic uncertainty in this extrapolation, simulations were performed using Geant4 and Penelope in the unmeasured region 50 eV to 20 keV. In the Geant4 simulation, the backscattered fraction was found to tend toward zero at small backscattered electron energy, in fair agreement with the linear extrapolation method. In the Penelope simulation, the backscattered fraction was found to rise steadily as the threshold of the simulation was reduced. In the range 100 eV to 10 keV, Penelope gave about 10% extra contribution to the integral compared to linear extrapolation, for beryllium, and an extra 4% for silicon. To average between the extrapolations implied by Geant4 and Penelope, an additional 5% contribution was added to the beryllium data, and an additional 2% contribution was added to the silicon data. The additional uncertainty in each case was taken to be the size of the correction.

The data were integrated over angle, using the same method described in relation to Fig. 5, to determine the total normal-incidence backscattered fraction. The results of this integration are shown in Fig. 7 and are compared with current integration measurements (described in section IV.2) and Monte Carlo simulations (described in Section V).

In Fig. 7, the total systematic uncertainty, including extrapolation to 50 V and extrapolation over unmeasured angles, is plotted.

### IV.2 Current Integration Mode

The results for $`\eta `$ based on our current integration measurements are also shown in Fig. 7.

The silicon detector measurements are found to be systematically higher than the current-mode measurements; however, the two methods of are found to agree within the systematic uncertainties. In the case of the current integration method, this systematic uncertainty is dominated by residual correction for secondary electron collection due to the penetration through the grid of the target rod. In the case of the silicon detector measurements, it is dominated by reproducibility of the measurements under varying conditions, and by uncertainties in alignment and solid angle effects.

The data are also compared with previous data on Be and Si targets due to Drescher et al. H. Drescher and Seidel (1970) and with data on Be due to Neubert et al. Neubert and Rogashewski (1980). Both groups used current integration techniques to arrive at their results. Neubert et al. Neubert and Rogashewski (1980) in particular used a second target apparatus to study the effects of secondary electrons, as opposed to the grid used in this work. Only the subsets of their data that overlap the region 43.5 to 124 keV are plotted. The Drescher data on Be are systematically higher than the Neubert data. However, due to the lack of additional data on this element, it is impossible to say which is more accurate. Our data tend to agree with the Neubert data, as do the data of Massoumi et al. Massoumi et al. (1992, 1993) taken below 40 keV.

Our data on Si is in good agreement with those of Drescher, and tend to suggest that there could be some systematic effect in either the current integration or silicon detector data for that element. Given the current level of the systematic uncertainties, it is difficult to make a firm statement.

## V Comparison with Monte Carlo

Monte Carlo simulations were conducted using the Geant4 Monte Carlo Agostinelli et al. (2003) and the Penelope Monte Carlo Sempau et al. (1997).

The version of Geant4 used was 4.4.0. In order to achieve a reasonable description of backscattering, it was found that three parameters had to be changed from their default values: the maximum step size, the threshold to create secondaries, and a parameter originally introduced into Geant4 to tune low-energy EM processes. For all our simulations, the maximum step size was set to the range of a 1 keV electron; the threshold to create secondaries was set to 1 $`\mu `$m; and the tuning parameter was set to zero. The step size and threshold parameters were chosen by reducing them until the integral backscattered fraction was relatively stable under variation of those parameters, and to be as small as possible for reasonable running time. The threshold parameter additionally was checked to give good results at low energy for thin targets. The tuning parameter had to be changed to zero, as it had been found to be erroneously set to 1.5 in this version of Geant4. Examples of how to make these modifications were supplied by the Geant4 Electromagnetic Physics and Low-Energy Electromagnetic physics groups Urban and Ivantchenko (2003).

The version of Penelope used was 2002b. Penelope was studied in detail under the variation of several simulation parameters and was found to be stable. The simulation parameters were therefore chosen to optimize simulation speed consistent with a full detailed simulation. For the comparison with our silicon detector data, the most suitable simulation parameters were found to be: $`E_{abs}=10\mathrm{keV}`$, $`W_{cc}=W_{cr}=5\mathrm{keV}`$, DSMAX = 0.005 cm, and $`C_1=C_2=0.05`$ Hoedl (2003). These parameters control energy cut-offs, the maximum step size, and the description of elastic and inelastic scattering in the medium. Particles identified as secondaries by Penelope were included. Secondaries were also studied in separate simulations related to our current-integration measurements, as mentioned earlier.

In both Monte Carlo models, backscattering is defined as any electron which exits the surface of the target. The Monte Carlo simulations included silicon detector response in a simple model including normal-incidence backscattering from the front face. This was validated by the fact that the detector response to mono-energetic electrons did not vary with beam position on the detector face, as discussed in Section III.2. Energy smearing based on a 2.5 keV energy-independent energy resolution of the silicon detector was also included. Chamber background effects were studied in separate Penelope-based simulations, as described earlier.

The main difference between Penelope and Geant4 relevant to backscattering is in the treatment of scattering from nuclei. Penelope treats these as Rutherford scattering events exactly using a relativistic screened Rutherford scattering cross section. Multiple electron-atom scattering events, with energy loss (dominated by electron-electron interactions in the target) are found to dominate backscattering at these energies in this model, for thick targets. Geant4, on the other hand, has no exact treatment of Rutherford scattering, relying on sampling from a multiple scattering distribution.

Fig. 4 compares the Geant4-based Monte Carlo with our silicon detector measurements at 124 keV. It can be seen that Geant4 somewhat overestimates the beryllium data, while having relatively good agreement with the silicon data. The beryllium data is globally overestimated by roughly 20%, hence a scale factor of 0.8 was applied to the Monte Carlo for this comparison so that the differences between the distributions can be more easily seen. In all cases, it is also apparent that the peak near $`q=0.95`$ is systematically underestimated. However the positions of the low-energy and elastic peaks are rather well-described by the Monte Carlo. In the case of a silicon target, the Geant4 low-energy peak is at a slightly higher $`q`$ than the data, while for beryllium this peak appears at slightly lower $`q`$ relative to the data.

Fig. 8 compares Penelope-based Monte Carlo with our measurements using the silicon detector.

A scale factor of 1.1 is applied to the Penelope simulation in the case of silicon. No scaling is applied to the beryllium simulation. The Penelope simulation in general somewhat underestimates the beryllium measurements by 5% and underestimates the silicon measurements by roughly 10%. When the Monte Carlo is rescaled, it is apparent that trends in both energy and angle are well represented by Penelope. Since the systematic uncertainty is larger than these discrepancies, the overall agreement is good.

Fig. 6 shows the Geant4 and Penelope simulations compared with data for the energy-integrated angular distributions ($`d\eta /d\mathrm{\Omega }`$). No Monte Carlo scale factors are applied for this comparison. The same discrepancies in overall magnitude of backscattering can be seen. As noted previously, the Penelope simulation tends to better describe the angular distribution (aside from the overall scale factor). The Geant4 distributions are somewhat narrower compared to the data and Penelope. Additionally, the Geant4 simulation gives systematically larger backscattering from each material than does the Penelope simulation.

Fig. 7 compares the integrated $`\eta `$ results for different beam energies with Geant4 and Penelope simulations. The same discrepancies in normalization are again observed. Both Penelope and Geant4 adequately describe the reduction of $`\eta `$ as the beam energy increases.

Comparisons of other Monte Carlo models to Geant4, Penelope, and the existing backscattering data were also carried out and were reported elsewhere Hoedl (2003).

## VI Conclusion

A detailed data set for normal incidence backscattering from beryllium and silicon bulk targets has been acquired for incident electron energies from 43.5 to 124 keV. Two methods of determining the total normal-incidence backscattered fraction were compared, and found to agree within systematic uncertainties. The data agree qualitatively with previous measurements of the doubly-differential distributions of backscatter done at lower energy. They also agree quantitatively with previous measurements using current integration techniques. The data agree well with models implemented in the Geant4 and Penelope Monte Carlo codes. In terms of overall normalization, Geant4 is found to give a good description of the data. Penelope is found to give good agreement in terms of both overall normalization and relative distributions of backscattered electrons in angle and energy.

In future measurements, we plan to investigate non-conducting targets such as scintillator (particularly interesting for a variety of nuclear physics applications) and active targets (such as scintillator and silicon detectors). We also plan to extend this data set to higher energy (up to 1 MeV) using the same setup with a higher energy electron source. This will allow further tests of the available Monte Carlo models. The comparison of these data to Monte Carlo models will provide an important constraint on a part of the systematic uncertainty for future high precision measurements of neutron beta decay. Future work will also include investigation of thin films and oblique-incidence backscattering which are both important in understanding systematic effects in these experiments.

We gratefully acknowledge the technical support of Robert Carr. This work was supported by the National Science Foundation.


format: text

text:
# Axisymmetric MHD simulations of the collapsar model for gamma-ray bursts

## 1 Introduction

The collapsar model is one of most promising scenarios to explain the huge release of energy in a matter of seconds, associated with gamma-ray bursts (GRBs; Woosley 1993; Paczyński, 1998; MacFadyen & Woosley 1999, hereafter MW; Popham, Woosley & Fryer 1999; MacFadyen, Woosley & Heger 2001). In this model, the collapsed iron core of a massive star accretes gas at a high rate ($`1\mathrm{M}_{}\mathrm{s}^1`$) producing a large neutrino flux, a powerful outflow, and a GRB. Despite many years of intensive theoretical studies of these events, basic properties of the central engine are uncertain. In part, this is because previous numerical studies of the collapsar model did not explicitly include magnetic fields, although they are commonly accepted as a key element of accretion flows and outflows.

In this letter we present a study of the time evolution of magnetohydrodynamic (MHD) flows in the collapsar model. This study is an extension of existing models of MHD accretion flows onto a black hole (BH; Proga & Begelman 2003, PB03 hereafter). In particular, we include a realistic equation of state (EOS), photodisintegration of bound nuclei and cooling due to neutrino emission. Our study is also an extension of MW’s collapsar simulations, as we consider very similar neutrino physics and initial conditions but solve MHD instead of hydrodynamical equations.

## 2 Method

To calculate the structure and evolution of an accreting flow, we solve the equations of MHD:

$$\frac{D\rho }{Dt}+\rho 𝐯=0,$$
(1)
$$\rho \frac{D𝐯}{Dt}=P\rho \mathrm{\Phi }+\frac{1}{4\pi }(\times 𝐁)\times 𝐁,$$
(2)
$$\rho \frac{D}{Dt}\left(\frac{e}{\rho }\right)=P𝐯+\eta _r𝐉^\mathrm{𝟐},$$
(3)
$$\frac{𝐁}{t}=\times (𝐯\times 𝐁\eta _r𝐉),$$
(4)
where $`\rho `$ is the mass density, $`P`$ is the total gas pressure plus radiation pressure, $`𝐯`$ is the fluid velocity, $`e`$ is the internal energy density, $`\mathrm{\Phi }`$ is the gravitational potential, $`𝐁`$ is the magnetic field vector, $`𝐉`$ is the current density, $`\eta _r`$ is an anomalous resistivity, and $``$ is the cooling rate due to neutrinos.

To compute resistivity, we follow Stone & Pringle (2001, see their equations 5 and A1). We perform simulations using the pseudo-Newtonian potential of the central mass $`\mathrm{\Phi }_{pw}=GM/(rR_S)`$, where $`R_S=2GM/c^2`$ is the Schwarzschild radius, introduced by Paczyński & Wiita (1980). We increase the mass of the BH during the calculation by the amount of baryonic rest mass accreted through the inner boundary.

Our calculations are performed in spherical polar coordinates $`(r,\theta ,\varphi )`$. We assume axial symmetry about the rotational axis of the accretion flow ($`\theta =0^{}`$ and $`180^{}`$). The computational domain is defined to occupy the radial range $`r_i=1.5R_Srr_o=1000R_S`$, and the angular range $`0^{}\theta 180^{}`$. The $`r\theta `$ domain is discretized on a non-uniform grid as in PB03, which yields $`\mathrm{\Delta }r/r=0.278`$ at the inner edge of the simulations.

We adopt a realistic EOS, which includes contributions from an ideal gas of nuclei, radiation, and electrons and positrons with arbitrary degrees of relativity and degeneracy (Blinnikov, Dunina-Barkovskaya & Nadyozhin 1996). To compute the cooling rate, we follow Itoh et al. (1989, 1990), taking into account thermal neutrino emission processes dominated by pair annihilation, as well as the capture of pairs on free nucleons.

Our calculations use the ZEUS-2D code described by Stone & Norman (1992a,b). We have extended the code to include the realistic EOS, artificial resistivity, photodisintegration and neutrino cooling. We included the terms due to the resistivity and cooling in an operator split fashion separately from the rest of the dynamical equations. For stability, these terms must be integrated using the time step computed based on the resistive and cooling time scales, respectively. We subcycle whenever either of these two time steps is smaller than the time step used in the MHD equations (see, e.g., Stone, Pringle, and Begelman 1999). Inclusion of a non-adiabatic EOS requires iterating over temperature, $`T`$ (see MW for details).

By including the neutrino cooling and realistic EOS, we consider very similar microphysics to that used in MW. Our simulations differ from those in MW in that we use the ZEUS MHD code whereas MW used the PROMETHEUS hydrodynamics code (Fryxell, M$`\ddot{\mathrm{u}}`$ller & Arnett 1989). This means that we can self-consistently calculate turbulent stresses generated by the magnetorotational instability (MRI) and thus include the outward transport of energy and angular momentum (e.g., Balbus & Hawley 1991). MW implemented the effects of viscosity in the disk using an alpha viscosity as prescribed by Shakura & Sunyaev (1973). We do not consider self-gravity and nuclear burning; however, as in MW our simulations track the photodisintegration of helium. Our simulations span the radial direction from 1.5 to 1000 $`R_S`$ whereas MW’s simulations span from 9.5 to 9500 $`R_S`$. Thus we can capture the innermost part of the flow near the BH but follow the evolution out to a smaller radius than MW’s simulations.

We adopt PB03’s boundary conditions and initial conditions for the magnetic field. In particular, the initial magnetic field is purely radial and weak ($`\beta 8\pi P/B^21`$ everywhere). For the initial conditions of the fluid variables, we follow MW and adopt the stellar model for the helium core of a 25 $`\mathrm{M}_{}`$ presupernova star (model S251S7B@14233 in Woosley & Weaver 1995). The masses of the helium and iron core derived from this star are 7.23 and 1.70 $`\mathrm{M}_{}`$, respectively. Similarly to MW’s simulations, our simulations start after the entire iron core is assumed to have collapsed to form a BH (with $`R_S=4.957\times 10^5`$ cm), but before the helium envelope has collapsed. The model predicts an inner radius of the helium envelope, $`R_{\mathrm{He}}`$, of 2.11$`\times 10^8`$ cm. Outside this radius, we adopt the radial velocity from the stellar model, set $`v_\theta =0`$ and assume a non-zero $`l`$. The angular momentum distribution is chosen such that the ratio between the centrifugal force and the component of gravity perpendicular to the rotational axis is 0.02 at all angles and radii, except where this prescription results in $`l>l_{max}=10^{17}\mathrm{cm}^2\mathrm{s}^1`$; then we set $`l=l_{max}`$. Inside $`R_{\mathrm{He}}`$, we set $`v_\theta =v_\varphi =0`$ and assume a free-fall radial velocity. We compute the initial density inside the helium envelope using the mass continuity equation and the mass accretion rate from the stellar model at $`R_{\mathrm{He}}`$, where $`\rho =1.16\times 10^7\mathrm{g}\mathrm{cm}^3`$ and $`v_r=8.81\times 10^7\mathrm{cm}\mathrm{s}^1`$.

## 3 Results

With these assumptions, there is only one free parameter which defines the strength of the initial magnetic field. In this letter, we present results from a single model for which $`\beta =10^6`$ at the outer boundary.

We find that after a transient episode of infall, lasting 0.13 s, the gas with $`l>2R_Sc`$ piles up outside the black hole and forms a thick torus bounded by a centrifugal barrier near the rotation axis. Soon after the torus forms (i.e., a couple of orbits at $`r=r_i`$), the magnetic field is amplified by both MRI and shear. We have verified that most of the inner torus is unstable to MRI, and that our simulations have enough resolution to resolve, albeit marginally, the fastest growing MRI mode. The torus starts evolving rapidly and accretes onto the black hole. Another important effect of magnetic fields is that the torus produces a magnetized corona and an outflow. The presence of the corona and outflow is essential to the evolution of the inner flow at all times and the entire flow close to the rotational axis during the latter phase of the evolution. We find that the outflow very quickly becomes sufficiently strong to overcome supersonically infalling gas (the radial Mach number in the polar funnel near the inner radius is $`5`$) and makes its way outward, reaching the outer boundary at $`t=0.25`$ s. Due to limited computing time, our simulations were stopped at $`t=0.28215`$ s, which corresponds to 6705 orbits of the flow near the inner boundary. We expect the accretion to continue much longer, roughly the collapse timescale of the Helium core ($`10`$ s), as in MW.

Figure 1 shows the time evolution of the mass accretion rate through the inner boundary (left panel), total magnetic energy (second left panel), neutrino luminosity (third left panel) and radial Poynting and kinetic flux along the polar axis at $`r=190R_S`$ (right panel). Unless otherwise stated, all quantities in this paper are in cgs units. Initially, during a precollapse phase, $`\dot{M}_a`$ stays nearly constant at the level of $`5\times 10^{32}\mathrm{g}\mathrm{s}^1`$. During this phase the zero-$`l`$ gas inside the initial helium envelope is accreted. Around $`t=0.13`$ s, $`\dot{M}_a`$ rises sharply as the gas from the initial helium envelope reaches the inner boundary. However, this gas has non-zero $`l`$ and a rotational supported torus and its corona and outflow form, causing a drop in $`\dot{M}_a`$ after it reaches a maximum of $`2\times 10^{33}\mathrm{g}\mathrm{s}^1`$ at $`t=0.145`$ s. The accretion rate reaches a minimum of $`6\times 10^{31}\mathrm{g}\mathrm{s}^1`$ at $`t0.182`$ s and then fluctuates with a clear long-term increase. This increase is caused by the contribution from gas with $`l<2R_Sc`$, which is directly accreted (without need to transport $`l`$) from outside the main body of the torus. The total mass and angular momentum accreted onto the BH during our simulation (0.3 s) are $`0.1\mathrm{M}_{}`$ and $`3\times 10^{39}\mathrm{g}\mathrm{cm}^2\mathrm{s}^1`$, respectively.

The time evolution of the total magnetic energy (integrated over the entire computational domain) is characteristic of weakly magnetized rotating accretion flows. Most of the magnetic energy is due to the toroidal component of the magnetic field. We note a huge increase of the toroidal magnetic field coinciding with the formation and development of the torus. Both $`B_\varphi `$ and $`B_\theta `$ are practically zero during the precollapse phase of the evolution. But at $`t=0.14`$ s the total energy in $`B_\varphi `$ equals that in $`B_r`$ and just 0.025 s later the $`B_\varphi `$ energy is higher than the $`B_r`$ energy by a factor of 50. At the end of simulations the total kinetic energy from the radial, latitudinal and rotational motion are $`4\times 10^{50}`$, $`6.5\times 10^{49}`$, and $`2.3\times 10^{51}`$ erg, respectively. These gross properties indicate that the magnetic energy is large enough to play an important role in the flow dynamics.

The time evolution of the neutrino luminosity, $`L_\nu `$, shows that the neutrino emission stays at a relatively constant level of $`3\times 10^{52}\mathrm{erg}\mathrm{s}^1`$ after the torus forms. We compute $`L_\nu `$ under the assumption that all the gas in the model is optically thin to neutrinos, so that $`L_\nu `$ is volume integrated $``$ over the entire computational domain. We note that $`L_\nu `$ is dominated by the neutrino emission due to pair capture on free nucleons (the so-called URCA cooling).

The last panel in Fig. 1 shows the area-integrated radial fluxes of magnetic and kinetic energy at $`r=190R_S`$ inside the polar outflow. The outflow is Poynting flux-dominated, with the Poynting flux exceeding the kinetic energy flux by up to an order of magnitude.

Our analysis of the inner flow shows that the outflow is magnetically driven from the torus. Soon after the torus forms, the magnetic field very quickly deviates from its initial radial configuration due to MRI and shear. This leads to fast growth of the toroidal magnetic field as field lines wind up due to the differential rotation. As a result the toroidal field dominates over the poloidal field and the gradient of the former drives an outflow. Figure 2 shows the flow pattern of the inner part of the flow at $`t=0.2735`$ s. The left and right panel show density and $`|B_\varphi |`$ maps, respectively. The maps are overlaid by the direction of the poloidal velocity. The polar regions of low density and high $`B_\varphi `$ coincide with the region of an outflow. We note also that during the latter phase of the evolution not all of the material in the outflow originated in the innermost part of the torus – a significant part of the outflow is “peeled off” the infalling gas at large radii by the magnetic pressure. Figure 2 illustrates that the inner torus and its corona and outflow cannot always prevent the low-$`l`$ gas from reaching the BH. Even the magnetic field cannot do it if the density of the incoming gas is too high. Therefore, we find that the outflow of the magnetic energy (mostly toroidal) from the innermost part of the flow does not always correspond to an outflow of gas (in other words, the Poynting flux and kinetic energy flux can be in opposite directions).

Figure 3 shows the radial profiles of several quantities in our run, angle-averaged over a small wedge near the equator (between $`\theta =86^{}`$ and $`94^{}`$), and time-averaged over 50 data files covering a period at the end of the simulations (from 0.2629 s through 0.2818 s). We indicate the location of the last stable circular orbit by the vertical dotted line in each panel.

The profiles of each variable are not simple power-laws but are rather complex. In particular, the density has a prominent maximum of $`4\times 10^{11}\mathrm{g}\mathrm{cm}^3`$ at $`r=5R_S`$. The gas plus radiation pressure is higher than the magnetic pressure. The rotational velocity is nearly Keplerian inside $`r=6R_S`$ and sub-Keplerian outside this radius.

We measure the Reynolds stress, $`\alpha _{gas}<\rho v_r\delta v_\varphi >/P`$, and the Maxwell stress normalized to the magnetic pressure, $`\alpha _{mag}<2B_rB_\varphi /B^2>`$. Note that Figure 3 shows only the magnitude, not the sign, of the normalized stresses. We find that except for $`r<2.5R_S`$ and $`10R_S<r<12R_S`$, the Maxwell stress dominates over the Reynolds stress in the inner flow. The last panel in Figure 3 shows that the toroidal component of the magnetic field is dominant for $`r<50R_S`$.

We have compared the cooling time scale and the advection time scale in the flow and found that overall the flow is advection-dominated expect for a small region inside the torus where the density reaches its maximum.

## 4 Discussion and Conclusions

We have performed time-dependent two-dimensional MHD simulations of the collapsar model. Our simulations show that: 1) soon after the rotationally supported torus forms, the magnetic field very quickly starts deviating from purely radial due to MRI and shear. This leads to fast growth of the toroidal magnetic field as field lines wind up due to the torus rotation; 2) The toroidal field dominates over the poloidal field and the gradient of the former drives a polar outflow against supersonically accreting gas through the polar funnel; 3) The polar outflow is Poynting flux-dominated; 4) The polar outflow reaches the outer boundary of the computational domain ($`5\times 10^8`$ cm) with an expansion velocity of 0.2 c; 6) The polar outflow is in a form of a relatively narrow jet (when the jet breaks through the outer boundary its half opening angle is $`5^{}`$); 7) Most of the energy released during the accretion is in neutrinos, $`L_\nu =2\times 10^{52}\mathrm{erg}\mathrm{s}^1`$. Therefore it is likely that neutrino driving can increase the outflow energy (e.g., Fryer & Mészáros 2003 and references therein).

Our simulations explore a relatively conservative case where we allow for neutrino emission but do not allow for the emitted neutrinos to interact with the gas or annihilate. The only sources of nonadiabatic heating in our simulations are the artificial viscosity and resistivity.

Our main conclusion is that, within the collapsar model, MHD effects are able to launch, accelerate and sustain a strong polar outflow. We believe that this conclusion will turn out to be largely independent of the initial magnetic field strength in the stellar core, because MRI can rapidly amplify weak fields until they are strong enough to drive a powerful outflow. Since our simulations are non-relativistic, and cover only the innermost region of the collapsing star, we cannot determine whether our outflows are sufficient to yield a GRB. Additional driving could also be necessary. We also find that the outflow is Poynting flux-dominated, and note that this provides favorable initial conditions for the subsequent production of a baryon-poor fireball \[e.g., Fuller, Pruet & Abazajian (2000); Beloborodov (2003); Vlahakis & K$`\ddot{\mathrm{o}}`$nigl (2003); Mészáros (2002)\], or a magnetically dominated “cold fireball” \[Lyutikov & Blandford (2002)\].

ACKNOWLEDGMENTS: DP acknowledges support from NASA under LTSA grants NAG5-11736 and NAG5-12867. MCB acknowledges support from NSF grants AST-9876887 and AST-0307502.

## REFERENCES

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Beloborodov, A.M. 2003, ApJ, 588, 931

Blinnikov, S.I., Dunina-Barkovskaya, N.V., & Nadyozhin, D. K. 1996, ApJS, 106, 171

Fuller, G. M., Pruet, J., & Abazajian, K. 2000, Phys. Rev. Lett., 85, 2673

Fryxell, B. A., M$`\ddot{\mathrm{u}}`$ller, E., & Arnett, W. D. 1989, MPA Rep. 449 (Garching: MPA)

Fryer, C. L., & Mészáros, P. 2003, ApJ, 588, L25

Itoh, N., Adachi, T., Nakagawa, M., Kohyama, Y., & Munakata, H. 1989, ApJ, 339, 354

Itoh, N., Adachi, T., Nakagawa, M., Kohyama, Y., & Munakata, H. 1990, ApJ, 360, 741

Lyutikov, M. & Blandford, R. 2002, APS, APR, 6008

MacFadyen, A., & Woosley, S.E. 1999, ApJ, 524, 262 (MW)

MacFadyen, A., Woosley, S.E., & Heger A. 2001, ApJ, 550, 410

Mészáros, P. 2002, ARA&A, 40, 137

Paczyński, B. 1998, ApJ, 494, L45

Paczyński, B., & Wiita, P. J. 1980, A&A, 88, 23 (PW)

Popham, R., Woosley, S.E., & Fryer, C. 1999, ApJ, 518, 356

Proga, D., & Begelman, M.C. 2003 ApJ, 592, 767 (PB03)

Shakura, N.I., & Sunyaev, R.A. 1973 A&A, 24, 337

Stone, J.M., & Norman, M.L. 1992a, ApJS, 80, 753

Stone, J.M., & Norman, M.L. 1992b, ApJS, 80, 791

Stone, J. M., & Pringle, J. E. 2001, MNRAS, 322, 461

Stone, J.M., Pringle, J.E., & Begelman, M.C. 1999, MNRAS, 310, 1002

Vlahakis, N., & K$`\ddot{\mathrm{o}}`$nigl, A. 2003, ApJ, submitted

Woosley, S.E. 1993, ApJ, 405, 273

Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181

Figure Captions

Figure 1 – The time evolution of the mass accretion rate (left panel), total magnetic energy due to each of the three field components (second left panel), neutrino luminosity (third left panel) and area-integrated radial Poynting and kinetic flux in the polar outflow at $`r=190R_S`$ (right panel). Formally, we define the polar outflow as the region where $`v_r>0`$ and $`\beta <1`$. Note the difference in the time range in the panel with the radial fluxes.

The last panel in Fig. 1 shows the area-integrated radial fluxes of magnetic and kinetic energy at $`r=190R_S`$ inside the polar outflow. Formally, we define the polar outflow as the region where $`v_r>0`$ and $`\beta <1`$. The outflow is Poynting flux-dominated, with the Poynting flux exceeding the kinetic energy flux by up to an order of magnitude.

Figure 2 – Color maps of logarithmic density and toroidal magnetic field overplotted with the direction of the poloidal velocity at $`t=0.2735`$ s. The length scale is in units of the BH radius (i.e., $`r^{}=r/R_S`$ and $`z^{}=z/R_S`$).

Figure 3 – Radial profiles of various quantities from our run, time-averaged from $`0.2629`$ through $`0.2818`$ s. To construct each plot, we averaged the profiles over angle between $`\theta =86^{}`$ and $`94^{}`$. The top left panel plots the density (solid line) and temperature (dashed line). The top middle panel plots the gas pressure (solid line) and magnetic pressure. The top right panel plots the rotational, radial, Keplerian, and Alfv$`\stackrel{´}{\mathrm{e}}`$n velocities (solid, dashed, dot-dashed, and dotted line, respectively), as well as the sound speed (triple-dot dashed line). The bottom left panel plots the angular velocity in units of $`2c/R_s`$. The bottom middle panel plots the Maxwell stress, $`\alpha _{mag}`$, and the Reynolds stress, $`\alpha _{gas}`$ (solid and dashed line, respectively). We calculate the Reynolds stress using eq. (15) in PB03 and show only its amplitude. The bottom right panel plots the radial, latitudinal and toroidal components of the magnetic field (dot-dashed, dashed, and solid line, respectively). The length scale is in units of the BH radius (i.e., $`r^{}=r/R_S`$).


text:
# Effelsberg HI Survey of Compact High-Velocity Clouds

## 1. Introduction

High-velocity clouds (HVCs) were discovered by Muller, Oort, & Raimond (1963) in the 21-cm line emission of neutral atomic hydrogen (HI). HVCs are gas clouds which are characterized by high radial velocities incompatible with a participation in Galactic rotation. Braun & Burton (1999) introduced a subclass of isolated HVCs with small angular sizes of less than $`2^{}`$ FWHM. They compiled a catalog of 66 of these so-called compact high-velocity clouds (CHVCs) from the Leiden/Dwingeloo Survey of Galactic neutral hydrogen (Hartmann & Burton 1997) and proposed that their statistical properties were consistent with a distribution throughout the entire Local Group with distances of the order of $`1\mathrm{Mpc}`$. The typical sizes of CHVCs would then be of the order of $`15\mathrm{kpc}`$ with typical HI masses of a few times $`10^7\mathrm{M}_{}`$.

We mapped 11 CHVCs from the Braun & Burton (1999) and de Heij, Braun, & Burton (2002) catalogs in 21-cm line emission with the 100-m telescope in Effelsberg to validate the results of Braun & Burton (1999). The high angular resolution of $`9^{}`$ HPBW allows us to extract the HI structure of these 11 CHVCs in much more detail than the previous Leiden/Dwingeloo Survey data. In section 2 we describe the sample selection and data acquisition, section 3 outlines our major results, and section 4 includes the summary and conclusions.

## 2. Sample selection and data acquisition

The 11 studied CHVCs were selected on the basis of earlier, less sensitive Effelsberg observations by two criteria. They still had to be classified as CHVCs in the improved catalog by de Heij et al. (2002), and the HI peak column density ratio between the Effelsberg data and the Leiden/Dwingeloo Survey had to be $`N_{\mathrm{Eff}}/N_{\mathrm{LDS}}3`$ so that only the most compact CHVCs were selected for reinvestigation.

In general, each object was mapped with $`11\times 11`$ spectra on a $`9^{}`$ grid (beam-by-beam sampling). We reach a sensitivity of $`\sigma _{\mathrm{rms}}50\mathrm{mK}`$ at $`2.6\mathrm{km}\mathrm{s}^1`$ velocity resolution. These maps allow us to investigate the morphology of our 11 CHVCs and the overall distribution of radial velocities and linewidths of the HI gas. In addition, we observed spectra along an appropriate axis of each CHVC with a longer integration time and a smaller separation of $`4.5^{}`$ or $`6.4^{}`$ between the spectra. Along these deep profiles we reach a sensitivity of $`\sigma _{\mathrm{rms}}30\mathrm{mK}`$, allowing us to obtain the column density profile in more detail.

## 3. Results

Table 1 summarizes the extracted physical parameters of the observed CHVCs, derived by fitting a Gaussian to the spectral lines. Radial velocities are mainly negative which is simply a selection effect. The mean linewidth of our 11 CHVCs is $`\mathrm{\Delta }v=23\pm 7\mathrm{km}\mathrm{s}^1`$, the mean peak brightness temperature is $`T_\mathrm{B}=1.6\pm 0.8\mathrm{K}`$, and the average peak column density is $`N_{\mathrm{HI}}=(5.3\pm 1.6)10^{19}\mathrm{cm}^2`$. In some cases, we find evidence for two distinct gas phases in the spectral lines, indicating the existence of a core of cold neutral gas embedded in a diffuse envelope of warm neutral gas.

Among the 11 investigated objects, only CHVC 148$``$82 appears to be spherically-symmetric (fig. 1 (a)). All other clouds have a more or less complex morphology. 4 CHVCs reveal a head-tail structure which indicates that their envelopes of diffuse, warm gas are stripped off by ram-pressure interaction with an ambient medium. The conception of an interaction process is also supported by the distribution of radial velocities and linewidths across all 4 clouds. Fig. 1 (b) shows CHVC 040$`+`$01 as an example for a pronounced head-tail structure.

Another example for a head-tail CHVC is given in fig. 2 (a). CHVC 017$``$25 is a relatively compact object, and its slightly asymmetric shape already indicates that gas might have been stripped off the cloud, creating a diffuse, faint tail in north-western direction. Furthermore, CHVC 017$``$25 shows a clear two-component structure in the HI line profiles. Spectral lines seem to disclose a superposition of a narrow Gaussian component of cold gas and a broad component of warm gas. The cold component shows line widths of about $`7\mathrm{km}\mathrm{s}^1`$ FWHM, resulting in an upper limit for the gas temperature of roughly $`1000\mathrm{K}`$. The warm component discloses much larger line widths around $`20\mathrm{km}\mathrm{s}^1`$ FWHM, indicating an upper limit for the gas temperature of about $`9000\mathrm{K}`$. Both gas components can be investigated separately by a Gaussian decomposition of the spectral lines. Fig. 2 (b) shows the distribution of HI column densities of the cold (open circles) and warm (filled circles) gas along the profile indicated by the dashed arrow in fig. 2 (a). One can clearly distinguish a compact cold core surrounded by an extended, diffuse envelope of warm gas. Furthermore, the cold core and the warm envelope are spatially separated from each other. At the south-eastern edge of the cloud, only the cold gas can be traced while the warm gas forms an extended, faint tail in north-western direction. These results suggest that the diffuse envelope of warm gas is currently being stripped off the compact cold core of CHVC 017$``$25 by ram-pressure interaction with an ambient medium.

Two CHVCs from our sample show up with a bow-shock shaped structure which again indicates the existence of an ambient medium. Fig. 1 (c) shows CHVC 172$``$60 as an example for a bow-shock shaped object. In both CHVCs, radial velocities indicate that the gas is decelerated at the presumable front of the cloud. At the same time, linewidths increase, indicating higher turbulence or heating of the gas. Again, these observations can best be explained by ram-pressure interaction with a medium surrounding the CHVCs.

The remaining 4 CHVCs from our sample reveal an irregular shape. Some of these irregular clouds have large radial velocity gradients, suggesting a rotation. In all cases, the irregular morphology again indicates a disturbance of the objects.

## 4. Summary and conclusions

We showed that 10 of our 11 investigated CHVCs reveal a complex structure and that many of them show signs for ram-pressure interaction with a surrounding medium. This is in opposite to the original idea of CHVCs being the gaseous counterparts of primordial Dark-Matter halos spread across the entire Local Group. Instead, the observed ram-pressure effects and different distance estimates (Westmeier 2003) indicate that CHVCs constitute a circumgalactic population with typical distances of the order of $`100\mathrm{kpc}`$. At such distances, the observed interactions could be caused by an extended Galactic halo gas. Typical HI masses of CHVCs would then be of the order of a few times $`10^5\mathrm{M}_{}`$ with typical sizes of about $`1\mathrm{kpc}`$.

### Acknowledgments.

This work is based on observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg.

## References

Braun, R. & Burton, W. B. 1999, A&A, 341, 437

Hartmann, D. & Burton, W. B. 1997, Atlas of Galactic Neutral Hydrogen, Cambridge University Press

de Heij, V., Braun, R., & Burton, W. B. 2002, A&A, 391, 159

Muller, C. A., Oort, J. H., & Raimond, E. 1963, C.R.Acad.Sci.Paris, 257, 1661

Westmeier, T. 2003, Diploma Thesis, Universität Bonn


text:
# A multi-particle model of the 3C 48 host

## 1 Introduction

3C 48 (Barkhouse & Hall 2001) was the first QSO to be discovered optically (Matthews et al. 1961; Matthews & Sandage 1963) and the first QSO to be directly identified with a host galaxy (Boroson & Oke 1982). Its basic properties are listed in Table 1. 3C 48 has attracted much attention regarding the proposed evolutionary sequence of active nuclei (Sanders et al. 1988). According to this scheme, interactions and mergers of galaxies trigger an evolution via ultra-luminous infrared galaxies (ULIRGs) to QSOs. The observational evidence, however, is hampered by the fact that many transitionary objects show only dubious indications of past or recent mergers. Clarification requires detailed multi-particle modeling which helps with disentangling the complex spatial structure of merger remnants.

3C 48 is an example of a transitionary object with prototypal properties in many respects: It has the typical far-infrared excess, originating from thermal radiation of dust which is heated by the quasar nucleus and by newly forming stars in the host galaxy (Neugebauer et al. 1985; Stockton & Ridgway 1991). Large amounts of molecular gas (Scoville et al. 1993; Wink et al. 1997) indicate the possibility of a young stellar population in the host. Finally, long-slit spectroscopy gives evidence for an ongoing starburst in the host which currently seems to be close to its maximum activity (Canalizo & Stockton 2000). But the merger scenario for 3C 48 is still unclear: Indeed, the host has a significant tail-like extension to the northwest whose tidal origin is rather compelling with regard to the kinematics and ages of its stars (Canalizo & Stockton 2000). However, the nature of the apparent second nucleus 3C 48A about $`1\mathrm{}`$ northeast of the QSO (Stockton & Ridgway 1991) and the location of the expected counter tidal tail remain an unsolved problem. 3C 48A could as well be due to the radio jet (Wilkinson et al. 1991) interacting with the dense interstellar medium (Chatzichristou et al. 1999). A feature at the southeast of the 3C 48 host, previously interpreted as a counter tidal tail (Boyce et al. 1999), has turned out to be a background galaxy (Canalizo & Stockton 2000). Instead, a counter tail extending from the southeast to the southwest is suspected (Canalizo & Stockton 2000) but not yet identified. Canalizo & Stockton (2000) suggest that such a location of the two tidal tails might be explicable by a certain projection of the merger scenario used to simulate the ”Antennae”.

This paper reports the first successful multi-particle model for the 3C 48 host. Suggesting simple solutions for the 3C 48A problem and the counter tail problem, the model largely resolves doubts about the merger hypothesis for 3C 48.

## 2 Methods

### 2.1 Numerical simulations

The stellar-dynamical 3-dimensional N-body simulations are performed with TREESPH (Hernquist & Katz 1989), a tree code used in its non-collisional mode. Stable initial particle distributions for the model galaxies are set up with BUILDGAL (see Hernquist 1993). The spatial orientation of the two disk galaxies is parameterized by their inclinations $`i`$ with respect to the orbital plane and their pericentric arguments $`\omega `$ as introduced by Toomre & Toomre (1972). Both galaxies have extended mass distributions so that the orbit of their encounter is not Keplerian but decaying. Two descriptions can be used to characterize these orbits: The first one is a pseudo-Keplerian description using the parameters of eccentricity $`e`$, pericentric distance $`r_{\mathrm{peri}}`$, and angle to pericentre $`\mathrm{\Omega }_{\mathrm{peri}}`$ of the corresponding Keplerian orbit for which the total mass of each galaxy is associated with a point mass at the respective centre of mass. The second one is a direct description of the decaying orbit using the true apocentric and pericentric distances $`r_{\mathrm{apo}}`$ and $`r_{\mathrm{peri}}`$ of the first passage to define an eccentricity in its generalized formulation $`e=(r_{\mathrm{apo}}r_{\mathrm{peri}})/(r_{\mathrm{apo}}+r_{\mathrm{peri}})`$. For convenience, the system of units, which remains intrinsically scale-free, is scaled to the system suitable for 3C 48.

The results of the simulations are analyzed as 2-dimensional projections. In order to mock the pixel array data of imaging observations, the particles are sorted into a $`512\times 512`$ grid. The virtual pixel values are computed by adding up all particles located in a grid cell along the line-of-sight. Without any special weighting of a nuclear component, the mock images are comparable to QSO-subtracted images of the 3C 48 host. Spectra for each grid cell are generated by sorting the particles into velocity channels according to their respective line-of-sight velocities. Thus, an average stellar line-of-sight velocity is assigned to each virtual pixel. The resulting data arrays are spatially smoothed by Gaussian convolution and converted into FITS format to facilitate the further data processing with standard astronomical software.

### 2.2 Observational data on 3C 48

The data presented by Canalizo & Stockton (2000) are used for comparing the simulations with observations. They provide information about the optical surface-brightness of the QSO-subtracted 3C 48 host (Fig. 1 therein) and about the stellar kinematics along the four slits A, B, C, G (Fig. 1 and Table 2 therein). In reference to these data, the basic proportions of the main body of the 3C 48 host are classified by dimensionless length ratios (left panel of Fig. 1 and left column of Table 3). Such a comparison is independent of the length scaling, in contrast to the comparison of line-of-sight velocities which requires a positioning of the four slits on the mock image.

Having determined the final physical length unit of the simulations, the curvature of the northwestern tidal tail is compared by using an angle-versus-distance plot (right panel of Fig. 1 and Fig. 5).

## 3 The 3C 48 look-alike

Different mass ratios of the initial galaxies, different snapshots during the merger process, and different projection angles of the merger remnants were probed in a still limited parameter study.

The nearest 3C 48 look-alike is found for the merger of two identical galaxies whose physical and numerical properties are given in Table 2. With these parameters the galaxies are similar to spirals of type Sb.

The experimental setup is the same as used for simulations of the ”Antennae” – i.e. both galaxies are symmetrically oriented with $`i_1=i_2=60\mathrm{°}`$ and $`\omega _1=\omega _2=30\mathrm{°}`$ (see Toomre & Toomre 1972; Barnes 1988). The model galaxies are initialized near the apocentre of the corresponding elliptical Keplerian orbit which is defined by the eccentricity of $`e=0.5`$, the pericentric distance of $`r_{\mathrm{peri}}=20`$ kpc, and the period of 1.2 Gyr. The time step in the simulations is fixed to 1 Myr which guarantees that energy is conserved to better than 1% during the merger. The true decaying orbit is characterized by the generalized eccentricity of $`e=0.8`$ and the pericentric distance of $`r_{\mathrm{peri}}=7`$ kpc. The outer regions of the two galaxy bulges begin to merge after $`272.5`$ Myr, just before pericentric passage. The merging is not a single process but the centres of the two bulges are repeatedly flung apart before they settle down in one common density peak.

The nearest 3C 48 look-alike emerges after 461.1 Myr. Two projections of this merger remnant are shown in Fig. 2. In the left panel (”Antennae” look-alike), the view is perpendicular to the orbital plane x-y. In the right panel, (3C 48 look-alike) the orbital plane is tilted southwards, westwards, and counterclockwise by $`120\mathrm{°}`$, $`160\mathrm{°}`$, and $`116\mathrm{°}`$, respectively. The proportions of the 3C 48 look-alike are listed in the right column of Table 3. The positions of the four slits A, B, C, G on the look-alike host and the resulting physical coordinate system are shown in Fig. 3.

In Fig. 4, the scaled velocities along the slits are compared to the confidence region of stellar line-of-sight velocities given for 3C 48 by Canalizo & Stockton (2000). The angle-versus-distance comparison for the curvature of the northwestern tidal tail is shown in Fig. 5.

## 4 Discussion

### 4.1 The nature of 3C 48A

Optical and near-infrared images of the 3C 48 host show two luminosity peaks at the positions of 3C 48 and 3C 48A, the latter being located about $`1\mathrm{}`$ northeast of 3C 48 (e.g. Stockton & Ridgway 1991; Canalizo & Stockton 2000). With softening lengths of $`0.25`$ kpc ($`0\stackrel{}{.}06`$) for the bulge and disk particles (see Table 2), the spatial resolution of the simulations is high enough to identify corresponding density peaks in the 3C 48 model. As shown in the small inset in Fig. 3, the centres of the bulge components of the two merging galaxies are still separated at the stage of the 3C 48 look-alike. Their distance of about $`0\stackrel{}{.}6`$ (2.5 kpc) and their relative positions on a southwestern to northeastern axis are similar to the observed configuration of 3C 48 and 3C 48A. Thus, a scenario with 3C 48 and 3C 48A being the two centres of merging galaxies is possible. However, the exact configuration of the density peaks in the simulations is very sensitive to the projection angle and to the time at which the snapshot is taken. About 20 Myr later, the two peaks have already merged into one. Furthermore, the purely stellar-dynamical model does not address the question of a possible nuclear activity at the positions of 3C 48 and 3C 48A. A detailed discussion of nuclear activity depends on whether or not a black hole exists at the mentioned positions and on the respective fueling rates.

### 4.2 Location of the counter tidal tail

Since each of the two merging disk galaxies forms a tidal tail, the missing second tidal tail has always been a caveat of the merger hypothesis for 3C 48. Here, the simulations suggest a simple solution: At the projection angle of the 3C 48 look-alike, the second tidal tail is mainly located in front of the body of the host and, therefore, severely foreshortened. It extends from the southwest towards the northeast, roughly along slit B, so that measurements along this slit trace a mixture of line-of-sight velocities from the body and from the tail. This could explain why the observed and the simulated line-of-sight velocity signatures along slit B (Fig. 4) are dominated by scattering around a mean velocity close to zero. Slits A, C, and G, in contrast, are characterized by large absolute line-of-sight velocities (up to $`200\mathrm{to}300\mathrm{km}\mathrm{s}^1`$) and strong variations along the slits. A counter tidal tail extending from the southwest towards the northeast in front of the main body of the 3C 48 host is a completely new alternative. Regarding the information about stellar kinematics, this location seems to be more likely than the two suggested tails in the southeast and from the southeast towards the southwest (Boyce et al. 1999; Canalizo & Stockton 2000) which have failed identification so far.

### 4.3 The evolutionary history of 3C 48

Conclusions about the orbital parameters for 3C 48 and the original parameters of the merging galaxies can only be tentative. The orbital period of the best fit model amounts to about 20% of the age of the universe at the redshift of 3C 48 ($`5.4`$ Gyr). A merger scenario with such an orbital period is plausible, assuming an initially highly eccentric orbit of the merging galaxies which is transformed into a bound orbit by dynamical friction of their dark matter halos (e.g. Jones & Stein 1989). It has been found that the morphology and the kinematics of tidal tails are very sensitive to the rotation curve of the interacting model galaxies (Dubinski et al. 1996; Mihos et al. 1998; Dubinski et al. 1999). Thus, instead of two identical galaxies, an alternative model for 3C 48 could start from two galaxies with different rotation curves so that only one of them forms an extended tidal tail. However, even in its generality the multi-particle model presented in this paper gives rather compelling evidence that the formation of 3C 48 is linked to a merger process. Therewith, 3C 48 ranks among these transitional objects which support the evolutionary scenario (Sanders et al. 1988) in its original merger-driven definition.

###### Acknowledgements.

Our special thanks go to Prof. Dr Lars Hernquist who kindly provided the codes TREESPH and BUILDGAL and gave helpful advice. This project was supported in part by the Deutsche Forschungsgemeinschaft (DFG) via grant SFB 494. J. Scharwächter is supported by a scholarship for doctoral students of the Studienstiftung des deutschen Volkes.


text:
# Chemical Evolution in Hierarchical Clustering Scenarios

## 1. Introduction

The current cosmological paradigm predicts a hierachical building up of the structure. This scenario has largely supported by observations in the local and high redshift Universe. However, few but important descrinpances remain to be explanied such as the inner distribution of dark matter in small haloes, the origin of the high baryonic content and metallicity of the intercluster medium, etc (Shanks 2004). Most of these open problems requires the comprehensive understanding of galaxy formation, including all the involved scales. Numerical simulations have played an important role in understanding and proving the connection between theory and observations at different scales. However, the modelling of small scale physical processes such as gas cooling, star formation and chemical evolution within a cosmological framework presents different difficulties because of our poor understanding of the relevant physics and numerical limitations. In particular, only recently chemical evolution has started to be treated within cosmological hydrodynamical simulations (Yepes et al. 1998; Cen & Ostriker 1999; Mosconi et al. 2001).

The dramatic building up of observational evidences on chemical properties of the stellar population and interstellar medium in the Universe has proved to have a significant impact on the study of galaxy formation. Chemical elements are produced in the stellar interior and ejected in a pacific or violent way after characteristic timescales determined by stellar evolution. Afterwards, these elements can be mixed up by mergers and interactions. It is also important to notice that the metal production is directly linked to the star formation history which can be also affected by the environment (Barton et al. 2000; Lambas et al. 2003; Balogh et al. 2004). As a consequence, it is expected that chemical patterns store information of the process of galaxy formation ( Freeman & Bland-Hawthorn 2002).

Another important motivation to include a chemical treatment in galaxy formation models is related to the fact that baryons cool according to their metallicity. Sutherland & Dopita (1993) showed that a gas cloud with solar metallicity cools at a rate of up to few orders of magnitude more efficiently than a primordial gas cloud. Hence, the enrichment of the interstellar medium can have dynamical consequences and affect the formation of the structure.

In this work, we present first results of detail treatement of chemical evolution in Gadget2 (Springel & Hernquist 2002) based on the implementation of Mosconi et al. (2001). In Section 2 we briefly summarize the main important aspects of the chemical code. Section 3 shows the first results and Section 4 presents the conclusions.

## 2. The model

The chemical model has been implemented within Gadget2 (Springel & Hernquist 2002) which is based on a modified version of SPH that conserves entropy (energy and angular momentum). Stars are formed according to a probability function defined by the relation between cooling and dynamical time (see also Lia et al. 2001). A standard Salpeter Initial mass function has been assumed. Supernovae II are considered to form from stars with $`M>8M_{}`$ and to produce metals according to Woosley & Weaver (1995) yields. Supernova Ia are assumed to originate from binary systems following the model of Thielemann et al. (1993). We have assumed a random sampling of the lifetime of binary systems in the range 0.1-5 Gyrs.

Chemical elements are distributed within gas neighbouring particles by using a spline kernel function. In order to estimate the correct environment of stars when SN explode we re-calculate the gaseous neighbours of the stellar particles each time elements are released. In this implementation we have not included the treatment of energy feedback.

We run cosmological tests of $`\mathrm{\Lambda }`$-CDM cosmologies consistent with the concordance model ($`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_m`$=0.3, $`H_o=70\mathrm{kms}^1\mathrm{Mpc}^1`$). The simulations correspond to a volume box of $`10h^1`$ Mpc ($`h=0.7`$) with $`2\times 50^3`$ particles. Three different runs were performed using different cooling functions: solar, primordial and metal-dependent.

## 3. Results

The cosmological simulations are analyzed from a global point of view and then, by identifying the gravitational bounded structures. In Figure 1 (left panel) we show the comoving star formation $`\rho _{\mathrm{SFR}}`$ and metal production $`\rho _\mathrm{Z}`$ rates for the three runs with solar (solid lines), primordial (dashed lines) and metal-dependent (dotted-dashed lines) cooling functions. We can see that the largest difference is found for the run with fixed solar cooling and for $`z3`$. The $`\rho _{\mathrm{SFR}}`$ of the metal-dependent cooling run can be up to $`50\%`$ lower than that of the solar cooling test, and up to $`20\%`$ higher than the primordial cooling run. Similar relations are found for the metal production rate $`\rho _\mathrm{Z}`$ . Figure 1 (right panel) shows the fraction of gas with $`T>10^5`$ K for the three runs, with the largest gas fraction associated to the primordial cooling run, as expected. the difference between primordial and metal-dependent cooling are localized in the densest regions.

At $`z=0`$ we identified galaxy-like objects (GLOs) at their virial radius by combining the friends-of-friends algorithm with a density contrast method in the three simulations (with solar, primordial and metal dependent cooling). We constructed the star formation history of the GLOs and analyzed the chemical properties of their stellar population and interstellar medium. In Figure 2a we show the SFR history of a typical GLO in the three runs. In the small box we have plotted the ratio of the SFR of the solar and metal-dependent cooling to that of the primordial one. We can see that the differences are important for $`z<3`$, but do not present a regular behavior. This suggests that the non-linear evolution of the structure determines complicate star formation history and chemical patterns that can produce non-linear responses to the characteristics of the cooling functions. Our model will allow us to try to understand the link between the dynamical evolution and the star formation history and chemical patterns, by studying the particular history of evolution of GLOs in different environments.

It was shown by Tissera, Lambas & Abadi (1997) that the combination of the history of star formation with models of population synthesis models results in a powerful tool to cnonect theory with observations. The treatment of the chemical enrichment in the simulations improve the technique since each simulated stellar population has its age and metallicity estimated consistently with the history of formation in a hierarchical scenario. With this new chemical code, we can also assess the effects of cooling the gas according to its metallicity. For this purpose, we constructed the synthetic spectral distribution of the GLOs but combining the SEDs of Bruzual & Charlot (2003, BC03) according to the metallicity and age of stars. We also estimated the simulated spectra (SS) of the GLOs by assuming fixed metallicity at $`Z=0.0001Z_{}`$ and $`Z=Z_{}`$ for the SEDs of BC03. In Figure 2b we show the SS for the GLOS of 2a for the three runs. Note that the SS has been calculated with the same SFR history ( 2a, solid line). From this figure, we can see that assuming a low fixed metallicity leads to a very different energy distribution. For this GLO, the real SS is more similar to that produced by assuming fixed solar metallicity.

We also analyzed the simulated spectra for the same GLO in the three runs (solar, primordial and metal dependent cooling), and estimated the SS by using the full information on age and metallicity. Hence, in this case we will looking at the impact of different gas cooling on the star formation and chemical history through the differences in the energy distribution. Figure 3 we show the evolution of the SS for the GLO of Figure 2 in the three runs. As it can be seen, the GLOs in the solar cooling run exhibits the more luminous spectra since because of the high efficient cooling, a large fraction of stars formed quickly before the gas has time to enriched. In the case of the corresponding GLO in the runs with primoridal and metal dependent cooling, the interstellar medium has time to get enriched by previous star generations. Hence, in this case, the SS of the GLO metal-dependent cooling run is more similar to that obtained from the primordial one.

## 4. Conclusions

We have incorporated a chemical evolution algorithm within Gadget2 which allows the description of the enrichment of baryons within a cosmological framework. A first analysis shows that in $`\mathrm{\Lambda }`$-CDM scenarios, the comoving star formation rate could be underestimated by up to $`20\%`$ if the primordial cooling function is adopted. However, within individual galactic systems it is complicated to disentangle the effects of cooling the gas consistently with its metallicity owing to the non-linear growth of the structure which affects the star formation, chemical histories and the mixing of metals in a non-linear way. First results show that the star formation history could be either underestimated or overestimated by an order of magnitude depending on the particular history path. We also show the simulated spectra obtained by combining the age and metallicity information of stars in the simulation with the SEDs of Bruzual & Charlot (2003).Future works will focus on the study of the link between the characteristic of the spectra and history of evolution of the galactic systems.

## Acknowledgments

We thank Simon White and Volker Springel for useful discussions, and V.Springel for making Gadget2 available for this work. We acknowledge the support of Fundanción Antorchas and DAAD that made this work possible. C. Scannapieco thanks the Alexander von Humboldt Foundation, the Federal Ministry of Education and Research, and the Programme for Investment in the Future (ZIP) of the German Government for partial support. This work was partially funded by CONICET.

## References

Balogh et al. 2004, submitted

Barton, E. J., Geller, M. J. & Kenyon, S. J., 2000, ApJ, 530, 660

Bruzual, G. & Charlot, S., 2003, in preparation

Freeman, K., & Bland-Hawthorn, J., 2002, ARA&A, 40, 487

Lambas, D. G., Tissera, P. B., Alonso, M. S., & Coldwell, G., 2003, MNRAS, in press

Lia, C., Portinari, L, & Carrari, G. 2002, MNRAS, 339, 821

Mosconi, M. B., Tissera, P. B., Lambas, D. G., & Cora, S. A. 2001, MNRAS, 325, 34

Springel, V., & Hernquist, L., 2002, MNRAS, 333, 649

Sutherland, R. S., Dopita, M. A., 1993, ApJS, 88, 253

Thielemann, F. K., Nomoto, K., & Hashimoto, M. 1993, Prantzos, N., Vangoni-Flam, E., Cassé N., eds., Origin and Evolution of the Elements, p.299

Tissera, P. B., Lambas, D. G., Abadi, M., 1997, MNRAS, 286, 384

Yepes, G., Elizondo, D., & Ascasibar, Y., 1998, Ap&SS, 263, 31

Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181