text:
# All Sky Doppler Extrasolar Planet Surveys with a Multi-object Dispersed Fixed-delay Interferometer

## 1. Introduction

Since the discovery of 51 Peg B (Mayor & Queloz 1995), the number of known planets has gone up dramatically and today more than 100 such companions are known. Most of these systems have been discovered using high resolution echelle spectrographs. These instruments are routinely achieving a radial velocity (RV) precision of $`\sigma `$ = 5-15 m/s and even as low as 3 m/s in the best cases (Butler et al. 1996; Baranne et al. 1996; Vogt et al. 2000). The echelle spectrographs themselves are large and expensive, with a complicated point spread function which has to be modeled appropriately. The high resolution along with the large wavelength coverage is paid for by the huge costs, the complexity of the instrument and the low throughput of around 1-4% (Vogt et al. 2000; D’Odorico et al. 2002). The observational challenge now lies in detecting as many planetary systems as possible to understand the statistical distribution of planetary masses and distances from the host star. Detecting a large sample of planets will also help in understanding the physical processes underlying planetary formation. A larger sample of stars needs to be surveyed with velocity precision less than 10 m/s to make this possible. Such a task will need large amounts of time on telescopes with current echelles and this is a significant challenge. The next generation of RV survey instruments must have higher efficiencies, should be able to do multi-object RV surveys and should extend the wavelength range to near IR and near UV. The development of our instrument, a multiple object dispersed fixed-delay interferometer, a new generation RV instrument, has been driven in part by these needs.

## 2. Principle and Unique Properties of Fixed-delay Interferometer

The use of a fixed-delay interferometer for Doppler RV measurements is completely different from the current echelle approach. Instead of measuring the absorption line centroid shifts in the echelle approach, the RV is measured through monitoring interference fringe shifts (Ge 2002). The original idea for using a fixed-delay interferometer for high precision Doppler RV measurements was proposed by two groups (Gorskii & Lebedev 1977; Beckers & Brown 1978). This interferometer with a narrow bandpass has been successfully used for very high Doppler precision measurements of the sun ($``$ 3 m/s, Kozhevatov et al. 1995, 1996; sub m/s precision for the GONG measurements, Harvey 2002 private communication). The concept of combining of a fixed-delay interferometer with a moderate resolution spectrometer, or a post-disperser, for broad band operations for high precision stellar Doppler measurements was proposed by Dave Erskine at LLNL in 1997. The initial lab experiments and telescope observing with a prototype demonstrated its feasibility (Erskine & Ge 2000; Ge et al. 2002a). A theory for this new instrument concept was developed by Jian Ge (Ge 2002). In this interferometer approach, the instrument response, determined by the two beam interference, is a simple and well-defined sinusoidal function. For comparison, the echelle response is considerably more complex due to the interference among thousands of divided beams from grating grooves.

In a fixed-delay interferometer (FDI), a fixed optical delay, $`d`$, is applied to one of the beams. Therefore, the interference happens at very high interference order, $`m`$, determined by $`m=\frac{d}{\lambda }`$, where $`\lambda `$ is the operating wavelength. The Doppler RV motion will shift the fringes of stellar absorption lines to neighboring orders. The corresponding Doppler velocity shift is

$$\mathrm{\Delta }\upsilon =\frac{c\lambda }{d}\mathrm{\Delta }m=\frac{c\lambda }{d}\frac{\mathrm{\Delta }\varphi }{2\pi },$$
(1)
where $`\mathrm{\Delta }\varphi `$ is the measured phase shift of a fringe. If the absorption line density is constant over the observed band, then the total observed Doppler error is

$$\sigma _{t,fringe,ob}\frac{1}{\sqrt{N_o}}\frac{1.1c\lambda }{D_ol_c\sqrt{F_o}}\sigma _{t,fringe,I}$$
(2)
where $`N_o`$ is the total number of absorption lines covered by the array, $`D_o`$ is the observed absorption line depth, $`l_c`$= $`\lambda `$2/$`\mathrm{\Delta }\lambda `$I is the coherence length of the interferometer beam for a bandwidth of $`\mathrm{\Delta }\lambda `$I, the intrinsic width of typical stellar absorption lines, $`F_o`$ is the observed flux within each observed fringe (or absorption line), and $`\sigma _{t,fringe,I}`$ is the intrinsic Doppler precision at infinite spectral resolution of a post-disperser (see Ge 2002 for details). This indicates that the Doppler sensitivity of the FDI is independent of the spectral resolving power of the post-disperser, contrary to echelle spectroscopy. In the echelle, the total observed Doppler error is

$$\sigma _{t,echelle,ob}\frac{1}{\sqrt{N_o}}\frac{1.1c\mathrm{\Delta }\lambda _o}{D_o\lambda \sqrt{F_o}}\left(\frac{\mathrm{\Delta }\lambda _o}{\mathrm{\Delta }\lambda _I}\right)\sigma _{t,echelle,I}\left(\frac{\mathrm{\Delta }\lambda _o}{\mathrm{\Delta }\lambda _I}\right)\sigma _{t,fringe,ob},$$
(3)
where the observed line FWHM is $`\mathrm{\Delta }\lambda _o=\sqrt{\mathrm{\Delta }\lambda _I^2+\mathrm{\Delta }\lambda _e^2},`$ and $`\mathrm{\Delta }\lambda `$e is the FWHM of the echelle response, $`\sigma _{t,echelle,I}`$ is the intrinsic Doppler precision at infinite spectral resolution of an echelle. The Doppler error in the echelle approach strongly depends on the echelle resolving power and stellar intrinsic line width (Bouchy et al. 2001; Ge 2002). For a solar type star with absorption lines of FWHM $``$ 5 km/s (Dravins 1987), at moderate resolution (such as $`\mathrm{\Delta }\lambda `$o$``$ 10$`\mathrm{\Delta }\lambda `$I, or $`R`$ 6000), the echelle approach has $``$ 10 times higher Doppler error than the interferometer approach. At high resolution, such as $`R`$ 60,000, being used for planet detection (e.g., Vogt et al. 1994; D’Odorico et al. 2000), the Doppler sensitivity for a solar type star is still $``$ 1.4 times worse than the interferometer.

The independence of Doppler sensitivity from the post-disperser resolving power in the interferometer approach opens up new possibilities for RV studies. The use of low resolution but high efficiency post-dispersers can significantly boost the overall detection efficiency, dramatically reduce the instrument size and cost and allow single dispersion order operations for multiple object observations. Full sky coverage for an RV survey for planets becomes possible with wide field telescopes. Multiple object capability is one of the most significant advantages for this interferometer approach. The simple and stable response function in the interferometer approach leads to potential low systematic errors, which may allow this approach to reach sub m/s Doppler precision.

Another exciting possibility with this interferometer technique is to extend RV surveys to wavelengths other than the visible, previously not covered by echelle surveys. Since the interferometer response is simple and stable, there is no need to calibrate the instrument response in contrast to the echelle; only wavelength calibration is required. Hence, reference sources with a lower line density than the iodine, which is popularly used in the echelle, can be used. Therefore, this instrument can be easily adapted to other wavelengths, in which more photon flux and stellar absorption lines are available for precision Doppler RV measurements. For instance, late M, L and T dwarfs have peak fluxes in the near-IR. Integration time can be significantly reduced if the IR spectra can be monitored. B and A main sequence stars and white dwarfs have very broad intrinsic absorption lines dominated by the Balmer series. The intrinsic Doppler error for each Balmer line is about the same as that for late type stars. Since there are only a few dozen broad lines that can be used for early type stars while $``$ 1000 lines can be used for late types, the overall observed Doppler error is about $``$ 10 times higher for the early types than late types for the same S/N data. Since early type stars are usually much brighter in the near UV than late types at the same astronomical distance (e.g. an A star is about 100 times brighter than a G type star in the visible), it is possible to achieve $``$ a few m/s Doppler precision by increasing S/N by a factor of $``$ 10.

## 3. Performance of a Prototype Dispersed Fixed-delay Interferometer

In 2000-2001, we developed a prototype dispersed fixed-delay interferometer, called Exoplanet Tracker (ET), similar to an earlier version built by Jian Ge and his collaborators at LLNL (Ge et al. 2002a). The instrument consists of a Michelson interferometer with 7 mm fixed optical path difference, a Czerny-Turner type spectrometer with a 100 mm diameter collimator beam and a 1k$`\times `$1k Photonics CCD camera with 24 $`\mu `$m pixel size. The first light stellar observations were conducted at the Hobby-Eberly 9 m telescope (HET) in October 2001 (Ge et al. 2002b). The spectral resolution is $`R`$ = 6700. The wavelength coverage is 140 Å with the CCD. The measured Doppler precision is approaching a photon-noise limit as demonstrated by a direct comparison between predicted theoretical RV errors and measured values shown in Figure 1(left). This figure also indicates that the iodine calibration contributes major errors in the measurements due to its very low fringe visibility compared to that for the star.

ET was further tested at the Palomar 5-m telescope in Dec. 2002. Our results demonstrate that multiple object observing is feasible with the interferometer. Figure 1(right) shows two adjacent fringe spectra obtained simultaneously with ET. One fringe data is from Aldebaran. The other is from a ThAr lamp. The two fringe spectra occupy only $``$ 1/8 CCD detector area. Therefore, we can simultaneously cover $``$ 15 fringe data from 15 stars if multiple stellar beams are available.

A modified version of ET was used for an engineering run at the KPNO 2.1-m telescope in August 2002 before we install a permanent one in 2003 for a long-term survey. The old $`f`$/10 spectrograph was replaced by an $`f`$/7.5 spectrograph. The KPNO 1k$`\times `$3k back-illuminated CCD with 15 $`\mu `$m pixels was used instead of the old 1k$`\times `$1k CCD. The wavelength coverage has been increased to 270 Å due to the faster instrument focal ratio and larger detector array. An $`f`$/8 telescope beam is fed into a 200 $`\mu `$m fiber, which matches a 2.5 arcsec stellar image. Due to the focal ratio degradation, the output focal ratio of the fiber is $`f`$/6, which is converted to $`f`$/7.5 to feed the spectrograph. The spectrograph entrance slit width was dialed to about 180 $`\mu `$m, which causes about 30% photon loss at the slit. The FWHM of each absorption line is sampled by 12 pixels. Under 1.5 arcsec seeing conditions, the total instrument throughput including the sky, telescope transmission, fiber loss, instrument and iodine cell transmission and detector quantum efficiency is 3.4%. This was achieved by feeding only one of the two interferometer outputs into the spectrograph. This allows us to routinely observe stars as faint as V = 7.6 during the whole run. Although the seeing conditions were never better than 1.5 arcsec, we still were able to continually monitor 6 stars (Arcturus, $`\eta `$ Cas, $`\upsilon `$ And, 51 Peg, 31 Aql, and HD 209458) over 8 nights when sky was relatively clear.

Part of the data have been reduced. Figure 2(left) shows a velocity curve for $`\eta `$ Cas within an hour. The measured RV values well match the predicted ones caused by the Earth’s motion. The RMS residual from this measurement indicates a Doppler precision of 2.9 m/s in a short period. Measurements from several other nights indicate we have reached a Doppler sensitivity of $``$ 3-8 m/s for $`\eta `$ Cas. Figure 2(right) shows a RV curve from 51 Peg after Earth’s motion is subtracted, superimposed with a predicted curve from previous echelle measurements. Our result is consistent with previous measurements, demonstrating that ET is capable of detecting extrasolar planets. Large long term RV measurement errors are mainly caused by data reduction and also instrument calibration. A better version of data reduction software package is being developed. Better long term RV precision is expected.

## 4. All Sky Doppler Surveys for Extra-solar Planets

Observations with ET have demonstrated that the dispersed fixed-delay interferometer is suitable for multiple object observing and can provide high throughput and high Doppler precision. Multi-object and high throughput Doppler surveys are important since they can significantly speed up the detection of extrasolar planetary systems and also reduce the cost. Currently, only about 4000 solar type stars are being searched by half a dozen telescopes using the echelles. Based on current planet detection rate of 10% of solar type stars, $``$ 400 planets will be discovered in the next $``$ 10 years. These surveys are time-consuming and very costly. This is because current echelles can only measure a single object per exposure, and because only relatively bright stars can be observed since the instrument has low detection efficiency, whereas high photon flux is required for precision measurements. Therefore, a multiple object RV survey can significantly improve the current situation.

Considering there are not many bright stars in the sky for the measurements, the multiple object survey must include relatively faint stars within the telescope field of view (FOV). Based on our estimation from previous star count surveys in the visible (Bahcall & Soniera 1980) and a dust extinction map (Schlegel et al. 1998), we find there are about a half million stars from A – M types brighter than V = 10, and about 4 millions of stars are brighter than V = 12. On average, about 100 stars with V magnitude brighter than 12 mag. are within a 1 deg field-of-the-view. Therefore, in order to conduct an efficient all sky survey for extrasolar planets with a multiple object dispersed fixed-delay interferometer, we need to reach about V = 12 with a wide field telescope.

This all sky survey is possible with modern wide field telescopes such as the Sloan 2.5-m and the WIYN 3.5m telescopes. The Sloan has a 7 degree FOV and the WIYN has a 1 degree FOV. Since both of the telescopes are moderate size, high throughput becomes critical to reach high Doppler precision for faint stars. High throughput can be achieved with an optimally designed instrument. The prototype at the KPNO 2.1-m has already demonstrated good throughput of 3.4% under 1.5 arcsec seeing. In the future, a factor of 4 times improvement should be achievable, including 2 times in the interferometer transmission by feeding both interferometer outputs to the spectrograph, 1.4 times in the étendu match through updating current spectrograph with a faster spectrograph (e.g., $`f`$/2 instead of $`f`$/7.5), and 1.5 times in the spectrograph transmission by using a higher efficiency grating (such as a VPH grating with $``$ 85% efficiency instead of current reflection grating with $``$ 55% efficiency). A total estimated detection efficiency of $``$15% can be expected. With this detection efficiency, a Doppler precision of $``$ 15 m/s can be reached for a V = 12 star within an hour integration and $``$ 1000 Å wavelength coverage. Therefore, most of the stars brighter than V = 12 can be surveyed with high Doppler precision for planet candidates at 2 m class wide field telescopes with multi-object RV instrument. The candidate systems can be further studied with a single object high throughput dispersed fixed-delay interferometer with higher precision at larger aperture telescopes to look for additional planet members.

In the future, in order to take full advantage of the potential of the multi-object dispersed fixed-delay interferometer for all sky RV surveys, we need to simultaneously feed multi-object fiber beams into three instruments which have best sensitivity in the near-IR, visible and near-UV, respectively. This design allows optimal match of the interferometer sensitivity with the peak of the spectral flux and lines. For instance, a late M type star has about 10 times more flux in the near-IR than in the visible and has many molecular and atomic lines for precision Doppler measurements (Kirkpatrick et al. 1993). Because the line width is very different from early type to very late types, the optical delay will be changed accordingly to minimize Doppler errors.

To fully achieve high Doppler precision with this instrument, a large waveband is very important since this will allow the capture of more photons from stars for the measurements. In the prototype, the wavelength coverage is about 200 Å. By simply increasing wavelength coverage to $``$ 1000 Å, more than a factor of 2 times better Doppler precision can result. Doppler precision can be possibly further reduced by using reference sources with higher visibility than the iodine absorption. Our measured fringe visibility for typical Thorium lines in our Palomar data is about 50%, resulting in intrinsic Doppler error well below 1 m/s. Current RV measurements using ThAr calibration in Dr. Mayor’s group have already achieved $``$ 5 m/s precision and with further improvement such as double fiber mode scrambling and vacuum operation, they believe that they can reach $``$ 1 m/s in the new HARPS high resolution echelle spectrometer (Queloz 2002 private communications; Pepe et al. 2000). We believe we may be able to achieve similar calibration precision by using similar procedures. If this technique works, the observing efficiency can be tripled over that of our current iodine-based technique since there will be no photon loss due to iodine absorption (typical loss by iodine absorption $``$ 30-40%) and also there is no need to create a separate stellar template for each observation. This will be a big plus for the survey.

In summary, to have a very successful all sky survey requires multi-object observation capability, high Doppler precision and high throughput. A wide field telescope with at least three multi-object dispersed fixed-delay interferometers optimized for near-UV, visible and near-IR wavebands can provide the best sensitivity for detecting thousands of planets around stars from very early types to very late types in the near future.

### Acknowledgments.

The authors are grateful to Larry Ramsey, Don Schneider, Eric Feigelson, Steinn Sigurdsson, Tom Soltysinski, Ron Reynolds, Fred Roesler, Didier Queloz, Harvey Moseley, Bruce Woodgate, Roger Angel, Mike Shao, Chas Beichman, Bill Cochran, Wes Traub, Ed Jenkins and Jim Gunn for useful discussions. We acknowledge support from the HET staff, the Palomar observatory staff and KPNO staff. This work is supported by NASA JPL and Penn State Eberly College of Science. A portion of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

## References

Bahcall, J.N. and R.M. Soniera. 1980, ApJS, 44, 73

Baranne, A. et al. 1996, A&AS, 119, 373

Beckers, J.M., & Brown T.M., 1978, Osser Mem Astrophys Obs Arcetri. No. 106, 189

Bouchy, F., Pepe, F., & Queloz, D. 2001, A&A, 374, 733

Butler, R.P., Marcy, G.W., Williams, E., McCarthy, C., Dosanjh,P., & Vogt, S.S. 1996, PASP, 108, 500

Dravins, D. 1987, A&A, 172, 200

D’Odorico, S. et al. 2000, Proc. SPIE, 4005, 121

Erskine, D.J., & Ge, J., 2000, in Proc. Imaging the Universe in Three Dimension, Edited by W. van Breugel and J. Bland-Hawthorn ASP Conference Series, 195, 501

Ge, J., 2002, ApJ, 571, L165

Ge, J., Erskine, D.J., & Rushford, M., 2002a, PASP, 114, 1016

Ge, J., et al. 2002b, Proc. SPIE 4835, in press

Gorskii, S.M., & Lebedev, V.P., 1977, Izv Krym Astrofiz Obs 57, 228

Kirkpatrick, D.J., et al. 1993, ApJ, 402, 643

Kozhevatov, I.E., Kulikova, E.Kh. & Cheragin, N.P. 1995, Astronomy Letters, 21, 418

Kozhevatov, I.E., Kulikova, E.H., & Cheragin, N.P. 1996, Solar Physics, 168, 251

Mayor, M., & Queloz, D. 1995, Nature, 378, 355

Pepe, F. et al. 2000, Proc. SPIE, 4008, 58

Schlegel, D.J., D.P. Finkbeiner and M. Davis. 1998, ApJ, 500, 525

Vogt, S.S. et al. 1994, Proc. SPIE, 2198, 362

Vogt, S.S., Marcy, G.W., Butler, R. P., Apps, K., 2000, ApJ, 536, 902


format: text

text:
# Searching for Bulges at the End of the Hubble SequenceBased in part on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal # 8599.

## 1 Introduction

The question whether the morphology of galaxies is imprinted by the initial conditions of their formation or rather determined by secular evolution remains a subject of intense debate. The existence of the Hubble sequence has for many years provided important constraints on this issue. In very simple terms the Hubble sequence tells us that galaxies are made up of two components: a bulge and a disk. The canonical view of these components has long been that bulges have $`R^{1/4}`$ surface brightness profiles (de Vaucouleurs, 1948) while disks have exponential surface brightness profiles. As one goes from early-type to late-type galaxies one goes from galaxies that are bulge-dominated to galaxies that are disk-dominated. While this simplistic interpretation of the Hubble sequence has definite value, reality is considerably more complicated.

In recent years, our views of the Hubble sequence have evolved and gained more nuance. For elliptical galaxies it has become clear that they are not necessarily pure bulge systems: many elliptical galaxies contain embedded disks. There is evidence from other information (e.g., kinematics) that elliptical galaxies form a heterogeneous class of galaxies that may have formed in different ways (e.g. Kormendy & Bender, 1996). For spiral galaxies a clearer understanding has developed of their bulge properties. High-resolution imaging - both from the ground (e.g. de Jong, 1996a), and with the Hubble Space Telescope (HST, Carollo et al., 1998) \- has shown that the central surface brightness profile (SBP) of many late-type spirals cannot be fit by the classical $`R^{1/4}`$ law that is well suited to describe the bulge profiles of early-type spirals. Instead, the SBPs of many late-type spirals rise above the extrapolation of the exponential disk in a way that can be well described by a second exponential (Andredakis & Sanders, 1994). This has led to the now popular view that spiral bulges come in two flavors: on the one hand, the classical $`R^{1/4}`$ bulges which are mostly observed in early-type spirals, and on the other the “pseudo-bulges” (Kormendy, 1993) or “exponential bulges” (Carollo et al., 1998) which are prevalent in later Hubble types. In reality there is probably a continuum of properties, instead of a dichotomy. When $`R^{1/n}`$ profiles (Sérsic, 1968) are fit to available SBPs, the profile shape parameter spans the full range of values $`1n6`$; the profile shape parameter correlates with both Hubble type and bulge-to-disk ratio of the galaxy111Obviously, Hubble type and bulge-to-disk ratio are themselves strongly correlated through the definition of the Hubble sequence., in the sense that spiral galaxies with earlier Hubble type have bulges with higher $`n`$ values (Andredakis, Peletier, & Balcells, 1995; Graham, 2001).

The existence of different types of bulges in disk galaxies can be plausibly explained in the context of popular scenarios for the formation and secular evolution of galaxies. The classical massive $`R^{1/4}`$ law bulges fit in with the “primordial collapse” formation scenario first suggested by Eggen, Lynden-Bell, & Sandage (1962), in which the bulge forms during the initial collapse of a galaxy-sized density perturbation, and later “acquires” a disk through accretion processes. By contrast, the pseudo-bulges may have formed by secular evolution of a pre-existing disk, so that they formed after the disk, out of disk material. Some support for this scenario comes from the fact that pseudo-bulges are dynamically similar to their host disks (Kormendy, 1993). Plausible secular evolution scenarios include the accretion of satellite galaxies (Aguerri et al., 2001), buckling instabilities in a stellar bar (Raha et al., 1991), and the disruption of a stellar bar through the accumulation of a central mass concentration (Norman, Sellwood, & Hasan, 1996). Many discussions of these and related topics can be found in the review by Wyse, Gilmore, & Franx (1997) and in the proceedings of the recent workshop on ‘The Formation of Galactic Bulges’ (Carollo, Ferguson, & Wyse, 1999).

In the present paper we study the presence and properties of bulges in the very latest-type spiral galaxies (Scd or later). This is an important topic for several reasons. First, these galaxies are generally classified as very late type spirals because they do not have a very prominent bulge. As a result, many observational studies of bulges have avoided these galaxies. Second, it has become clear from recent work with HST that the majority of spiral galaxies contain a central star cluster. In the very latest-type spiral galaxies we find that $`75`$% of the galaxies contain such a star cluster (Böker et al., 2002, hereafter paper I). In late Hubble types, these clusters are easily mistaken for a small bulge when observed from the ground, even in good seeing conditions. So the bright, compact “bulges” in late-type spirals which were used as a classification criterion in the original work of Hubble (1926, 1936) may in fact be dense star clusters occupying the photocenter of the galaxy. The purpose of this paper is to shed some light on these issues. In particular, we investigate whether the very latest-type spirals are completely bulgeless, whether they show excess light above the constant scale-length disk, and if so, whether this in fact implies the presence of a separate entity which could rightfully be called a bulge.

HST resolution is needed to separate the luminous nuclear star cluster from a putative bulge. Our I-band snapshot survey of late-type spiral galaxies conducted with the Wide Field and Planetary Camera 2 (WFPC2) and discussed in Paper I therefore forms the basis of our analysis. We complement the HST observations with ground-based data that extends to larger radii. The paper is organized as follows: in § 2, we describe the data and the analysis methods that form the basis of our work. The results of our analysis are summarized in § 3. We discuss the implications of our findings, and present our conclusions in § 4.

## 2 Data

### 2.1 HST Data

The target list and selection criteria for the HST sample are described in paper I. In brief, we selected 113 nearby ($`\mathrm{v}_{\mathrm{helio}}<2000\mathrm{km}/\mathrm{s}`$) galaxies with low inclination ($`\mathrm{R}_{25}\mathrm{log}(\mathrm{a}/\mathrm{b})<0.2`$) and Hubble type between Scd and Sm ($`6<\mathrm{T}<9`$). 77 galaxies have been observed to date with WFPC2 on board HST in the F814W filter which is similar to Johnson I-band. In all cases, the integration time was $`600s`$. After standard data reduction to correct for instrumental effects and cosmic rays, we performed an isophotal analysis using the ellipse task in IRAF. For all galaxies studied in the present analysis, the ellipse center, ellipticity, and position angle were allowed to vary freely from isophote to isophote.

The WFPC2 SBPs have been published in Figures 1 and 2 of paper I. In this study, we discuss the 19 galaxies for which we were able to find ground-based, wide-field surface brightness data for comparison. Our focus lies on the structural properties of the disks rather than the photometry of the nuclear clusters (which was the subject of paper I).

### 2.2 Ground-Based Surface Photometry

Due to the short integration time ($`600\mathrm{s}`$) of the WFPC2 snapshot observations, the signal-to-noise ratio is generally insufficient to make use of the WF chips in our analysis. It is therefore prudent to ask whether the rather limited field-of-view of the PC chip yields a representative view of the galaxy disk. In order to address this question, we searched the literature for published wide-field observations of our sample galaxies. In general, these are scarce, most likely due to the less than spectacular star formation activity and rather low surface brightness of the latest Hubble type spirals. In this context, it is worth pointing out that all-sky surveys such as 2MASS do not reach the required sensitivity to map the outer disk of most of our sample galaxies, and hence are not useful for our purposes.

We found published wide-field SBPs for 19 of our sample galaxies which are listed in Table 1. All but one (NGC 4904) of these 19 galaxies harbor a nuclear cluster as described in paper I. The wide-field surface brightness data for 18 galaxies are based on images in the Photometric Atlas of Northern Bright Galaxies (PANBG, Kodaira et al., 1990). The digitized profiles, which have been studied in detail by Baggett, Bagett, & Anderson (1998), were kindly provided to us in machine readable form by M. Hamabe. The PANBG images were taken in the V-band, with SBPs derived from cuts along the major and minor semi-axes of the galaxy, typically extending out to radii between $`100^{\prime \prime }`$ and $`250^{\prime \prime }`$. For comparison to the high-resolution WFPC2 data, we use (in most cases) the average of the two semi-major axis cuts. For some galaxies, visual inspection of the images and profiles in the PANBG showed that one or both semi-major axis cuts appeared affected by large-amplitude, small-scale variations. In these cases, we either used the remaining semi-major axis cut, or the average of the semi-minor axes, scaled by the ellipticity of the envelope ellipse (the elliptical fit to the isophote contour at $`\mu _\mathrm{V}=25`$) listed in column 13 of the PANBG. The remaining galaxy (NGC 2805) was part of the sample of Kent (1984), who used a red passband centered at 6500Å.

Figure Searching for Bulges at the End of the Hubble Sequence11affiliation: Based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal # 8599. shows the available data for the 19 galaxies. For each galaxy, the left panel contains the WFPC2 I-band SBPs; the right panel contains the PANBG data for each galaxy. The ordering of the galaxies in the figure and the solid lines in the panels are discussed in § 3.1 below. In order to normalize both datasets, we shifted the (V-band) PANBG profiles by the amount needed to match the average I-band surface brightness in the region between $`8^{\prime \prime }`$ and $`12^{\prime \prime }`$. Typical offsets were in the range $`1.0VI1.5`$. This is not unreasonable, given that the typical integrated color of galaxies of these Hubble types is $`VI=1.0\pm 0.2`$ (e.g. de Jong, 1996b). The colors that we infer appear to be somewhat too red, probably indicating a systematic offset between the two datasets in their photometric calibration. However, in this study, we are interested mainly in the structural properties of the light distribution and this problem therefore does not affect our analysis.

There are a few limitations to the combined data set thus obtained. For one, the PANBG profiles were derived from simple cuts along the galaxy major or minor axis, rather than from azimuthally averaged isophote fitting methods (which is what we used for the HST data). The PANBG profiles are therefore affected by spiral arms, dust lanes, star forming regions, or other asymmetries, even though we made every attempt to minimize these by selecting the least-affected semi-axes. Also, possible mass-to-light variations or color gradients could complicate the direct comparison of the two SBPs. However, the results of de Jong (1996b) and Matthews & Gallagher (1997) have shown that color gradients in late-type disks are generally small, with variations in (B-I) of $`<0.2`$mag over a few scalelengths. It thus appears that - except for the very center - the stellar populations of late-type disks are rather uniform, and our approach of profile matching seems justified. In summary, the PANBG profiles should yield a fairly robust impression of the SBP of the outer disk. They are certainly useful for comparison to our WFPC2 images, especially to test whether the central $`15^{\prime \prime }`$ are representative for the profile of the stellar disk.

## 3 Results

### 3.1 Surface Brightness Profile Shapes

Although the galaxies in our sample were selected as a fairly homogeneous group with respect to Hubble type, distance, and inclination, the SBPs show a wide range of shapes. To address this quantitatively, we started by fitting for each galaxy an exponential profile to the I-band SBPs inside the PC field. We will refer to this as the ‘inner exponential fit’. In this fit we always excluded the nuclear cluster, as well as any ‘central excess emission’ (the definition and nature of which is discussed below), if present. The fits are shown as solid lines in the left panel of each plot. We show the same exponential fit also in the right panel of each plot (so the solid line in the right panel is not a fit to the PANBG data, but only the continuation of the fit in the left panel).

The extent to which the inner exponential fit describes the data at the different radii differs considerably for the different galaxies. To illustrate this point, we have ordered the galaxies in Figure Searching for Bulges at the End of the Hubble Sequence11affiliation: Based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal # 8599. in a rough sequence, according to the following scheme. The sequence starts with galaxies for which there is no evidence for a bulge of other inner component in the SBP, except for the nuclear star cluster (which has a typical HWHM radius of $`510\mathrm{pc}`$; paper I). This group itself is somewhat heterogeneous. For some galaxies the inner exponential fit yields a good description of the large-scale stellar disk as measured by the PANBG profiles (e.g., NGC 428). These galaxies are fairly well described by a single exponential. For other galaxies (e.g., NGC 4299), the outer disk falls below the inner exponential fit, which implies that the SBP at small radii falls below the inward extrapolation of the outer disk. This is exactly opposite to what one would expect if a bulge were present (see § 3.2 below). We interpret both types of galaxy as pure disk systems. Interestingly, these pure disks have SBPs that can differ from a single exponential.

Progressing along the sequence, the profiles show an increasing amount of light in excess of the inner exponential fit on radial scales of a few hundred pc (the affected radii were excluded from the inner exponential fit). We refer to this light as ‘central excess emission’. This emission has a distinct profile shape: while the overall galaxy profile steepens towards the center, the excess component flattens towards the center (see, e.g., NGC 275, NGC 2139, or NGC 3346). This is different from what one would expect for a bulge. Bulges tend to have $`R^{1/n}`$ profiles with $`n>1`$ (Andredakis, Peletier, & Balcells, 1995; Graham, 2001). Such profiles steepen towards the center, when plotted as a function of (linear) $`r`$. Consequently, we interpret the central excess emission as a phenomenon that is different from a central bulge. We discuss this component further in § 3.3.

Approximately halfway along the sequence the outward extrapolation of the inner exponential fit begins to underpredict considerably the true brightness of the outer disk (e.g., NGC 5068). Conversely, this implies that there is excess light in the central few kpc over the inward extrapolation of the outer disk. This is the component that is traditionally called a bulge. We discuss the nature of this component further in § 3.2.

Given the narrow range of Hubble types of our sample (between Scd and Sm), and the uniformity in inclination (less than $`40\mathrm{deg}`$ from face-on), the variety of profile shapes is somewhat surprising. In addition, there seems to be no clear trend of position along our sequence with either Hubble type or galaxy luminosity. For example, the profiles of NGC 275 and NGC 1042 are very different, even though both galaxies have nearly identical Hubble type (SB(rs)cd vs. SAB(rs)cd) and total blue magnitude ($`\mathrm{M}_\mathrm{B}=18.8`$). Clearly, late-type spirals are not simple systems, their morphologies (and possibly evolutionary states) can differ drastically.

### 3.2 Disks and Bulges

A considerable source of confusion in discussions of bulges is that different definitions are often used when referring to a bulge. Modern theorists tend to think of a bulge as a kinematically hot component with an extended three-dimensional structure. However, the observational definition of a bulge is often a different one, and relies on the assumption that bulges have different SBPs than disks. The standard wisdom is that a disk can be well fit by an exponential profile, $`I(r)e^{r/r_d}`$. By contrast, bulges can be described empirically with either a second exponential with different scale length (e.g. de Jong, 1996a), a de Vaucouleurs $`R^{1/4}`$ profile (e.g. Boroson, 1981; Baggett et al., 1998), or most generally, a Sérsic $`R^{1/n}`$ profile (e.g. Graham, 2001).

If one defines a bulge as ‘the central concentration of mass \[or light\] in excess of the inward extrapolation of the outer, constant scale-length, exponential disk’ (Carollo, Ferguson, & Wyse, 1999), one must make sure that exponential profiles indeed provide a highly accurate description of galaxy disks. Here one is on shaky ground. The exponential model has been widely used to describe galaxy disks because it provides a fairly good fit to many observed galaxy profiles (e.g. Freeman, 1970; Kent, 1984; de Jong, 1996a). However, the fits are generally restricted to the region outside the central few kpc, so that it is actually very hard to know whether disks remain exponential all the way into the center. The disks of many late-type spirals in fact show sharp, well-defined “breaks” in their brightness profiles, which separate sections of the disk that follow exponential profiles with different scale lengths (Bosma & Freeman, 1993; Matthews & Gallagher, 1997). So observationally, it is not at all clear how photometry alone can distinguish between an exponential bulge and a somewhat steeper disk section in the inner part of the galaxy.

Theoretical considerations also do not clinch the argument. There is no theory that unequivocally predicts that disks must be exponential. As an ‘a posteriori’ explanation, (semi-)analytical models that invoke angular momentum redistribution via viscosity-driven radial gas flows can provide a plausible mechanism for building profile shapes close (but not identical) to an exponential. This has been demonstrated in the initial work by Lin & Pringle (1987) and Yoshii & Sommer-Larsen (1989), and recently confirmed by the refined study of Ferguson & Clarke (2001). An alternative scenario which invokes self-propagating stochastic star formation in a disk of constant density atomic hydrogen gas (Seiden, Schulman, & Elmegreen, 1984) also results in brightness profiles that approximate the exponential shape. However, all these models are dependent on a number of parameters, which can be fine-tuned to reproduce a range of profile shapes. In fact, most model profiles presented by Yoshii & Sommer-Larsen (1989) show significant curvature at all radii. It is therefore not obvious that a simple exponential is a good analytic model for disks to begin with.

In order to demonstrate the ambiguity of profile decompositions, and the ensuing interpretation of the derived “bulge” properties, we have compared two different fits to the combined (WFPC2 plus PANBG) datasets, excluding again the nuclear cluster. The two models are: (a) a single Sérsic profile; and (b) the sum of two exponentials of different scale lengths. Interestingly, we find that most galaxies in the second half of our sequence are equally well fit with either model, thus questioning the need for a “bulge” component. This is demonstrated in Figure Searching for Bulges at the End of the Hubble Sequence11affiliation: Based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal # 8599. which compares these two fits for five galaxies with apparent disk-bulge transitions. In most cases the two profiles provide equally adequate fits, and for one galaxy, NGC 5668, the Sérsic profile provides an obviously better fit. Only for one galaxy, NGC 1042, does a sum of two exponentials provide the better fit.

Overall, we find that roughly half of our sample galaxies show excess light (over the inward extrapolation of the outer exponential disk) that can be fit by a second exponential. This is similar to the ‘pseudo bulges’ or ‘exponential bulges’ that have been reported by other authors, mostly in earlier type spirals. However, the extent to which this excess light is associated with a bona-fide bulge component remains an open question. The SBPs of these galaxies are generally well fit with Sérsic profiles with $`n1`$$`2.5`$ (as listed in the panels of Figure Searching for Bulges at the End of the Hubble Sequence11affiliation: Based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal # 8599.). In the absence of information on kinematics or three-dimensional structure it is difficult to rule out that we are dealing merely with non-exponential disks, or possibly small-scale stellar bars which are known to also have exponential-like SBPs (e.g. Elmegreen & Elmegreen, 1985).

### 3.3 The central light excess

As described in § 3.1, several galaxies of our sample show evidence for central excess emission on scales of a few hundred pc which does not show the steepening characteristic of bulges. In order to investigate the nature of this emission, we have measured the amount of excess light for four example cases (NGC 275, NGC 2139, NGC 3346 and NGC 5584). We first determined by eye the radius $`r_{ap}`$ at which the SBP starts to deviate from the inward extrapolation of the exponential fit. We then analyzed the WFPC2 I-band images presented in paper I with the IRAF photometry package apphot in order to measure the total flux within a circular aperture of radius $`r_{ap}`$, centered on the nuclear star cluster. To give a better visual impression of the four case studies, Figure Searching for Bulges at the End of the Hubble Sequence11affiliation: Based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal # 8599. reproduces their images presented in paper I. Here, we have overlaid circles with radius $`r_{ap}`$ to indicate the aperture over which we measure the excess emission. Also, we have slightly modified the grey scale stretch to emphasize the disk structure rather than the nuclear cluster.

From the total flux inside $`r_{ap}`$, we subtract the contribution of the galaxy disk which is calculated by integrating the exponential fit (the solid lines in Figure Searching for Bulges at the End of the Hubble Sequence11affiliation: Based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal # 8599.) inside $`r_{ap}`$. We also subtract the light from the nuclear cluster as listed in paper I, which usually constitutes only a small correction. Table 2 summarizes the results of this analysis. For the four galaxies of Figure Searching for Bulges at the End of the Hubble Sequence11affiliation: Based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal # 8599., the apparent magnitude of the excess emission is in the range $`15.2m_I17.4`$ which corresponds to luminosities between $`210^7\mathrm{L}_{}`$ and $`210^8\mathrm{L}_{}`$. This is between 10 and 100 times brighter than the median luminosity of the nuclear star cluster in our sample ($`\mathrm{M}_\mathrm{I}=11.5`$ or $`1.610^6\mathrm{L}_{}`$, paper I), but constitutes only between 1% and 5% of the total galaxy luminosity. It is interesting that with absolute luminosities of $`16.7M_I14.2`$ (using the distances in Table 2), this excess emission appears to be a faint end continuation of the distribution of spiral bulge luminosities shown in Figure 14 of Graham (2001), if one assumes a mean color of $`\mathrm{B}\mathrm{I}1`$.

The images in Figure Searching for Bulges at the End of the Hubble Sequence11affiliation: Based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal # 8599. reveal a variety of morphologies for the excess emission. For NGC 275, NGC 5584, and NGC 3346, the emission appears fairly smooth and follows the larger scale disk structure. However, this is clearly not the case for NGC 2139, for which the observed excess emission (which is also by far the brightest in our sample) is due to a prominent, bright star forming region which is elongated in the E-W direction. Inspection of the WF chips shows that this structure is the inward continuation of at least one spiral arm (the other direction is off the WFPC2 field of view). This underlines the danger of relying solely on one-dimensional surface brightness data for evaluating galaxy morphology.

## 4 Summarizing Discussion

We have presented an investigation into the structural properties of 19 spiral galaxies with Hubble type between Scd and Sm. From a combination of our high-resolution HST data and wide-field ground-based images, we obtain surface brightness profiles (SBPs) for the sample galaxies that cover a large dynamic range in galactic radius. We use these profiles to study quantitatively the presence and properties of bulges in spiral galaxies of the very latest Hubble types. Of course, by the very definition of the Hubble sequence, one does not expect very prominent bulges in these galaxies. However, it has been somewhat of an open question whether bulges are present at all in these galaxies. It has been realized only recently that space-based resolution is required to properly address this question. The nuclear morphologies of most spiral galaxies are complex, and a large fraction of spiral galaxies has a nuclear star cluster. These clusters are easily mistaken for “compact bulges” when observed with ground-based resolution. Previous studies with the HST have focused in majority on earlier Hubble types, and the present study is the first to focus exclusively on the very latest Hubble types.

Approximately 30% of the sample galaxies seem to be more-or-less “pure” exponential disks, without any type of a central bulge. Our sample was selected to focus on face-on galaxies, but studies of edge-on galaxies support the view that disk galaxies can indeed be completely bulgeless. Matthews, Gallagher, & van Driel (1999) and Matthews (2000) have studied edge-on “super-thin” galaxies such as UGC 7321, and have demonstrated that in these galaxies there is no evidence for a spheroidal component.

Despite being disk-dominated systems, most galaxies in our sample have SBPs that cannot be well fit by a single exponential, in the sense that the surface brightness in the central few kpc exceeds the inward extrapolation of the outer exponential disk. This has generally been found for other samples of spirals as well, and has generally prompted the addition of a bulge component to analytic models of the SBP. In particular, numerous studies have shown that the SBPs of intermediate- to late-type spirals can be well fit by a sum of two exponentials. The inner exponential has traditionally been interpreted as an ‘exponential’ bulge or ‘pseudo’ bulge. Such bulges can be (qualitatively) explained theoretically as a result of secular evolution of the disk. While this may be correct, we point out that this is not an unambiguous interpretation. In the absence of information on three-dimensional structure or dynamics there is no guarantee that one is dealing with a bona-fide bulge component. We explicitly illustrate this point for the galaxies in our sample that are not well fit by a single exponential. Indeed, the SBPs of these galaxies are generally well fit by a sum of two exponentials. However, most profiles can be described at least equally well with a Sérsic-type $`R^{1/n}`$ model over the entire radial range (outside the nuclear star cluster). The shape parameter $`n`$ is in the range $`1<n<2.5`$, which is not unrealistically large. As we have discussed in §3.2, there is no a priori theoretical reason for the SBPs of disk galaxies to be pure exponentials; Sérsic-type profiles are in some sense equally arbitrary as the model of choice. So it may well be that we are simply dealing with non-exponential disks in most of our sample galaxies, and that bona-fide bulges are rare at these Hubble types.

We have found that a number of late-type galaxies show central excess emission on spatial scales of a few hundred parsec that can be attributed neither to the nuclear cluster, nor to the Sérsic-type description of the stellar disk, nor to what one would generally consider to be a bulge component. The origin of this light component remains unclear.

One of the interesting findings from our work is that, despite the narrow range in Hubble type, the SBPs of the sample galaxies are far from uniform. Our study finds no systematic trends in the structural properties with morphological type. The exact Hubble type of spirals between Scd and Sm appears somewhat arbitrary, in the sense that it provides little information about the presence and relative importance of galaxy bulges. It is quite possible that this may be due to the presence of nuclear clusters, which may have played an important role in the morphological classification in photographic catalogs. In general our results fit in well with the picture that emerged from a ground-based imaging study of 49 late-type spirals by Matthews & Gallagher (1997). They found that late-type spirals exhibit a diverse array of structural properties and morphologies, even in galaxies with otherwise similar parameters, and they concluded that bulges are often very weak or non-existent. In studies of this kind one does need to be concerned about selection effects and observational bias. The Matthews & Gallagher (1997) sample was selected for low galaxy luminosity ($`M_V>18.8`$), and one might worry that it did not provide a representative view of the family of very late-type spirals. Our galaxy sample, on the other hand, has been selected only for Hubble type, distance, and inclination. It thus includes more luminous (and presumably more massive) galaxies than the Matthews & Gallagher (1997) sample. The blue absolute luminosity of the galaxies in our sample ranges from $`M_B=17.4`$ (NGC 2552, NGC 4701) to $`M_B=20.5`$ (NGC 2805). However, we find no systematic correlation between absolute luminosity and position along the sequence of SBP shapes discussed in §3.1. This suggests that more luminous late-type spirals are not systematically different from their faint-end relatives, and that presumably the Matthews & Gallagher (1997) results are valid for most late-type spirals.

Late type spiral galaxies are in most ways “normal” spiral galaxies. Their angular momenta and rotation velocities are not atypical, as demonstrated by a recent survey of optical rotation curves (Matthews & Gallagher, 2002). They are in some sense the dynamically simplest type of disk galaxies. They are (mostly) disk-dominated and often have only faint spiral arm structure or even no detectable density perturbations at all. Yet, we are clearly a long way from understanding their formation and evolution. The latest-type disk galaxies provide the most stringent observational constraints on the well-known angular momentum problem in cold dark matter (CDM) galaxy formation models (e.g. Navarro & Steinmetz, 2000). Because of their unevolved disks and apparent history of only modest star formation, they are a challenge for proposed solutions for the angular momentum problem which rely on energy feedback from supernova explosions (Silk, 2001). To gain a better understanding of these issues it will be important to continue to improve our knowledge of the structure of late-type spiral galaxies. Observational studies such as the one presented here will continue to be essential.

We are grateful to M. Hamabe for providing us with the digital surface brightness data from the PANBG, to Roelof de Jong for helpful discussions, and to the anonymous referee for useful comments. R. S. acknowledges a grant from the STScI summer student program. Support for proposal # 8599 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association for Research in Astronomy, Inc., under NASA contract no. NAS 5-26555. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. It has also benefited greatly from use of the Lyon-Meudon Extragalactic Database (LEDA, http://leda.univ-lyon1.fr).


text:
# Particle acceleration in ultra-relativistic parallel shock waves

## 1 Introduction

There are three distinct astrophysical situations where the bulk plasma flow is extremely relativistic with Lorentz factors $`\mathrm{\Gamma }=(1V^2/c^2)^{1/2}10`$. These are some Active Galactic Nuclei (AGN) jet and hot spot sites, Gamma Ray Burst (GRB) fireballs and pulsar polar winds. Within these flows relativistic shocks appear and diffusive particle acceleration takes place. In each case, there is evidence for energetic particle acceleration to high $`\mathrm{\Gamma }`$ factors but the upper limit to the possible energy attained becomes less certain with increasing bulk flow velocity. For AGN jets, $`\mathrm{\Gamma }10`$ plausibly results in the known $`10^{19}`$ eV cosmic rays (Quenby and Lieu, 1989) via diffusive shock acceleration at the termination shocks. However, the more direct synchrotron/self-Compton modeling only reveals 10 GeV gamma-rays (e.g. Maraschi et al., 1992). For GRB fireballs, $`\mathrm{\Gamma }10010^3`$ appears to eventually produce gamma-rays at least up to 100 MeV (Fenimore et al., 1993) but at such flow speeds, it is not clear whether repeated shock crossings are possible and the predictions of $`10^{19}`$ eV neutrinos (Vietri, 1998), which would be direct evidence of a diffusive like process have yet to be verified. Similarly, pulsar winds could result in a $`\mathrm{\Gamma }10^4`$ flow encountering the nebula envelope (Hoshino et al., 1992), but there is no evidence for more than TeV acceleration anywhere in the system. However Lucek and Bell (1994) show that combined motion across a near perpendicular shock with a non-uniform field and in the electromagnetic field of the nebula may produce much higher energies.
Past work has suggested that if the shocks are close to parallel so that an $`𝐄=0`$ frame can be defined, diffusive shock acceleration is in principle possible and provided large angle scattering occurs in the plasma frame, each shock crossing cycle yields a $`\mathrm{\Gamma }^2`$ energy increase (Quenby and Lieu, 1989). This result has been exploited by Vietri (1995,1998) modeling GRB fireball particle shock acceleration. It is the purpose of this paper to present simulation results on particle acceleration with $`\mathrm{\Gamma }`$ values extended up to $``$ $`10^3`$. Here we concentrate on the parallel shock configurations ($`\psi =0^o`$). Our results, measured in the shock frame, will be compared with the standard spectral prediction of non-relativistic theory, where the momentum spectrum is $`n(p)p^\alpha `$ with $`\alpha =(r+2)/(r1)`$, for shock compression ratio $`r=V_1/V_2`$, independent of scattering details, and $`\alpha =2`$ for strong unmodified shocks with $`r=4`$. Analytically it has been shown that over a limited momentum range, the non-relativistic time scale for acceleration is $`T=[c/(V_1V_2)][\lambda _1/V_1+\lambda _2/V_2]`$ where, ’1’ and ’2’ refer to upstream and downstream and $`\lambda `$ is the scattering mean free path. The shock jump solutions of Ballard and Heavens (1991) and the review of Kirk and Duffy (1999) show that, a range of shock frame $`r`$ values largely limited between 3 and 4, typically arise for oblique, relativistic shocks, whereas $`r`$=4 is the extreme relativistic hydrodynamic limit. Hence for simplicity we adopt $`r`$=4 but, in the companion paper dealing with oblique shocks we demonstrate the insensitivity of our chief conclusions on the exact value. For our parameters, $`T=(8/3)(nr_{g,1}/c)`$ where, $`\lambda =nr_g`$ and $`r_g`$ is Larmor radius.
The first relativistic correction noticed to the non-relativistic theory, was the spectral flattening seen in parallel shocks (Kirk and Schneider, 1987a,b), followed by the discovery of a speed-up in the acceleration rate (Quenby and Lieu, 1989; Ellison et al., 1990). These results have subsequently been extended to non-parallel and non-linear sub-luminal shocks (Lieu, et al., 1994 and Bednarz and Ostrowski, 1996 who use $`T=r_g/c`$ as the unit of acceleration time), where there is large angle scattering or pitch angle diffusion. Most recently, much interest has been focused on differences occurring between small angle and large angle scattering models for the test particle-turbulence interactions in the fluid frame. Gallant and Achterberg (1999) suggest that for a field model consisting either of randomly orientated uniform field cells or a uniform field, the scattering is limited to $`\delta \theta <1/\mathrm{\Gamma }`$ in the upstream fluid frame. Generally, such a model yields a test particle differential spectral slope $`2.2`$ (Baring, 1999). Furthermore, on this model, Gallant and Achterberg (1999) provided an analytical demonstration that in the high $`\mathrm{\Gamma }`$ limit, only the first crossing cycle exhibited a $`\mathrm{\Gamma }^2`$ energy increase, while in all other crossings the energy is only doubled and hence one may expect the relativistic speed-up to be limited to this crossing. The computations of Baring (1999) also showed a much reduced energy enhancement in subsequent cycles for the small angle scattering case. Achterberg et al. (2001) extend the conclusions of Gallant and Achterberg (1999) to a wider variety of upstream field models while, Kirk et al. (2000) employ an eigenfunction method to deal with the anisotropies in diffusion when the scattering depends on the angle with the shock normal. A recent general review of particle acceleration and relativistic shocks is given by Kirk and Duffy (1999). However, the detailed contrast between the effects of various particle scattering models and the behaviour of these models at $`\mathrm{\Gamma }`$ factors $`10^3`$ remain to be elucidated, both for parallel and quasi-perpendicular shocks. Hence this paper starts a series of investigations in the ultra-relativistic limit with a consideration of parallel shocks. In section 2 we present the numerical method used, section 3 presents the results while in section 4 conclusions and consequences are given. A further paper seeks to extend this work to non-parallel shocks (Meli and Quenby, 2002b).

## 2 Numerical Method

The actual purpose of these simulations is to find a solution to the particle transport equation for highly relativistic parallel flow velocities. The appropriate time independent Boltzmann equation is given by the following as we assume steady state,

$$\mathrm{\Gamma }(V+\upsilon \mu )\frac{f}{x}=\frac{f}{t}|_c$$
(1)
where $`\mathrm{\Gamma }`$ is the Lorentz factor of the fluid frame, $`V`$ is the fluid velocity, $`\upsilon `$ the velocity of the particle, $`\mu `$ the cosine of the particle’s pitch angle and $`\frac{f}{t}|_c`$ the collision operator. The reference frames we use during the simulations are the upstream and downstream fluid frames and the shock frame. Pitch angle is measured in the local fluid frame, while $`x`$ is in the shock rest frame and the model assumes variability in only one spatial dimension. In order for the above equation to be solved by a simple Monte Carlo technique, we assume that the collisions represent scattering in pitch angle with no cross-field diffusion, the scattering is elastic in the fluid frame and a phase averaged distribution function may be employed. We will discuss these assumptions in more detail in due course. The collision operator may be characterized by two types of scattering that the particles suffer, large angle and small angle scattering, which are respectively represented by the two following equations:

$$\frac{f}{t}|_c=\frac{<f>_\mu f}{\tau }$$
(2)
where $`\tau `$ is the mean time between collisions, related to the mean free path by the relation $`\lambda =\upsilon \tau `$ and,

$$\frac{f}{t}|_c=\frac{}{\mu }(D_{\mu \mu }\frac{f}{\mu })$$
(3)
where $`D_{\mu \mu }`$ is the pitch angle diffusion coefficient. As mentioned previously, Gallant and Achterberg (1999) have demonstrated that small angle scattering with $`\delta \theta <1/\mathrm{\Gamma }`$ applies with $`\theta `$ measured in the upstream fluid frame for scattering in a uniform field or a randomly orientated set of uniform field cells. This arises because particles attempting to penetrate upstream from the shock are swept back into the shock before they can move far in pitch. To establish the conditions for large angle scattering, we demand the particle penetrating upstream can move over about a Larmor radius, $`r_g`$, on encountering the scattering centre while still upstream. A suitable scattering centre would be a rotational discontinuity or a much enhanced field region. Let $`\upsilon `$ be the particle velocity in the shock frame, directed perpendicular to the shock in the upstream direction, the most favourable case for re-entry upstream. If $`\upsilon c`$, transforming to the upstream frame via $`\upsilon _1=(\upsilon V_1)/(1\frac{\upsilon V_1}{c^2})`$ and allowing $`V_1=c\delta c`$ yields $`\upsilon c`$. Now in the upstream frame, the shock velocity is $`c+\delta `$, if the positive spatial coordinate, $`x`$, is towards the shock. If $`\lambda _1`$ is the mean free path for scattering (encountering a large angle scatterer) in the upstream frame, the average time to reach this scatter centre after leaving the shock is $`\mathrm{\Delta }t_1=\lambda _1/c`$. Meanwhile, the shock moves in the upstream frame a distance $`(c\delta )\mathrm{\Delta }t_1`$. Hence we require $`r_{g,1}<\lambda _1(c\delta )\mathrm{\Delta }t_1`$ to allow space for the upstream large angle scatter to occur. If $`\mathrm{\Gamma }_1=[1(V_1^2/c^2)]^{0.5}`$, we find $`r_g<\lambda _1/2\mathrm{\Gamma }_1^2`$ as the condition for large angle scattering.
The most favourable site to find high $`\mathrm{\Gamma }`$ large angle scattering is in GRB which may arise from collapse resulting in a rapidly rotating black-hole torus system (van Putten, 2001). Rotation allows a $`1/r`$ field and we start from a torus of about one Schwarzchild radius, depending on spin, of a $`20M/M_{}`$ hypernova collapse and the presence of a field near the torus of about $`B_{}=10^{14}`$ gauss. Although Brainerd (1992) shows that for typical GRB parameters and a cosmological jet model, cyclotron and synchrotron radiation are only possible for fields below about $`10^{10}`$ gauss, his arguments apply only to the regions from which emission emerges. Lieu et al. (1989) are among authors pointing out the extra quantum mechanical processes arising at $`\mathrm{\Gamma }B>10^{14}`$, thus suggesting that the maximum possible field strength, which could arise in equipartition with the burst binding energy, $`10^{18}`$ gauss, will rapidly reduce in strength as discussed by Lerche and Schramm (1977). Blobs, rather than a steady ’wind’ arise in some models (e.g. Blackman et al., 1996). Meszaros et al. (1993) suggest a comoving fireball size during the main acceleration phase $`r_{b,1}10^{14}10^{16}`$ cm and a mean fireball field just due to equipartition with internal shocks of $`B_{a,1}10^2`$ gauss. The suffix ’1’ is used for comoving quantities and this suffix is omitted for fixed frame quantities. The transverse comoving frame ’bullet’ fields which are likely to be effective for our scattering are estimated by $`B_s=B_{}r_{}/r`$ where $`B_1=B/\mathrm{\Gamma }`$ and $`r_1=\mathrm{\Gamma }r`$. Let $`\lambda _{a,1}`$ be the parallel mean free path in the ambient or general shocked field outside the ’bullets’ with, $`\lambda _{a,1}=nr_{g,a,1}`$ where, $`r_{g,a,1}`$ is Larmor radius and $`n>10`$ as suggested by detailed computations in the turbulent interplanetary medium (Moussas et al., 1992 and references therein). Diffusive shock acceleration can hardly be expected unless the mean free path is significantly smaller than the dimensions of the region. Also $`r_g\gamma `$, the particle Lorentz factor. Hence we may write $`\lambda _{a,1}=min[r_{b,1}/10,n\rho _{1,a}]`$. Then inserting this expression in our large angle scatter limit requires either $`B_{a,1}/B_{s,1}<n/2\mathrm{\Gamma }^2`$ independent of $`\gamma `$ or,

$$p_1(eV)<\frac{300B_{s,b,1}r_{b,1}}{10\times 2\mathrm{\Gamma }^2}(gauss,cm)$$
(4)
where, the particle momentum $`p_1`$ and field $`B_{s,b,1}`$ are measured near the fireball boundary. Taking the extreme $`\mathrm{\Gamma }=10^3`$, with $`n`$ = 40, the first inequality is marginally satisfied if $`r_{b,1}10^{14}`$ cm while the second inequality is satisfied at a moving frame momentum of $`<10^{16}`$ eV or rest frame momentum of $`<10^{19}`$ eV. Thus, it is possible to envisage GRB models with $`\mathrm{\Gamma }<10^3`$ where large angle scattering is occurring, especially if a larger torus or torus field is allowed.
Here we consider one-dimensional parallel shocks ($`\psi =0^{}`$, where $`\psi `$ is the angle between the shock normal and the magnetic field), either because that is the field configuration, or because turbulence removes ’reflection’ at the interface even though the magnetic field does not lie along the $`x`$-axis. A test particle approximation is used for simplicity to allow a step by step approach to understanding the full problem, though eventually non-linear effects must be included. Relativistic particles of initial $`\gamma (\mathrm{\Gamma }+10)`$ are injected upstream towards the shock and are allowed to scatter in the respective fluid frames with their basic motion described by the guiding centre approximation. One justification is we note that Bell (1987) and Jones and Ellison (1991) have shown that ’thin’ sub-shocks appear even in the non-linear regime, so at some energy above the plasma value, the accelerated particles may be dynamically unimportant while they recross the discontinuity. Another way of arriving at the test particle regime is if particles are injected well above the plasma particle energy, remaining dynamically unimportant and thus requiring the seed particles to have been already pre-accelerated. In AGN jets, traveling shocks superimposed on the relativistic flows accelerated by pressure at the jet base, could provide the seed for the terminal ’hot spot’ acceleration. Such seed particles appear in the neutron star binary merger scenario for GRB (Narayan et al., 1992), which includes the presence of pre-accelerated particles before any terminal shock acceleration phase, possibly due to explosive reconnection. For isolated pulsars, we note that there is strong observational evidence (e.g. PSR1913+16, PSR2127+11C) for the continuous injection of relativistic particles into the surrounding medium, over their lifetime. While this plasma will probably be predominantly in the form of $`e^+e^{}`$ pairs, created in the pulsar magnetosphere, it has been argued that pulsar winds must also contain ions in order to account for the electric currents in the Crab Nebula (Hoshino et al., 1992; Gallant and Arons, 1994). These authors employ a wave acceleration process at the shocked termination of the flow, to provide a non-thermal pair spectrum to yield the X-ray and gamma-ray output. However their input distribution function is also not mono-energetic and therefore containing ’seed particles’ above the bulk flow energy. The supposed nature of the slot/gap electric field acceleration of the pair plasma is unlikely to give a mono-energetic output.
A relativistic transformation is performed from the local plasma frames to the shock frame to check for shock crossings. In this parallel shock case, there is no reflection allowed and on crossing the computation is continued in the other plasma frame. We change units so that $`mc=1`$. The particle is made to leave the system if it ’escapes’ far downstream at the spatial boundary, or if it reaches a well defined maximum energy $`E_{max}=E_o10^{14}`$, for computational convenience, even though other physics describing particle escape or energy loss would probably need to be taken into account in realistic situations. The downstream spatial boundary required can be estimated from the solution of the convection-diffusion equation in a non-relativistic, large-angle scattering approximation in the downstream plasma which gives the chance of return to the shock, $`exp(V_2r_b/D)`$, yielding a probability of return of $`210^3`$ if $`r_b=8\lambda _{}`$. In fact, runs are performed with different spatial boundaries to investigate the effect of the size of the acceleration region on the spectrum, as well as to find a region where the spectrum is size independent. In the small-angle scatter case, the inherent anisotropy due to the high downstream sweeping effect may greatly modify this analytical estimate.
The particles ($``$$`10^5`$) of a weight ($`w_p`$) equal to 1.0, are injected far upstream and they are allowed to move towards the shock, along the way colliding with the scattering centers. Provided multiple scattering between the upstream and downstream regions of the shock can occur and they gain energy in each crossing, diffusive shock acceleration is simulated. The compression ratio ($`r`$) is kept at the value of 4 for simplicity.

A splitting technique is applied, similar to Bednarz and Ostrowski’s (1998) in order to obtain statistical accuracy over a wide range of particle Lorentz factors. The splitting technique helps to avoid the consequence that in highly relativistic flow environments, only a few high energy accelerated particles, which remain in the acceleration process, dominate the recorded distribution function, thereby limiting the statistical accuracy to an energy range of only two or three decades above injection energy $`E_o`$. By introducing more particles, but of corresponding lower weight when some energy well above injection is attained, a greater sampling of the possible regions within the accessible phase space of the model is achieved.
The mean free path $`\lambda `$ for particle scattering is calculated in the respective fluid frames and if it is assumed that it is dependent on the particle’s energy as given by the formula, $`\lambda =\lambda _opcos\theta `$, where $`\theta `$ is the particle’s pitch angle, so as to allow for a particle ’current’ and phase space weighting, the probability that a particle will move a distance $`x`$ along the field lines at pitch angle $`\theta `$ before a scattering is given by the expression, $`p(x)exp(x/\lambda )`$. At the scattering centers the energy (momentum) of the particle is kept constant and only the direction of the velocity vector $`\upsilon `$ is randomized, using a computational random number generator. The downstream mean free path is taken as a factor 4 less than upstream while the absolute value is arbitrary since the model boundaries are specified in units of the mean free path.
To model the large angle scatter case, the new pitch angle is chosen at random (phase angle is neglected in this guiding centre approximation but will need to be chosen to specify the scattering). For the picture of pitch angle diffusion case, we assume that the tip of a particle’s momentum vector undergoes randomly a small change $`\delta \theta `$ in its direction on the surface of a sphere and within a small range of polar angle (after a small increment of time) (Ellison et al., 1990). If the particle had an initial pitch angle $`\theta `$, we calculate its new pitch angle $`\theta ^{}`$ by the simple trigonometric formula,

$$cos\theta ^{}=cos\theta \sqrt{1sin^2\delta \theta }+sin\delta \theta \sqrt{1cos^2\theta }cos\varphi $$
(5)
where, $`\varphi (0,2\pi )`$ is the azimuth angle with respect to the original momentum direction and must be chosen at random in this case. For the case of highly relativistic flows, we limit particles pitch angle diffusion to angles chosen at random, up to an angle $`1/\mathrm{\Gamma }`$, where $`\mathrm{\Gamma }`$ is the upstream gamma, measured in the shock frame (Gallant and Achterberg, 1999) in order to be consistent with these authors’ model for upstream re-entry but some runs are performed with the scattering angle $`0.1/\mathrm{\Gamma }`$ to correspond to the later work of Achterberg et al. (2001) and Protheroe et al. (2002) -private communication.
A consequence of this limitation in upstream scatter for the Gallant and Achterberg (1999) model follows if we calculate the ratio of energies crossing down-upstream measured in the respective fluid frames as $`\mathrm{\Gamma }(1+\beta _r\mu _u^{^{}})`$ where $`\beta _r`$ is the relative velocity of the two streams as a fraction of $`c`$ and $`\mu _u^{^{}}`$ is cosine pitch angle for down to upstream crossing in the downstream frame. This angle only needs to exceed $`1/41/3`$ for kinematic reasons. Hence the energy gain can still be $`\mathrm{\Gamma }`$ in contrast to the gain ratio for the up to down transition which is $`\mathrm{\Gamma }(1\beta _r\mu _d)`$ and where, $`\mu _d`$ the upstream frame cosine pitch angle on up to down transmission is nearly unity due to the limitation on the scattering permissible.
This prediction of limited energy gain on all crossings subsequent to the first cycle, will be investigated in our simulations, as it is important to establish the dependence of the model results on whether or not large angle or pitch angle diffusion operate in highly relativistic flows.

## 3 Results

First, test runs to verify the validity of the codes which have been performed in the mildly relativistic limit (e.g. $`V_1=0.1c0.6c`$) tend to excellent agreement with the non-relativistic theory, giving smooth spectra shapes and spectral indexes close to -2.2. In addition, similar runs show that as the flow becomes more relativistic, the spectrum flattens, which is in agreement with the analytical and numerical work of Kirk and Schneider (1987a,b). We also compare our results with similar studies for $`\mathrm{\Gamma }5`$ (Baring, 1999) and find agreement concerning the spectral index and the spectral shape for both large angle scattering and pitch angle diffusion mechanisms.
As we have already mentioned a compression ratio ($`r`$) of 4 is used for simplicity and also to allow immediate comparison with non-relativistic unmodified strong shock results, but our preliminary code runs show that the qualitative trend of the results is insensitive to the exact compression ratio value, due to the mildly-relativistic nature of the downstream plasma. Initial Lorentz factors of the flow investigated are 5, 50, 200, 500 and 990.

All results are given at the shock, in the shock frame. The particle distribution function is obtained in a particular energy and $`\mu `$ space cell by recording the passage of each simulation particle crossing the cell, weighted by the time taken in cell crossing.
In figure 1 we show the ratio of the computational time constant to the non-relativistic analytical acceleration time constant, as defined in the ’Introduction’ section, calculated as a function of the upstream Lorentz factor $`\mathrm{\Gamma }`$ flow for large and small angle scattering. Note that we are effectively measuring ’speed-up’ of the acceleration time in units of a few times $`r_g/c`$.

We obtain the time constant $`t_{acc}`$ from the computation by recording the energy increment for each complete cycle, up-down-up and the time taken for the cycle at a particular energy level, $`\mathrm{\Gamma }`$ and finding the mean increments, $`\mathrm{\Delta }\mathrm{\Gamma }`$ and $`\mathrm{\Delta }t`$. Then $`t_{acc}/\mathrm{\Delta }t=\mathrm{\Gamma }/\mathrm{\Delta }\mathrm{\Gamma }`$. We find a ’speed-up’ of about a factor of 5 for large angle scattering and a considerable speed-up of a factor $``$ 20 for pitch angle diffusion. In this latter case, the parallel mean free path is estimated for the analytical formula from the time to multiple scatter in pitch angle through $`\pi /2`$. That is $`\lambda _{}=(\pi ^2/4)(c\lambda /\delta \theta ^2)`$ where $`\lambda \gamma `$ is the step length between collisions yielding $`\delta \theta `$.

Particle energy is specified in units of $`\gamma `$ and hence applies both to electron and ion acceleration. Figures 2 and 3 show the logarithm of the particle’s mean-energy against the number of shock crossings for the two scattering models, at $`\mathrm{\Gamma }`$=50 and 990. Here crossings 1 to 3 (downstream $``$ upstream$``$ downstream) represent one cycle.
The energy is measured immediately after the particle has crossed the shock front in the shock frame. In all cases, the energy gain on the first complete cycle, specified by crossing numbers 1 and 3 is $`\mathrm{\Gamma }^2`$ but subsequent crossings show a reduced gain, passing through a region where the increase is ($`\mathrm{\Gamma }`$) and tending towards the factor $`2`$ predicted by Gallant and Achterberg (1999).

However, the gain reduction is slower than expected by these previous authors and the qualitative idea expressed above. The discrepancy could lie in details of the scattering model in that the computation assumes negligible scattering before the first scatter centre is reached, equivalent to almost undeviated guiding centre motion along a field line, whereas the Gallant and Achterberg (1999) approach takes the opposite extreme of continued deflection in a plane at a large angle to a uniform field.
In figure 4 the angular distribution is measured in the shock frame just downstream for the upstream to downstream transition and just upstream for the down to up transition, for large angle scattering and upstream flow $`\mathrm{\Gamma }=200`$. Similar plots for small angle scattering (pitch angle diffusion case) are also presented. We only show the contribution to the distribution function of particles which have just crossed, not the complete time averaged distribution.
Similar results were obtained for $`\mathrm{\Gamma }=10990`$. It is evident that there is extreme peaking in the angular distribution towards smaller pitch angles in all but one case, continuing the trend noticed by Quenby and Lieu (1989) and Ellison et al. (1990) at lower flow $`\mathrm{\Gamma }`$, who found that as the velocities become more relativistic the anisotropies in pitch angle become greater. If the distribution function at the shock is isotropic, the number of particles between $`\mu `$ and $`\mu +\mathrm{\Delta }\mu `$ is proportional to $`\mu `$, taking into account solid angle weighting.

Both scattering cases reveal highly peaked distributions pointing downstream on up to down transition, due to the beaming of the upstream distribution function as seen by the shock.

Downstream to up for large angle scattering produces near triangular distributions in the upstream directed half space, due to the very efficient scattering in the model. For small angle scattering, highly peaked distributions are found due to the lack of time for significant deflection from nearly zero pitch angle before downstream re-entry and presumably for those particles which can re-enter upstream.

In figures 5,6 and 7 we show the spectral shapes obtained at the shock just downstream, for large and small angle scattering (pitch angle diffusion) for upstream flow $`\mathrm{\Gamma }`$ values of 5, 50 and 990 respectively. We investigated the sensitivity of these spectral shapes for small angle scattering to downstream boundary condition, and found that the shapes did not change significantly for the boundary distance, $`r_b210^3\lambda `$. Probably the small chance of return to the shock except for a relatively small subset of downstream pitch angle particle ’histories’, produced this result and the anisotropy in figure 4 is seen for these returning particles. Note the relatively smooth spectra for the relativistic flow with $`\mathrm{\Gamma }=5`$ (fig. 5) which becomes plateau-like at the more extreme value of $`\mathrm{\Gamma }=990`$ (fig. 7) where the effects of individual acceleration cycles are clearly evident in the large angle case, but not for small angle scattering, where the spectral shape remains rather ’smooth’. Features can be seen developing in the lower $`\mathrm{\Gamma }`$ simulations of Quenby and Lieu (1989) and Ellison et al. (1990) while Protheroe (2001) shows a similar contrast in behaviour, between large and small angle scattering models up to $`\mathrm{\Gamma }=20`$. On any individual cycle, greater energy gain is available for the large angle case because of the greater upstream scattering, as demonstrated in the simple analytical approach mentioned above.
Hence, the plateaus develop for the large angle case with flatter segments, than the smoother pitch angle case spectrum. An opposing effect, however, is that because the accelerated particles re-entering downstream have pitch angles more favourable to eventual loss downstream, the drop between plateaus is greater and so the overall spectral slopes in the two cases are rather similar. The actual mean slope of the differential number spectrum is well above -2, as illustrated by the lines of slope -1.45 and -1.52 shown for $`\mathrm{\Gamma }`$=5. This relativistic flattening effect is consistent with previous studies by Baring (1999) while, the mechanism of plateau development as an acceleration cycle effect is implicit in figure 6 of Protheroe (2001). It is perhaps not surprising to observe a step-like behaviour in the accelerated spectrum, clearly related to groups of particles undergoing different numbers of acceleration cycles, because the simulations also reveal that about $`90\%`$ of all particles in the simulation box are lost downstream each cycle. One might also expect that small angle scattering has the greater chance of smoothing the inherent structure in the acceleration process. Clearly, smoothing is much more efficient at lower flow $`\mathrm{\Gamma }`$ where the gain per cycle is less.
One can infer that at extreme relativistic flow velocities and different scatter models, it is difficult to assume that the spectral shape of the accelerated particles follows a smooth power-law shape, as in the classic non-relativistic first order Fermi acceleration mechanism, under all circumstances.

## 4 Conclusions

To what extent is first order Fermi acceleration universally applicable in shocked, turbulent, MHD flows ? Aspects of the problem still needing elucidation include the spectrum produced, the spectral index, the angular distribution of the particles and non-linear effects in relativistic flows.
We have presented a detailed Monte Carlo approach for the highly relativistic diffusive acceleration in parallel shocks under the test particle approximation, extending previous studies to more extreme values of the upstream flow $`\mathrm{\Gamma }`$ and carefully distinguishing the effects of different assumptions concerning the particle scattering. Comparison is made with results under non-relativistic conditions. It is found that the spectral shape in the very relativistic flow regime depends on the scatter model (large angle or pitch angle diffusion/small angle scattering), with the former model exhibiting a structure in the spectrum at the shock related to the cycle of acceleration, while the later one showing a smoother shape. One general indication which supports the validity of our computational approach is that the spectral index becomes flatter as the flow becomes more relativistic, a result to be compared with the conclusions of Kirk and Schneider (1987a,b).
Extreme anisotropy in the shock frame for particles emerging downstream irrespective of scattering model is observed, due to the beaming of the transmitted flux. High anisotropy also occurs in the flux emerging upstream in the small angle scattering case, due to the lack of time to isotropize the particles but in the large angle scattering model there is a large enough probability of scatter before shock re-entry for near isotropy to set in. This means that for any given acceleration cycle, the spectrum is much flatter in the large angle case (plateaus) because the gain per cycle can be greater. However, evidently because of an enhanced loss downstream of the more scattered particles of the large angle scattering case, the average spectral slope for the two mechanisms is similar.
Most important, a significant acceleration rate increase (’speed up’ effect) has been found for both large angle and pitch angle diffusion. This must be seen in conjunction with a ’gamma-squared’ energy boosting of the particle noted in the first shock cycle, which is in accordance with theoretical indications (e.g. Gallant and Achterberg, 1999) and past simulations (e.g. Lieu and Quenby, 1989). There is a gain by a significant factor in subsequent crossings which is $`\mathrm{\Gamma }`$ or even more for several cycles. Regarding the acceleration time decrease, computations show that the comparison of the simulation (experimental) acceleration time with the analytic non-relativistic time constant, indicates a decrease of about 5 with increasing flow gammas (10-$`10^3`$). This speed-up is clearly directly related to the boosting of the energy gain.
In the theory of GRB production (Meszaros and Rees, 1993) by cosmic fireballs, the $`\mathrm{\Gamma }10^3`$ flows are expected to produce some shock acceleration of electrons and protons. Vietri (1995) suggests protons of $`\gamma 10^{11}`$ produce photons of 300 GeV by proton synchrotron radiation and neutrinos of $`10^{19}`$ eV via photopion production and subsequent cascading. This requires only two shock acceleration cycles at flow $`\mathrm{\Gamma }10^3`$, a suggestion well supported by the results obtained here for both scattering models. However, at these high $`\mathrm{\Gamma }`$ values, extreme ’billard ball’ scattering is required over distances much less than $`\lambda `$ in the large angle case. What finally limits the proton maximum energy is escape from the shock turbulence region, as the Larmor radius of the accelerated particles increases in size. Despite the more favourable radiative process, electrons are not able to produce gamma-ray synchrotron energies higher than those due to these protons because the competing synchrotron energy loss severely limits the diffusive shock gain (e.g. Gallant et al., 2000) so, the net result is a similar upper photon cutoff. Because the accelerated particle spectra at the shock as here are not necessarily smooth, it might be expected that a detailed experimental study of radiation from this region of a GRB may reveal a related structure. In turn, such a study may reveal new insight into the magnetic field regime of the GRB environment.
The suggestion by Vietri (1995) that the efficient production in GRB of $`10^{20}`$ eV protons constitutes a cosmic ray source, at least in previous epochs, is also supported. This argument depends on the similarity of the burst gamma ray energy output and very high energy cosmic ray energy flux (Vietri, 1995; Waxman, 1995).
Finally, we may conclude that first order Fermi acceleration in the sense of additive energy change over a number of cycles of magnetic discontinuity crossing, appears to be applicable to $`\mathrm{\Gamma }10^3`$ flows where, there is no test particle reflection at the discontinuity. However, the large gradients in particle distribution function appearing in the simulations argue against a spatial diffusion approximation description of the phenomenon.

Acknowledgments

We wish to thank Prof. Drury and the unknown referee for their valuable comments and suggestions regarding this work.