text:
# Chandra observations of RXJ1347.5-1145: the distribution of mass in the most X-ray luminous galaxy cluster known

## 1 Introduction

The launch of the Chandra Observatory (Weisskopf et al. 2000) in 1999 July has provided the first opportunity for detailed, spatially-resolved X-ray spectroscopy of clusters of galaxies at moderate to high redshifts. The Advanced CCD Imaging Spectrometer (ACIS) on Chandra permits direct, simultaneous measurements of the X-ray gas temperature and density profiles in clusters and, via the hydrostatic assumption, the mass distributions, spanning scales from $`r10`$ kpc in cluster cores out to $`r1`$ Mpc at the detection limits.

In this paper, we present the first results from Chandra observations of RXJ1347.5-1145, the most X-ray luminous galaxy cluster known (Schindler et al. 1995). This cluster has been the subject of several previous X-ray, optical and Sunyaev-Zeldovich (SZ) effect studies (e.g. Schindler et al. 1995, 1997; Fischer & Tyson 1997, Allen 1998; Komatsu et al. 1999, 2001; Pointecouteau et al. 1999, 2001). Here, we present the first measurements of the X-ray temperature structure within the cluster and place tight constraints on the total mass distribution. We compare our results with those from strong and weak lensing analyses. A comparison of the Chandra and SZ results, and a determination of the Hubble Constant from the combined data sets, is presented elsewhere (Schmidt & Allen 2002, in preparation).

Except where stated otherwise, the cosmological parameters $`H_0`$=50 $`\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, $`\mathrm{\Omega }=1`$ and $`\mathrm{\Lambda }=0`$ are assumed. At the redshift of RXJ1347.5-1145 ($`z=0.451`$), an angular scale of 1 arcsec corresponds to physical size of 6.81 kpc in this cosmology.

## 2 Observations

The Chandra observations of RXJ1347.5-1145 were carried out using the ACIS on 2000 March 05 and 2000 April 29. The target was observed in the back-illuminated S3 detector and positioned near the centre of node-1 on CCD 7. The light curves for both observations were of high quality with no strong background flares. The net good exposure times were 8.9 and 10.0 ks, respectively, giving a total good exposure time of 18.9 ks. The focal plane temperature at the time of both observations was -120C.

We have used the CIAO software and the level-2 events files provided by the standard Chandra pipeline processing for our analysis. Only those X-ray events with grade classifications of 0,2,3,4 and 6 were included in our final cleaned data sets.

## 3 X-ray imaging

### 3.1 X-ray morphology

The raw $`0.37.0`$ keV image of the central $`6\times 6`$ arcmin2 region of the cluster, from the combined 18.9 ks data set, is shown in Fig. 1(a). The pixel size is $`1.97\times 1.97`$ arcsec2, corresponding to $`4\times 4`$ raw detector pixels. Fig. 1(b) shows an adaptively smoothed contour plot of the same data, using the smoothing algorithm of Ebeling, White & Rangarajan (2002). The images reveal a number of notable features. Firstly, we see a sharp central surface brightness peak at a position 13h47m30.63 -11d45m09.3s (J2000.), in good agreement with the optical centroid for the dominant cluster galaxy, 13h47m30.5 -11d45m09s (J2000.; Schindler et al. 1995). Secondly, we identify a region of enhanced emission $`20`$ arcsec to the southeast of the X-ray peak (see below). Thirdly, on large scales ($`r80`$ arcsec) we detect an extension of the X-ray isophotes to the south.

Fig. 2 shows the data for the central regions of the cluster on a finer $`0.492\times 0.492`$ arcsec2 scale (corresponding to $`1\times 1`$ raw detector pixels). From this figure it appears that the region of enhanced emission to the southeast of the X-ray peak has a roughly circular morphology with a flattened northwest edge. Beyond this edge, at a radius $`r10`$ arcsec, the X-ray subclump appears to be separated from the main cluster core by a valley of reduced emission. The spectral analysis in Section 4.2 shows that the X-ray subclump has a hotter temperature than the gas at the same radius in other directions from the X-ray peak and probably contains shocked gas resulting from a recent subcluster merger event.

Komatsu et al. (2001) report the detection of a region of enhanced SZ effect at a position coincident with the X-ray subclump. The enhanced SZ effect is consistent with the higher temperature and density for this region measured from the Chandra data.

No strong X-ray point sources are detected in the Chandra S3 field. (The brightest point source has a flux of $`3.1\times 10^{14}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$in the $`0.57.0`$ keV band; Gandhi et al. , in preparation.)

### 3.2 Surface brightness profiles

The azimuthally-averaged, $`0.37.0`$ keV X-ray surface brightness profile for RXJ1347.5-1145 for position angles of $`18090`$ degrees (i.e. excluding the southeast quadrant) and for the disturbed, southeast quadrant are shown in Figs. 3(a,b), respectively. The profiles have been flat-fielded and background subtracted using an on-chip region free from cluster emission. The bin-size is 2 detector pixels (0.984 arcsec).

We see that the data for position angles of $`18090`$ degrees appear regular and, for $`r<500`$ kpc, can be approximated ($`\chi ^2=102`$ for 72 degrees of freedom) by a $`\beta `$-model (e.g. Jones & Forman 1984) of the form $`S_\mathrm{X}(r)=S(0)\left[1+(r/r_\mathrm{c})^2\right]^{1/23\beta }`$ with a core radius $`r_\mathrm{c}=29.2\pm 0.7`$ kpc and a slope parameter $`\beta =0.535\pm 0.003`$ ($`1\sigma `$ errors). Between radii of $`0.11.0`$ Mpc we find evidence for a steepening of the surface brightness profile with increasing radius: for $`0.1<r<1.0`$ Mpc the data can be described ($`\chi ^2=135`$ for 131 degrees of freedom) by a simple broken power-law model, with a break at a radius of $`487_{16}^{+23}`$ kpc, and slopes in the regions internal and external to the break radius of $`2.25\pm 0.03`$ and $`3.66\pm 0.18`$, respectively. A similar result was obtained for Abell 2390 by Allen, Ettori & Fabian (2001b), although in that case the slope interior to the break radius was flatter.

The data for the southeast quadrant containing the X-ray subclump reveal a clear excess of emission relative to other directions. Enhanced emission is detected between radii of $`70300`$ kpc, with a maximum enhancement of a factor $`3`$ at a radius $`r180`$ kpc.

## 4 Spatially-resolved spectroscopy

### 4.1 Method

For our spectral analysis, we have separated the southeast quadrant containing the X-ray subclump from the rest of the cluster. The data for position angles of $`18090`$ (the relaxed part of the cluster) were divided into the annular regions detailed in Table 1.111For the southeast quadrant, we have extracted a single spectrum covering radii $`60195`$ kpc (18-58 pixels), where the emission enhancement in Fig. 2 is most obvious.  A single spectrum was extracted for each region in 1024 Pulse Invariant (PI) channels. The spectra were re-grouped to contain a minimum of 20 counts per PI channel thereby allowing $`\chi ^2`$ statistics to be used.

Background spectra, appropriate for the regions studied, were extracted from the ACIS-S3 blank-field data sets provided by Maxim Markevitch (which are available from the Chandra X-ray Center). Separate photon-weighted response matrices and effective area files were constructed for each region using the available calibration and response files appropriate for the -120C focal plane temperature.222For each $`32\times 32`$ pixel2 sub-region of the S3 chip, a spectral response (.rmf) and an auxiliary response (.arf) matrix were created using the CIAO tools $`mkrmf`$ and $`mkarf`$, respectively. For each of the annular regions studied, the number of source counts in each $`32\times 32`$ pixel2 sub-region was determined. The individual .rmf and .arf files were then combined using the FTOOLS programs $`addrmf`$ and $`addarf`$ (using a script provided by Roderick Johnstone) to form a counts-weighted spectral response and auxiliary response matrix appropriate for each annulus.

The data from both Chandra observations were modelled simultaneously using the XSPEC code (version 11.01; Arnaud 1996). We have limited our analysis to the $`0.57.0`$ keV energy band, over which the calibration of the back-illuminated CCD detectors is currently best understood.

The spectra have been modelled using the MEKAL plasma emission code of Kaastra & Mewe (1993; incorporating the Fe L calculations of Liedhal, Osterheld & Goldstein 1995) and the photoelectric absorption models of Balucinska-Church & McCammon (1992). We have fitted each annular spectrum using a simple, single-temperature model with the absorbing column density fixed at the nominal Galactic value ($`N_\mathrm{H}=4.85\times 10^{20}`$atom cm-2; Dickey & Lockman 1990). The free parameters in this model were the temperature ($`kT`$), metallicity ($`Z`$; measured relative to the solar photospheric values of Anders & Grevesse 1989, with the various elements assumed to be present in their solar ratios) and emission measure. We note that including the absorbing column density as an additional free parameter did not result in significant improvements in the goodness of fit.

### 4.2 Results

The best-fit parameter values and $`1\sigma `$ ($`\mathrm{\Delta }\chi ^2=1.0`$) confidence limits determined from the fits to the annular spectra (excluding the southeast quadrant) in the $`0.57.0`$ keV band are summarized in Table 1. The projected temperature profile determined with this model is shown in Fig. 4. The temperature rises from a mean value of $`7.0_{0.6}^{+0.8}`$ keV within $`r=24`$ kpc to $`kT=15.8_{1.2}^{+1.6}`$ keV over the $`0.10.5`$ Mpc region.

A fit with the same model to the data for the southeast quadrant between radii of $`60195`$ kpc (where the X-ray enhancement is most obvious in the image) gives a best-fit temperature $`kT=18.0_{2.3}^{+2.7}`$ keV. (A similar result, $`kT=18.2_{2.9}^{+3.8}`$ keV, is obtained using a circular region of radius 10 arcsec centred on the enhancement.) In other directions, the mean temperature between radii of $`60195`$ kpc is $`kT=12.7\pm 1`$ keV. The higher temperature and enhanced X-ray surface brightness in the southeast quadrant are consistent with the X-ray gas in that region having undergone shock compression (Section 8).

We find marginal evidence for a metallicity gradient in the cluster, with a mean metallicity for the central 100 kpc radius of $`Z=0.55\pm 0.09`$ solar, as opposed to $`Z=0.34\pm 0.12`$ solar for the $`100740`$ kpc region.

### 4.3 Spectral deprojection analysis

The results discussed above are based on the analysis of projected spectra. We have also carried out a simple deprojection analysis of the Chandra spectral data using the method described by Allen et al. (2001b).

For this analysis we have used the same annular regions (covering position angles of $`18090`$ degrees) and have assumed that the emission from each spherical shell (the shells are defined by the same inner and outer radii as the annular regions) is isothermal and absorbed by the Galactic column density. The fit to the outermost annulus is used to determine the temperature and emission measure in the outermost spherical shell. The contribution from that shell to each inner annulus is then determined by purely geometric factors (e.g. Kriss, Cioffi & Canizares 1983). The fit to the second annulus inward is used to determine the parameters for the second spherical shell, and so forth, working inwards.

The data for all eight annular spectra were fitted simultaneously in order to determine the parameter values and confidence limits. The metallicity was linked to take the same value at all radii. The temperature profile determined with the spectral deprojection method is shown in Fig. 5.

### 4.4 Comparison with previous work

The mean emission-weighted temperature and metallicity for RXJ1347.5-1145 of $`kT=12.2\pm 0.6`$ keV and $`Z=0.42\pm 0.07`$ solar, respectively, measured in the $`0.00.74`$ Mpc range over the full 360 degrees, are in good agreement with the values reported by Allen & Fabian (1998; $`kT=12.5_{0.8}^{+0.9}`$ keV, $`Z=0.38_{0.10}^{+0.11}`$ solar) using ASCA data and a similar spectral model. Excluding the data for the southeast quadrant, we obtain only small changes in these results: $`kT=12.0\pm 0.6`$ keV and $`Z=0.41\pm 0.07`$ solar. The Chandra results are also in good agreement with those reported by Schindler et al. (1997) using the same ASCA data and a 3 arcmin radius aperture ($`kT=11.8_{1.0}^{+1.6}`$ keV, $`Z=0.39_{0.13}^{+0.12}`$ solar). The emission-weighted temperature measured with Chandra is slightly cooler than the value of $`14.5_{1.5}^{+1.7}`$ keV measured with BeppoSAX by Ettori, Allen & Fabian (2001), as can be expected given the observed temperature profile (Figs. 4,5) and the different response characteristics of the instruments. (The data for the southeast quadrant are not excluded from the BeppoSAX or ASCA studies.)

The ‘ambient’ emission-weighted temperature measured with Chandra in the $`0.10.5`$ Mpc range (i.e. excluding the data for the southeast quadrant, the central, cool region, and the regions at large radii where the temperature may start to decline again) of $`kT=15.8_{1.2}^{+1.6}`$ keV is in good agreement with the value of $`15.9_{2.6}^{+6.5}`$ keV estimated from BeppoSAX data by Ettori et al. (2001) using a spectral model which included a constant-pressure cooling flow component. The Chandra result is also consistent with the ambient temperature of $`kT=26_{12}^{+8}`$ keV estimated from ASCA data using a similar constant-pressure cooling flow model (Allen & Fabian 1998), and $`kT=18.6_{3.0}^{+4.1}`$ keV from a re-analysis of the ASCA data by the same authors using the more appropriate ‘isothermal’ cooling flow models of Nulsen (1998).

The emission weighted ambient temperature determined from the Chandra data in the $`0.10.5`$ Mpc region is in good agreement with the mean gas mass-weighted temperature of $`16.1_{2.7}^{+5.3}`$ keV measured within $`r_{2500}`$ ($`r=0.72`$ Mpc) by Allen, Schmidt & Fabian (2001c). We note that this temperature is also consistent with the predicted value of $`kT=16.8`$ keV, using the cooling-flow corrected $`kT_\mathrm{X}/L_{\mathrm{Bol}}`$ relation of Allen & Fabian (1998; we assume a bolometric luminosity of $`2.2\times 10^{46}`$ $`\mathrm{erg}\mathrm{s}^1`$ as measured by ASCA).

## 5 Measurement of the cluster mass profile

### 5.1 Method

The observed X-ray surface brightness profile (Fig. 3a) and deprojected, spectrally-determined temperature profile (Fig. 5) may together be used to determine the X-ray gas mass and total mass profiles in the cluster. For this analysis, we have used an enhanced version of the image deprojection code described by White, Jones & Forman (1997), and have followed the methods outlined by Allen et al. (2001b) and Schmidt et al. (2001).

A variety of simple parameterizations for the cluster mass distribution were examined to establish which could provide an adequate description of the Chandra data (see below). The best-fit parameter values and confidence limits were determined by examining parameter grids and evaluating $`\chi ^2`$ for each set of parameters.333The observed surface brightness profile (Fig. 3a) and a particular parameterized mass model are together used to predict the temperature profile of the X-ray gas. (We use the median temperature profile determined from 100 Monte-Carlo simulations in which the surface brightness profile is perturbed according to the statistical uncertainties at each radius. The outermost pressure is fixed using an iterative technique which ensures a smooth pressure gradient in these regions.) The predicted temperature profile is rebinned to the same binning as the spectral results and the $`\chi ^2`$ difference between the observed and predicted, deprojected temperature profiles is calculated. The parameters for the mass model are stepped through a regular grid of values in the $`r_\mathrm{s}`$-$`\sigma `$ plane (see text) to determine the best-fit values and 68 per cent confidence limits. The gas mass profile is determined to high precision at each grid point using the observed surface brightness profile and model temperature profile. Spherical symmetry and hydrostatic equilibrium are assumed throughout.

### 5.2 NFW mass models

We find that a good fit ($`\chi _{\mathrm{min}}^2=5.2`$ for 6 degrees of freedom, hereafter DOF) to the Chandra data for RXJ1347.5-1145 can be obtained using a Navarro, Frenk & White (1997, hereafter NFW) model:

$$\rho (r)=\frac{\rho _{\mathrm{crit}}(z)\delta _\mathrm{c}}{(r/r_\mathrm{s})\left(1+r/r_\mathrm{s}\right)^2},$$
(1)
where $`\rho (r)`$ is the mass density, $`\rho _{\mathrm{crit}}(z)=3H(z)^2/8\pi G`$ is the critical density for closure and

$$\delta _\mathrm{c}=\frac{200}{3}\frac{c^3}{\left[\mathrm{ln}(1+c)c/(1+c)\right]}.$$
(2)
The best-fit scale radius, $`r_\mathrm{s}=0.40_{0.12}^{+0.24}`$ Mpc and the concentration parameter, $`c=5.87_{1.44}^{+1.35}`$ (68 per cent confidence limits). The normalization of the mass profile may also be expressed in terms of an effective velocity dispersion, $`\sigma =\sqrt{50}H(z)r_\mathrm{s}c=1450_{200}^{+300}`$$`\mathrm{km}\mathrm{s}^1`$(with $`r_\mathrm{s}`$ in units of Mpc). The deprojected X-ray gas temperature profile predicted by the best-fitting NFW mass model (given the observed surface brightness profile) is shown overlaid on the observed spectral results in Fig. 5.

Fig. 6 shows a contour plot of the $`\chi ^2`$ values obtained for the full range of NFW models studied. The minimum $`\chi ^2`$ value is marked with a cross. Contours have been drawn at intervals of $`\mathrm{\Delta }\chi ^2=`$ 2.30, 6.17 and 11.8, corresponding to formal confidence limits of 68.3 (1 $`\sigma `$), 95.4 (2 $`\sigma `$) and 99.73 (3 $`\sigma `$) per cent for two interesting parameters. Fig. 7 shows the best-fitting NFW mass profile, which has a virial radius $`r_{200}=cr_\mathrm{s}=2.35_{0.33}^{+0.49}`$ Mpc and an integrated mass within the virial radius, $`M_{200}=2.29_{0.82}^{+1.74}\times 10^{15}`$$`\mathrm{M}_{}`$.

### 5.3 Other parameterized mass models

We have also examined a variety of other popular parameterized mass models. Firstly, the model of Moore et al. (1998):

$$\rho (r)=\frac{\rho _{\mathrm{crit}(\mathrm{z})}\delta _{\mathrm{c}_\mathrm{m}}}{(r/r_\mathrm{s})^{1.5}\left(1+r/r_\mathrm{s}\right)^{1.5}},$$
(3)
where

$$\delta _{\mathrm{c}_\mathrm{m}}=100\frac{c_m^3}{\mathrm{ln}(1+c_m^{1.5})}.$$
(4)
For this model we obtain $`\chi _{\mathrm{min}}^2=8.6`$ (6 DOF), with best-fit values for the scale radius $`r_\mathrm{s}=5.0`$ Mpc (the maximum allowed value in our grid), concentration parameter $`c_m=0.69`$ and effective velocity dispersion $`\sigma =2125`$$`\mathrm{km}\mathrm{s}^1`$(where $`\sigma =\sqrt{50}H_0(z)r_\mathrm{s}c_m`$).

Secondly, we examined a simple, non-singular isothermal sphere model:

$$\rho (r)=\frac{\sigma _{\mathrm{iso}}^2}{2\pi \mathrm{G}}\frac{1}{r^2+r_\mathrm{c}^2},$$
(5)
for which we obtain $`\chi _{\mathrm{min}}^2=4.9`$ (6 DOF) with $`r_\mathrm{c}=45\pm 10`$ kpc and $`\sigma _{\mathrm{iso}}=1590\pm 150`$ $`\mathrm{km}\mathrm{s}^1`$.

We thus find that all of the two-parameter mass models described above provide acceptable descriptions of the Chandra data for RXJ1347.5-1145. We note, however, that a singular isothermal sphere \[$`\rho (r)r^2`$\] does not provide an acceptable fit ($`\chi _{\mathrm{min}}^2=33.1`$ for 7 DOF).

## 6 Comparison with gravitational lensing results

### 6.1 Weak lensing

Fischer & Tyson (1997) present a detailed weak lensing study of RXJ1347.5-1145. They measure a mass profile that can be parameterized by an NFW model with $`r_{200}=2.4`$ Mpc, or a singular isothermal sphere with $`\sigma =1500`$$`\mathrm{km}\mathrm{s}^1`$, in good agreement with the Chandra results.

The azimuthally-averaged, projected surface mass density profile for RXJ1347.5-1145 determined by Fischer & Tyson (1997) is shown in Fig. 8. We have excluded the innermost point from their data, which extends into the strong lensing regime (see below; $`r_{\mathrm{max}}=2.72`$ Mpc in our adopted cosmology). The same quantity, calculated for all NFW mass models within the 68 per cent confidence contour determined from the Chandra X-ray data (Fig. 6) has been overlaid. We see that the Chandra and weak-lensing results exhibit good agreement on all spatial scales studied.

The agreement between the independent lensing and X-ray mass measurements confirms the validity of the hydrostatic assumption used in the X-ray analysis (having excluded the southeast quadrant) and suggests that the mass profile in the cluster has been robustly determined.

### 6.2 Strong lensing

Schindler et al. (1995) report the discovery of two bright gravitationally-lensed arcs in RXJ1347.5-1145 using ground based optical imaging. Sahu et al. (1998) present Hubble Space Telescope Imaging Spectrograph (STIS) observations of the cluster which reveal a number of additional strongly-lensed features. The STIS image for RXJ1347.5-1145 is shown in the left panel of Fig. 9. Sahu et al. (1998) also present ground-based spectroscopy of the bright northern arc in the cluster (arc 1) at a radius of $`34.9`$ arcsec (240 kpc) from the dominant cluster galaxy, for which they measure a redshift $`z_1=0.806`$.

Following the methods outlined in Schmidt et al. (2001), we have examined the constraints that the observed strong lensing configuration can place on the mass distribution in the cluster core. We first examined a simple, spherically-symmetric mass model, centred on the dominant cluster galaxy, with the best-fit parameter values determined from the Chandra data for the relaxed regions of the cluster (Sect. 5.2). We find that such a mass model does not provide a sufficient central mass density to explain the strong lensing data.

We next examined refined lensing models in which a second mass component, centered on the second brightest galaxy was introduced, together with appropriate ellipticities for the matter distributions. (The introduction of such a mass clump is well-motivated by the Chandra data for the southeast quadrant of the cluster. The second brightest galaxy is likely to have been the dominant galaxy of the recently-merged subcluster; Section 8). Both mass components have been modelled by elliptical NFW models (see Schmidt et al. 2001 for definitions). The parameters for the main mass component were fixed to the values determined from the Chandra data for the relaxed regions of the cluster in Sect. 5.2, i.e. scale radius $`r_\mathrm{s}=0.4`$ Mpc, velocity dispersion $`\sigma =1450`$$`\mathrm{km}\mathrm{s}^1`$ and concentration parameter $`c=5.87`$. The axis ratio and position angle were fixed to the values measured for the dominant galaxy from the STIS data: axis ratio $`q=0.76`$ and position angle $`\theta =13.1^{}`$. The centre of the main mass component was fixed to the centre of the dominant galaxy.

The parameters for the second mass component have been determined by requiring that the overall potential is able to produce the northern arc (arc 1) at the measured redshift, $`z_1=0.806`$. For an assumed scale radius $`r_{\mathrm{s},\mathrm{clump}}=0.25`$ Mpc, we obtain a best-fitting velocity dispersion $`\sigma _{\mathrm{clump}}=815`$$`\mathrm{km}\mathrm{s}^1`$ and a concentration parameter $`c_{\mathrm{clump}}=5.3`$. Such parameters are typical for mid-temperature galaxy clusters ($`kT`$ 5 keV. e.g. Allen et al. 2001c). The position, ellipticity and position angle of the second mass component were fixed to the values determined for the second brightest galaxy from the STIS image: $`q_{\mathrm{clump}}=0.79`$ and $`\theta _{\mathrm{clump}}=28.1^{}`$.

The redshift measurement for the northern arc constrains the projected mass within the arc radius. It is possible to vary the scale radius for the second mass component and to adjust the concentration parameter accordingly, so long as the projected mass within the arc radius remains constant. The absence of multiple arc images within the cluster core suggests $`r_{\mathrm{s},\mathrm{clump}}>0.1`$ Mpc.

Having adjusted the parameters for the second mass component to produce the northern arc, we then examined whether the resulting mass model can explain the other arcs observed by Sahu et al. (1998). Unfortunately, there are no measured redshifts for the other arcs. (Sahu et al. 1998 note a possible lower limit for arc 4 of $`z1`$ based on the non-detection of \[Oii\] in their ground-based spectra, although they caution that the arc may simply be a galaxy without a strong emission-line feature). The details of our model predictions are shown in the right panel of Fig. 9. The results on the predicted arc redshifts and the assumed sizes of the sources are listed in Table 2. The predicted redshifts are cosmology dependent and would be modified slightly if, for example, an accelerating universe with a different Hubble Constant were assumed. For the mass model used here, the predicted redshifts should be accurate to $`\pm 0.05`$, with the exception of arc 4 which is determined more precisely. We have modelled the lensed galaxies using circular Gaussian brightness distributions with the full-width-at half-maximum (FWHM) values listed in Table 2. The FWHMs for the arcs were adjusted so that their 30 per cent brightness contours provide a good fit to the STIS data. The effects of using different source sizes is illustrated by the different brightness contours in Fig. 9.

Although we include arc 2 in our modelling, we caution that this source appears quite red in the colour image published by Fischer & Tyson (1997) and that it is possible that the source is simply an edge-on spiral galaxy. The predicted redshift for arc 2 is lower than $`z_1`$, since otherwise the lens model would predict a much stronger distortion. Arcs 3, 4 and 5 are consistent with being at the same redshift. The projected separation of the sources producing arcs 4 and 5 is 43.4 $`h_{50}^1`$ kpc. Although arcs 4 and 5 have a similar colour in the data of Fischer & Tyson (1997), it is not possible to explain them as images of the same source using our mass model. In exploring different redshifts for arc 3, we found that for redshifts higher than $`z_31`$ it is possible to produce a counter image on the other side of the central galaxy. This alerted us to the presence of a blue object (identified from the colour image of Fischer & Tyson 1997) to the south-east of the merging subclump (see the arrow in the left panel of Fig. 9). However, both the appearance in the STIS image and our lensing analysis suggest that these objects are not images of the same source.

We note that circular NFW models for one or both of the main mass components cannot explain the observed arc shapes and orientations well. In particular, the orientation of arc 4 reflects the orientation of the critical curve (Fig. 9) in its vicinity. A small $`5^{}`$ difference between the orientation of the observed and predicted arc 4 remains in our best-fit model, which points to a small imprecision of the lensing model.

The grey circle in Fig. 8 shows the azimuthally-averaged, surface mass density contrast at $`r=240`$ kpc (centred on the main mass component) for the best-fit strong lensing model. Note that the strong lensing results should lie above the X-ray results at the same radii, since the Chandra analysis excludes the regions of the cluster affected by the second mass subclump. The surface mass density contrast rises sharply between the innermost weak lensing point and the strong lensing result, suggesting that the second mass component is centrally concentrated and does not make a major contribution to the total mass beyond the strong lensing regime.

In conclusion we find that a simple two-component mass model, with ellipticities and orientations for the mass components matching those of the dominant cluster galaxies, provides a reasonable description for the overall arc geometry in RXJ1347.5-1145. Such a model is consistent with the Chandra X-ray and weak lensing results.

## 7 The properties of the cluster gas

### 7.1 Electron density and cooling time profiles

The results on the electron density and cooling time as a function of radius, determined from the image deprojection analysis using the best-fit NFW mass model, are shown in Fig. 10. Within the central 500 kpc radius, the electron density profile can be parameterized ($`\chi ^2=43.2`$ for 34 degrees of freedom) by a $`\beta `$-profile with a core radius, $`r_c=27.8\pm 1.6`$ kpc, $`\beta =0.521\pm 0.012`$ and a central density, $`n_\mathrm{e}(0)=2.3\pm 0.1\times 10^2`$ cm-3 ($`1\sigma `$ errors). We measure a central cooling time (i.e. the mean cooling time within the central 1.97 arcsec ($`13.4`$ kpc) bin of $`t_{\mathrm{cool}}=4.1_{0.5}^{+0.6}\times 10^8`$ yr (uncertainties are the 10 and 90 percentile values from 100 Monte Carlo simulations).

### 7.2 The X-ray gas mass fraction

The X-ray gas-to-total-mass ratio as a function of radius, $`f_{\mathrm{gas}}(r)`$, determined from the Chandra data is shown in Fig. 11. We find that the best-fit $`f_{\mathrm{gas}}`$ value rises rapidly with increasing radius within the central $`r<40`$ kpc radius and then remains approximately constant, or rises slowly, out to the limits of the data at $`r=0.74`$ Mpc where we measure $`f_{\mathrm{gas}}=0.141_{0.027}^{+0.035}`$.

Following the usual arguments, which assume that the properties of clusters provide a fair sample of those of the Universe as a whole (e.g. White & Frenk 1991; White et al. 1993; White & Fabian 1995; Evrard 1997; Ettori & Fabian 1999), we may use our results on the X-ray gas mass fraction to estimate the total matter density in the Universe, $`\mathrm{\Omega }_\mathrm{m}`$. Assuming that the luminous baryonic mass in galaxies in RXJ1347-1145 is $`0.134h_{50}^{0.5}`$ of the X-ray gas mass (e.g. White et al. 1993; Fukugita, Hogan & Peebles 1998) and neglecting other possible sources of baryonic dark matter in the cluster, one can show that $`\mathrm{\Omega }_\mathrm{m}=\mathrm{\Omega }_\mathrm{b}/1.134f_{\mathrm{gas}}`$, where $`\mathrm{\Omega }_\mathrm{b}`$ is mean baryon density in the Universe. For $`\mathrm{\Omega }_\mathrm{b}h_{50}^2=0.0820\pm 0.0072`$ (O’Meara et al. 2001) we obtain $`\mathrm{\Omega }_\mathrm{m}=0.51\pm 0.12`$. For a $`\mathrm{\Lambda }`$CDM cosmology with $`H_0=70`$$`\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, our results for RXJ1347.5-1145 would imply $`\mathrm{\Omega }_\mathrm{m}=0.33\pm 0.08`$.

### 7.3 Cooling flow models

Under the assumption that the X-ray gas in the core of RXJ1347.5-1145 is in a steady-state cooling flow, we can parameterize the luminosity profile in terms of an equivalent mass deposition rate (e.g. White et al. 1997). The results of this calculation are shown in Fig. 12. We see that the integrated mass deposition rate (MDR) rises with increasing radius within the central 100 kpc, to a maximum value of $`1300`$ $`\mathrm{M}_{}\mathrm{yr}^1`$.

We have examined the spectrum for the central 100 kpc (radius) region using a model appropriate for a cooling flow with distributed mass deposition. (We use the model of Nulsen 1998 and assume that the integrated MDR within radius $`r`$, $`\dot{M}r`$.) The normalization of the cooling-flow component is parameterized in terms of a mass deposition rate, $`\dot{M}`$. The cooling flow was also allowed to be absorbed by an intrinsic column density, $`\mathrm{\Delta }N_\mathrm{H}`$, of cold gas (with solar metallicity) located at the redshift of the cluster. Both $`\dot{M}`$ and $`\mathrm{\Delta }N_\mathrm{H}`$ were free parameters in the fits. Accounting for projection effects, we find that the $`0.57.0`$ keV spectrum for the central 100 kpc radius region is consistent with a cooling flow of $`<600`$ $`\mathrm{M}_{}\mathrm{yr}^1`$, intrinsically absorbed by a column density of $`14\times 10^{21}`$atom cm-2, with a mean mass-weighted flow temperature of $`10`$ keV.

Fitting the mass deposition profile in Fig. 12 with a broken power-law model, we determine a break radius of $`31\pm 2`$ kpc, and slopes internal and external to the break radius of $`2.7\pm 0.3`$ and $`0.5\pm 0.1`$, respectively ($`1\sigma `$ errors). Allen et al. (2001a) argue that such breaks may mark the outer edges of the present-day cooling flows in clusters. The predicted MDR for the central $`r20`$ kpc radius, where a steady cooling flow with gas cooling to zero degrees might be expected to occur, is $`200`$ $`\mathrm{M}_{}\mathrm{yr}^1`$. The cooling time of the X-ray gas at the break radius is $`6\times 10^8`$ yr.

## 8 Discussion

As discussed in Section 4.2, the region of enhanced X-ray emission (and enhanced SZ effect; Komatsu et al. 2001) to the southeast of the main X-ray peak is significantly hotter ($`kT=18\pm 3`$ keV) than at the same distance in other directions ($`kT=12\pm 1`$ keV). Following Markevitch et al. (1999), we can model the properties of this region as a one-dimensional shock. In this case, the relative velocity, $`v`$, of the infalling subcluster and the shock compression factor, $`\rho _1/\rho _0`$, can be estimated from the pre-shock ($`kT_0`$) and post-shock ($`kT_1`$) temperatures. For $`kT_0=12`$ keV and $`kT_1=18`$ keV, we obtain $`v=2250`$ $`\mathrm{km}\mathrm{s}^1`$ and $`\rho _1/\rho _0=1.7`$. If we instead adopt $`kT_0=5`$ keV, a temperature consistent with the predicted mass of the merging subcluster from the strong lensing analysis discussed in Section 6.2 (using the mass-temperature relation of Allen et al. 2001c), then for $`kT_1=18`$ keV we obtain $`v=4550`$ $`\mathrm{km}\mathrm{s}^1`$ and $`\rho _1/\rho _0=3.0`$. We note that the roughly spherical appearance of the shocked gas suggests that this material has expanded following the initial shock. Since the material will cool as it expands, it is likely that the initial post-shock temperature exceeded 18 keV and, for $`kT_0=5`$, that the relative collision velocity exceeded 5000$`\mathrm{km}\mathrm{s}^1`$. Such values are broadly consistent with the peak X-ray surface brightness enhancement of a factor $`3`$ observed at a radius of $`180`$ kpc.

The second brightest galaxy in the cluster lies approximately 18 arcsec to the east and 2 arcsec to the south of the dominant cluster galaxy. This galaxy has an extended, diffuse halo and appears to lie at the centre of a significant sub-concentration of galaxies within the cluster (Fig. 9a). These properties are consistent with it having been the dominant galaxy of the recently merged subcluster. The region of shocked X-ray gas lies approximately 11 arcsec to the south and 6.5 arcsec to the west of the second brightest galaxy. Assuming that the dark matter and galaxies associated with the subcluster have moved on ahead of the shock, it seems plausible that the subcluster travelled in from a PA of $`210`$ degrees. This is consistent with the extension of the X-ray emission observed on large scales to the south (Fig. 1b) and the weak lensing mass map presented by Fischer & Tyson (1997). Given the separation of the shocked gas and second brightest galaxy, we estimate that the shock probably occurred a few $`10^7`$ years ago. This timescale is consistent with the spatial extent of the shocked gas, assuming that it has expanded at the sound speed, $`c_\mathrm{s}1500(T/10^8\mathrm{K})^{1/2}`$$`\mathrm{km}\mathrm{s}^1`$. Note also that this timescale is comparable to the time required for the electrons and ions to reach equipartition, so it is possible that the electron and ion populations in the shocked gas have slightly different temperatures.

Despite the evidence for shocked gas in the southeast quadrant, the Chandra data for the rest of the cluster place RXJ1347.5-1145 on the $`M_{2500}kT_{2500}`$ and $`kT_{2500}L_{2500}`$ relations for relaxed, massive clusters discussed by Allen, Schmidt & Fabian (2001). As discussed in Section 4.4, the shocked gas produces only a small rise in the mean, emission-weighted temperature measured within the central 0.74 Mpc radius. The excess bolometric luminosity of the southeast quadrant between radii of 60-195 kpc (where the shock is most apparent) is $`10^{45}`$ $`\mathrm{erg}\mathrm{s}^1`$, which corresponds to only $`5`$ per cent of the total cluster luminosity (Section 4.4). It therefore appears that the overall temperature and luminosity of the cluster have not been boosted substantially by the merger event. This is supported by the consistent X-ray and weak lensing mass results discussed in Section 6.1.

Ricker & Sarazin (2001) and Ritchie & Thomas (2002) discuss simulations of mergers between clusters with similar central densities and argue that such mergers can lead to significant short-term increases in the overall X-ray temperatures and luminosities of clusters during the periods of closest approach. These effects are most pronounced and long-lived in clusters with relatively low central densities ($`n_\mathrm{e}<`$ few $`10^3`$cm-3; Ritchie & Thomas 2002). The fact that large boosts in the overall temperature and luminosity are not observed in RXJ1347.5-1145 may be related to the very high central density in the main cluster core ($`n_\mathrm{e}>10^1`$cm-3; Fig. 10a) and could indicate that the recently-merged subcluster had a relatively low central gas density and interacted only weakly with the main cluster core.

Recently, Cohen & Kneib (2002) have published spectroscopic redshift measurements for 47 cluster members for which they determine a velocity dispersion of $`910\pm 130`$$`\mathrm{km}\mathrm{s}^1`$. Inspection of their Fig. 2 suggests that a significant fraction of the galaxies with measured redshifts may have been associated with the recently merged subcluster. This could explain the relatively low velocity dispersion, in comparison to the X-ray and lensing results.

The agreement between the Chandra X-ray and gravitational lensing mass measurements for RXJ1347.5-1145 reinforces the results and conclusions drawn from previous studies of other relaxed lensing clusters e.g. Abell 2390 (Allen et al. 2001b) and Abell 1835 (Schmidt et al. 2001); see also the results for MS1358.4+6245 by Arabadjis et al. (2002). The close agreement of the independent X-ray and lensing masses indicates that the mass measurements are robust and limits the contributions from non-thermal sources of pressure support in the X-ray gas, such as bulk and/or turbulent motions and magnetic fields, in the relaxed regions of the cluster.

## 9 Conclusions

The main conclusions from this work may be summarized as follows:

(i) We have reported Chandra observations of RXJ1347.5-1145, the most X-ray luminous cluster of galaxies known. We have identified a region of shocked gas (enhanced X-ray brightness and temperature) to the southeast of the main X-ray peak, at a position consistent with the region of enhanced SZ effect reported by Komatsu et al. (2001). The shocked gas probably results from recent subcluster merger activity. The merger appears not to have boosted the overall luminosity and temperature of the cluster substantially.

(ii) Excluding the data for the southeast quadrant, we have measured the density, temperature and mass profiles for the cluster. The mass profile can be parameterized by an NFW model with a scale radius $`r_\mathrm{s}=0.40_{0.12}^{+0.24}`$ Mpc and a concentration parameter, $`c=5.9\pm 1.4`$ (68 per cent confidence limits). The normalization of the mass profile may also be expressed in terms of an effective velocity dispersion, $`\sigma =\sqrt{50}H(z)r_\mathrm{s}c=1450_{200}^{+300}`$$`\mathrm{km}\mathrm{s}^1`$.

(iii) The best-fit Chandra mass model for RXJ1347.5-1145 is in good agreement with independent measurements from weak gravitational lensing studies (Fischer & Tyson 1997). The observed strong lensing configuration in the cluster core can also be explained with the introduction of an additional mass clump centred on the second brightest galaxy, which is likely to have been the dominant galaxy of the recently merged subcluster.

## Acknowledgements

We thank R. Johnstone and S. Ettori for supplying a number of scripts used in the X-ray reduction/analysis and Philippe Fischer for communicating his weak lensing results in a machine readable form. SWA and ACF acknowledge the support of the Royal Society.


format: text

text:
# Statistical Properties of Fano Resonances in Atomic and Molecular Photoabsorption

## Abstract

Statistical properties of Fano resonances occurring in photoabsorption to highly excited atomic or molecular states are derived. The situation with one open and one closed channel is analyzed when the classical motion of the excited complex in the closed channel is chaotic. The closed channel subspace is modeled by random matrix theory. The probability distribution of the Fano parameter is derived both for the case of time reversal symmetry (TRS) and broken time reversal symmetry. For the TRS case the area distribution under a resonance profile relevant for low resolution experiments is discussed in detail.

Introduction — Resonance phenomena are ubiquitous in highly excited quantum systems. Typical examples are auto-ionizing atomic resonances Fri98  or photodissociation of rovibrational molecular states Schi93 . When the density of resonances is high a common and natural way to extract information about resonance properties is to resort to statistical methods Haake01 ; Stoeck99 ; Guhr98 . Typically as a consequence of the high excitation the resonant complex behaves chaotically from a classical point of view. It is then appropriate to model the excited complex by means of random matrix theory BGS84 ; Alhas92 ; Alhas98 .

Here I address a problem frequently encountered in photoabsorption processes of atomic and molecular systems: What are the statistical properties of line shapes when photoabsorption into a highly excited resonant state takes place? Since the seminal work of Fano Fano61  it is known that when at least two pathways exist for an atomic or molecular system to decay after having absorbed a photon the line shapes of resonances differ in general considerably from the Breit-Wigner form. In photodissociation of molecules or photoexcitation of an atom to an auto-ionizing resonance interference between the indirect and the direct decay process results in a Beutler-Fano resonance profile of the photoabsorption cross section Fri98 ; Conn88 ; Durand01  (abbreviated as Fano profile in the following). The shape of the resonant part of the cross section is determined by a single parameter, the Fano parameter $`q`$. This contribution is aimed at such a statistical theory of Fano resonances.

It is worthwhile to stress that atomic and molecular systems are ideally suited to test the predictions on Fano resonances as neither Coulomb or dephasing effects which may become relevant in mesoscopic structures have to be taken into account Clerk01 .

Theory — For photodissociation of a molecule Fig. 1 schematically depicts the situation envisaged footnt1 . The molecule is coherently laser excited from the ground state or a low lying state $`|0`$ of energy $`E_0`$ onto two electronic potential surfaces. The transition is accomplished by the dipole operator $`D`$. Electronic surface S1 is the open channel — the energy $`E`$ of the molecular complex is above the dissociation threshold — while the motion on the potential surface S2 is bound (closed channel). The dipole matrix element for transition to the open channel is given by $`d_10|D|\psi _E^{(\mathrm{reg})}`$, where $`|\psi _E^{(\mathrm{reg})}`$ is the regular continuum solution at energy $`E`$ in the open channel, normalized in energy. When $`E`$ coincides with the energy $`E_n`$ of a bound state $`|\varphi _n`$ in the channel 2 the closed channel carries the transition amplitude $`0|D|\varphi _n`$ in the absence of coupling to the continuum. In the following it is assumed that transitions between the two excited manifolds are possible. In the diabatic representation the coupling between the two surfaces is given by a non-diagonal potential $`V`$. Then the eigenstate $`|\varphi _n`$ turns into a resonance with width $`\mathrm{\Gamma }_n=2\pi |\psi _E^{(\mathrm{reg})}|V|\varphi _n|^2`$.

In the regime $`\mathrm{\Gamma }_n\mathrm{\Delta }`$ of isolated resonances ($`\mathrm{\Delta }`$ being the mean energy spacing of resonances) the oscillator strength $`Df(E)`$ for the dipole transition is given by Fri98

$$Df(E)=Df^{(\mathrm{bg})}\frac{|q_n+ϵ_n|^2}{1+ϵ_n^2},$$
(1)
where $`Df^{(\mathrm{bg})}2\mu (EE_0)|d_1|^2/\mathrm{}^2`$ is the background oscillator strength for the transition to the open channel in the absence of coupling to the closed channel, $`\mu `$ is the reduced mass of the excited complex and $`ϵ_n2(EE_n\mathrm{\Delta }_n)/\mathrm{\Gamma }_n`$ the reduced energy. The energy shift $`\mathrm{\Delta }_n`$ of the resonance is neglected in the following. The line shape of the resonance is parameterized by the Fano parameter $`q_n`$. For a non-degenerate continuum it is given by Fri98

$$q_n=\frac{0|D|\varphi _n}{\pi d_1\psi _E^{(\mathrm{reg})}|V|\varphi _n}+\overline{q},\overline{q}=\frac{0|D|\psi _E^{(\mathrm{irr})}}{0|D|\psi _E^{(\mathrm{reg})}}.$$
(2)
The irregular continuum solution $`|\psi _E^{(\mathrm{irr})}`$ appears in the expression for $`q_n`$ because near resonance the wave function in the open channel is a superposition of the regular and the irregular solution. Both $`d_1`$ and $`\overline{q}`$ depend only weakly on $`E`$ and are assumed to be constant within the energy window over which a sample of resonances is taken. From (2) it is seen that the distribution of a set of Fano parameters $`\{q_n\}`$ taken over a stretch of the energy spectrum is determined by the statistical properties of the set of the eigenstates $`\{|\varphi _n\}`$ of the closed channel.

The distribution of the Fano parameter is derived under the two following conditions: Firstly the classical motion of the excited molecular complex on the electronic surface S1 of the closed channel is completely chaotic. This ensures that the statistical properties of generic wave functions in chaotic systems apply and that the closed channel subspace can be modeled by random matrix theory Haake01 . Secondly the excitation process and the coupling between the two electronic surfaces are assumed to be spatially well separated: the overlap of the initial wave packet $`𝐫|D|0`$ and the coupling potential $`V(𝐫)`$ in coordinate representation is negligible. Therefore $`x0|D|\varphi _n/d_1`$ and $`y\pi \psi _E^{(\mathrm{reg})}|V|\varphi _n`$ can be taken statistically independent random variables footnt2 .

Distribution of the Fano parameter — With these ingredients the calculation of the probability distribution $`P(q)`$ of the Fano parameter is straightforward. (The index $`n`$ of the resonance is omitted in the following.) In the case of time reversal symmetry the wave functions and viz. $`x`$ and $`y`$ can be chosen real. The closed channel subspace is modeled by $`N\times N`$ matrices taken from the Gaussian orthogonal ensemble (GOE). All results are understood in the limit $`N\mathrm{}`$ where $`x`$ and $`y`$ are Gaussian random random variables with zero mean and variances $`\sigma _x^2`$ and $`\sigma _y^2`$. The probability distribution is given by

$$P_{\mathrm{GOE}}(q)=\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑x\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑yP_{\sigma _x}(x)P_{\sigma _y}(y)\delta \left(q\overline{q}\frac{x}{y}\right)$$
(3)
and $`P_{\sigma _x}(x)`$ is the probability distribution of $`x`$ (and likewise for $`y`$). The integral is most easily evaluated by using the Fourier integral representation of the delta function and the result is

$$P_{\mathrm{GOE}}(q;s)=\frac{1}{\pi }\frac{s}{s^2+(q\overline{q})^2},s\sigma _x/\sigma _y.$$
(4)
The probability distribution of the Fano parameter in the GOE case thus turns out to be Lorentzian with mean value $`\overline{q}`$ and width $`s`$.

The width $`s`$ of the probability distribution is related to the coupling strength $`V`$ between the closed channel and the continuum channel and the ratio of the dipole transition matrix elements to both channels. Assume $`V`$ can be written in the form $`V=\lambda V_0`$ where $`\lambda `$ characterizes the coupling strength and $`V_0`$ is fixed. Then $`\sigma _y\lambda `$ and therefore $`s1/\lambda `$. For strong coupling to the continuum the Lorentz distribution acquires a small width centered around $`\overline{q}`$. The same holds if direct photoexcitation dominates over the indirect process since then $`\sigma _x`$ becomes small.

In the case of broken time reversal symmetry the Hamiltonian of the closed channel is modeled by $`N\times N`$ matrices from the Gaussian unitary ensemble (GUE). In this case $`x`$ and $`y`$ are complex Gaussian random variables with independent real and imaginary parts and the Fano parameter is in general complex. Most conveniently the distribution of $`q`$ is characterized by the probability distribution of the phase $`\phi _q`$ and the modulus $`r_q`$ of the quantity $`qq`$. The phases $`\phi _x`$ of $`x`$ and $`\phi _y`$ of $`y`$ are uniformly distributed $`\mathrm{mod}(2\pi )`$ and the same holds for $`\phi _q=\phi _x\phi _y`$. Denoting the modulus of $`x`$ by $`r_x`$ and of $`y`$ by $`r_y`$ the probability distribution $`P_{\mathrm{GUE}}(r_q)`$ is given by

$$P_{\mathrm{GUE}}(r_q)=\underset{0}{\overset{\mathrm{}}{}}𝑑r_x\underset{0}{\overset{\mathrm{}}{}}𝑑r_yP(r_x)P(r_y)\delta \left(r_q\frac{r_x}{r_y}\right),$$
(5)
where $`r_x`$ (and likewise $`r_y`$) has the probability distribution $`P(r_x)=\sigma _{r_x}^2r_x\mathrm{exp}(r_x^2/2\sigma _{r_x}^2)`$. \[Notice that $`P(r_x^2)\mathrm{exp}(r_x^2/2\sigma _{r_x}^2)`$ has the form of a Porter-Thomas distribution for GUE.\] Proceeding as in the GOE case the integrals turn out to be more complicated as they involve error functions. The final result is nevertheless simple and reads

$$P_{\mathrm{GUE}}(r_q)=\delta (r_q)+\frac{s_r^2r_q}{\left(s_r^2+r_q^2\right)^2},(r_q0),$$
(6)
where $`s_r\sigma _{r_x}/\sigma _{r_y}`$. The distribution is delta peaked at $`r_q=0`$ and has a local maximum at $`r_q^{(0)}=s_r/\sqrt{3}`$ with $`P(r_q^{(0)})=3\sqrt{3}/(16s_r)`$. Again, as discussed for the GOE case the width $`s_r`$ of the probability distribution $`P_{\mathrm{GUE}}(r_q)`$ is determined by the strength of the coupling between the two channels and the ratio of the strength between the dipole transition to the closed and to the open channel.

Peak height distribution for GOE — For the TRS case the maximum of the Fano profile is given at $`ϵ=1/q`$ with peak height $`z[DfDf^{(\mathrm{bg})}]/Df^{(\mathrm{bg})}=q^2`$ relative to the background oscillator strength. The distribution of peak heights can easily be calculated from (4) and is given by

$$P_{\mathrm{GOE}}(z)=\frac{1}{2\sqrt{z}}\left[P_{\mathrm{GOE}}(\sqrt{z};s)+P_{\mathrm{GOE}}(\sqrt{z};s)\right]$$
(7)
with $`z>0`$. In Fig. 2 $`P_{\mathrm{GOE}}(z)`$ is plotted in terms of the scaled variable $`z/\overline{q}^2`$ and different values of the parameter $`\beta =\overline{q}/s`$. It has a local maximum at a value $`0<z/\overline{q}^2<1`$ for $`\beta 1.95517`$. Therefore for a local maximum to occur the width $`s`$ of the Lorentz profile (4) must be smaller than approximately twice the average value $`\overline{q}`$ of the Fano parameter. For $`\beta \mathrm{}`$ the position of the maximum approaches $`z=\overline{q}^2`$.

Profile area distribution for GOE — In low resolution experiments the relevant experimental quantity is the profile area distribution under the resonance rather than the peak height distribution. To be more specific assume that the band width $`\delta E`$ of the exciting laser beam is larger than the width $`\mathrm{\Gamma }`$ of the resonance but still smaller than the mean spacing $`\mathrm{\Delta }`$ between adjacent resonances ($`\mathrm{\Gamma }<\delta E<\mathrm{\Delta }`$). Assuming a rectangular laser profile the excess oscillator strength averaged over the resonance profile with respect to the background strength is given by Bohm86

$$\overline{Df}^{(\mathrm{res})}(E)=\frac{1}{\delta E}\underset{E\delta E/2}{\overset{E+\delta E/2}{}}Df^{(\mathrm{bg})}(E^{})\left[F(q;ϵ)1\right]𝑑E^{}$$
(8)
where $`F(q;ϵ)|q+ϵ|^2/(1+ϵ^2)`$ is the profile function of the resonance \[cf. (1)\]. Since $`\mathrm{\Gamma }<\delta E`$ the range of integration in (8) can be extended to infinity and $`\overline{Df}^{(\mathrm{res})}(E)`$ is given by

$$\overline{Df}^{(\mathrm{res})}(E)=Df^{(\mathrm{bg})}\frac{\mathrm{\Gamma }}{2\delta E}\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{q^2+2qϵ1}{1+ϵ^2}𝑑ϵ,$$
(9)
which performing the integration gives $`\overline{Df}^{(\mathrm{res})}/Df^{(\mathrm{bg})}=\pi \mathrm{\Gamma }(q^21)/2\mathrm{\Theta }`$. In the following the statistical distribution of the observable $`\mathrm{\Theta }`$ is discussed. It can be written as $`\mathrm{\Theta }=x^2+2\overline{q}xy+(\overline{q}^21)y^2f(x,y)`$ in terms of the matrix elements $`x`$, $`y`$ and the mean value $`\overline{q}`$ of the Fano parameter. The probability distribution of $`\mathrm{\Theta }`$ is given by

$$P(\mathrm{\Theta })=\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑xP_{\sigma _x}(x)\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑yP_{\sigma _y}(y)\delta \left(\mathrm{\Theta }f(x,y)\right).$$
(10)
Introducing the scaled area $`\widehat{\mathrm{\Theta }}a\mathrm{\Theta }`$ with $`a=(\sigma _x\sigma _y)^1\sqrt{1/4+b^2}`$ and $`b=\left[s+s^1(\overline{q}^21)\right]/4`$ the distribution $`P(\widehat{\mathrm{\Theta }})`$ can be written as

$$P(\widehat{\mathrm{\Theta }})=\frac{1}{\pi \sqrt{1+4b^2}}\mathrm{exp}(\mu \widehat{\mathrm{\Theta }})K_0(|\widehat{\mathrm{\Theta }}|),$$
(11)
where $`K_0`$ is the MacDonald function. Additionally the asymmetry parameter $`\mu =b[\frac{1}{4}+b^2]^{1/2}`$ has been introduced which determines how much $`P(\widehat{\mathrm{\Theta }})`$ deviates from a symmetric distribution. $`P(\widehat{\mathrm{\Theta }})`$ diverges logarithmically as $`\widehat{\mathrm{\Theta }}`$ approaches zero.

Fig. 3 demonstrates the dependence of $`P(\widehat{\mathrm{\Theta }})`$ at fixed $`\overline{q}=0`$ in the three characteristic regimes $`\mu <0`$, $`\mu =0`$ and $`\mu >0`$ by varying $`s`$. For $`\mu =0`$ the probability distribution is symmetric with respect to $`\widehat{\mathrm{\Theta }}=0`$ (solid line). For $`\mu <0`$ the area distribution is asymmetric with the flat side at negative values of $`\widehat{\mathrm{\Theta }}`$ (dashed line). The opposite holds for $`\mu >0`$ (dot-dashed line).

As $`\mu `$ depends both on $`s`$ and $`\overline{q}`$ a contour plot of $`\mu `$ as a function of $`(\overline{q},1/s)`$ is presented in Fig. 4. The variable $`1/s`$ instead of $`s`$ has been chosen for matters of discussion as $`1/s\mathrm{}`$ corresponds to the limit of strong coupling between the closed and the continuum channel. The range of $`\mu `$ is $`1<\mu <1`$. The white dashed line marks the subspace of parameters where $`P(\widehat{\mathrm{\Theta }})`$ is symmetric $`(\mu =0)`$, given by $`\overline{q}^2+s^2=1`$. If $`|\overline{q}|1`$ the asymmetry parameter $`\mu `$ is always positive regardless of $`1/s`$. For fixed $`|\overline{q}|<1`$ and weak coupling ($`1/s`$ small) the asymmetry parameter $`\mu `$ is positive, too. If the coupling is enlarged $`\mu `$ becomes negative above a certain value of $`1/s`$ which depends on $`\overline{q}`$.

The expectation value of the scaled area distribution is given by $`\widehat{\mathrm{\Theta }}=\mu (1+4b^2)`$. Thus $`\widehat{\mathrm{\Theta }}`$ is always positive when $`|\overline{q}|1`$. Note that the area under an individual resonance is positive if for its Fano parameter $`|q|>1`$ holds. When $`|q|<1`$ the area associated with an individual resonance is negative. In contrast $`\widehat{\mathrm{\Theta }}`$ can either be positive or negative for $`|\overline{q}|<1`$ depending on the coupling strength to the continuum. Strong coupling to the continuum ($`1/s\mathrm{}`$) results in a negative value of $`\widehat{\mathrm{\Theta }}`$. In the limit $`\mu +1`$ Breit-Wigner resonances dominate the spectrum. For fixed $`|\overline{q}|>1`$ this limit is reached both for weak and large coupling to the continuum. For $`|\overline{q}|<1`$ the Breit-Wigner limit is reached only for weak coupling ($`1/s0`$). For $`|\overline{q}|<1`$ and strong coupling to the continuum ($`1/s\mathrm{}`$) the expectation value of the area distribution is negative and window resonances with negative area dominate the spectrum in the limit $`\mu 1`$.

Summary — Statistical properties of Fano resonances were derived for the situation of competing direct and indirect photoabsorption on the basis that the classical motion in the closed channel is chaotic. This situation is frequently encountered in photodissociation of highly excited rovibrational states of molecules or photoabsorption of atoms to auto-ionizing resonances. The closed channel subspace is modeled by random matrix theory. The distribution of the Fano parameter $`q`$ was derived for both the GOE and the GUE cases. The peak height distribution and the resonance area distribution were discussed in detail for the GOE.

###### Acknowledgements.

Acknowledgments — I am grateful to John S. Briggs for pointing out to me the problem of Fano resonances in classically chaotic systems. I am indebted to Thomas Seligman for discussions and his hospitality during a visit at CIC, Cuernavaca, Mexico, where part of this work was completed. Discussions with J. Flores, T. Gorin, B. Mehlig and M. Müller are appreciated. Financial support from SFB 276 (“Correlated dynamics of highly excited atomic and molecular systems”) is acknowledged.


text:
# Crystal structure of solid Oxygen at high pressure and low temperature

## Abstract

Results of X–ray diffraction experiments on solid oxygen at low temperature and at pressures up to 10 GPa are presented. A careful sample preparation and annealing around 240 K allowed very good diffraction patterns in the orthorhombic $`\delta `$–phase to be obtained. This phase is stable at low temperature, in contrast to some recent data \[Y. Akahama et al., Phys. Rev. B64, 054105 (2001)\], and transforms with decreasing pressure into a monoclinic phase, which is identified as the low pressure $`\alpha `$–phase. The discontinuous change of the lattice parameters, and the observed metastability of the $`\alpha `$–phase increasing pressure suggest that the transition is of the first order.

The magnetic interaction between molecules is of fundamental importance for determining the structural properties of solid oxygen, and, in particular, is responsible for the stabilization of the anti-ferromagnetic $`\alpha `$ phase. Indeed, calculations demonstrate that, in absence of magnetic interaction, the lowest energy structure at zero temperature and ambient pressure would be a planar structure, consisting of an hexagonal close packing of parallel oxygen molecules with centers of mass on the planes, and orientation perpendicular to the planes. This structure corresponds to the actual oxygen’s $`\beta `$–phase (rhombohedral, space group $`R\overline{3}m`$), which however is stable only above 23.9 K (up to 43.8 K). At lower temperature the stable structure ($`\alpha `$–phase) is instead monoclinic (space group $`C2/m`$) , and presents an anti-ferromagnetic order of the molecular spins, which lay in the ab plane. The $`\alpha `$–phase has the same planar arrangement of parallel molecules. At the $`\alpha `$$`\beta `$ transition, therefore, both structural and magnetic changes occur, and the transition has been described as being of first order, and driven by the magnetic interaction. The distortion of the regular hexagons formed by the centers of mass of the molecules in the $`\beta `$–phase is due to the magnetic interaction and determines an increase of the $`b/a`$ ratio, which is $`1/\sqrt{3}=0.577`$ in the $`\beta `$–phase and $`0.63`$ in the $`\alpha `$–phase at low pressure.

Three different phases have been detected by Raman spectroscopy in the pressurized solid at room temperature, and their structure has been determined by X–ray diffraction several years ago. Fluid oxygen crystallizes at room temperature into a rhombohedral structure (stable between 5.5 and 9.6 GPa) which coincides with the low temperature $`\beta `$–phase, and the phase stable between 9.6 and 9.9 GPa ($`\delta `$–phase) is orthorhombic, with space group $`Fmmm`$. The structure of the $`ϵ`$–phase, stable above 9.9 GPa at room temperature, is monoclinic , with space group $`C2/m`$, but the positions of the molecules in the unit cell have not been determined.

The structure of the low temperature monoclinic $`\alpha `$–phase and of the orthorhombic $`\delta `$–phase are very closely related and may in principle transform one into the other by a second order transition consisting in a sliding of the basal planes (see fig. 1 of ref. ). A lattice dynamic calculation predicts instead a first order transition, occurring at about 2.4 GPa. Experimental information on the oxygen phases at low temperature has been derived only indirectly, and it has been summarized by Desgreniers et al.. Indications of the $`\alpha `$$`\delta `$ transition at about 1.0 GPa have been initially reported from low temperature Raman experiments but other authors describe a phase transition around 0.5 GPa ($`\alpha \alpha ^{}`$) and another one around 3 GPa ($`\alpha ^{}\delta ^{}`$). Yen and Nicol observe indications of a phase transition at higher pressure and temperature ($`5`$ GPa and 120 K).

Direct X–ray diffraction structural studies at low temperature have not been attempted until recently. The authors, increasing pressure along the isotherms at 19, 180 and 240 K, observed that at all the three temperatures the monclinic $`\alpha `$–phase transforms directly into the $`ϵ`$–phase, and therefore confined the range of existence of the $`\delta `$—phase in the unexplored region above 240 K. The measured monoclinic distortion is however very small, and possible hysteresis effect may have prevented the determination of the $`\alpha `$$`\delta `$ transition.

We have recently performed infrared spectroscopy studies of oxygen at low temperature and in the pressure range up to 8 GPa. From the detection of an infrared absorption peak in the frequency region of the fundamental vibration of the oxygen molecule, and, in particular, from its evolution with temperature, we demonstrated that solid oxygen is anti-ferromagnetic up to high pressure. This result has been confirmed recently by a study of the low energy electronic transitions.

In this paper we present X–ray diffraction results for pressure up to 10 GPa from which the problem of the structure at low temperature can be settled. Because of a careful annealing procedure, very good diffraction patterns in the orthorombic $`\delta `$–phase were obtained. The main conclusion of this study is the observation of the orthorhombic $`Fmmm`$ phase at temperatures lower than 300 K. This phase is stable decreasing $`T`$, and transforms to the monoclinic $`\alpha `$–phase decreasing $`P`$.

We have used a diamond anvil cell (DAC) of the membrane type, loaded condensing high purity liquid oxygen in a sealed vessel. The cell is fixed to the cold finger of a He flux cryostat. Temperature is determined by a silicon diode, placed very close to one diamond, with a good reliability ($`\pm 1`$ K). Pressure is measured by the shift of the ruby fluorescence wavelength, using the calibration given in ref. , corrected for measurements at low temperature. Sample dimensions ranged from 150 to 170 $`\mu `$m, with a thickness of about 50 $`\mu `$m.

Angle dispersive X–ray diffraction has been measured at the ID09 beamline of the European Synchrotron Radiation Facility (ESRF, Grenoble), with a monochromatic beam, ($`\lambda =0.4171`$Å) and image–plate detection. The beam reaches the sample through one diamond, forming a square spot with dimension $`40\times 40\mu `$m. We used a beryllium seat for the second diamond, to obtain a maximum diffraction angle $`2\theta 25`$ degrees. The diffraction patterns were analyzed and integrated by means of the FIT2D computer code, to obtain the one dimensional intensity distribution as a function of the $`2\theta `$ scattering angle.

Pressurizing fluid oxygen at room temperature into the solid $`\beta `$–phase usually results in the formation of large crystals in the DAC. To produce a fine–grained polycrystalline powder, we have caused an abrupt transition of the fluid sample into the $`ϵ`$–phase, with a very rapid pressure increase up to about 40–50 GPa. Subsequently, the $`\delta `$$`ϵ`$ phase transition, where solid oxygen experiences a large volume change has been crossed a few times in both directions, at a temperature of about 270 K. With this procedure the preferential orientation of the powder in the $`ϵ`$–phase is greatly reduced, as monitored by the almost uniform angular distribution of the intensity along the diffraction rings.

Starting from the $`ϵ`$–phase at a temperature of about 270 K, we have slowly decreased the pressure until the $`ϵ`$$`\delta `$ phase transition was observed. The diffraction pattern obtained after the transition is represented in fig. 1(a). This can be unequivocally indexed on the basis of an orthorhombic cell. Decreasing slowly the temperature and the pressure, along a stepwise path made up of subsequent isotherms and isobars, we have found that the orthorhombic structure is maintained in a quite large region of the phase diagram. Indeed, a monoclinic distortion ($`\beta 90^{}`$) would imply the splitting of the 111, 113, 311, 204, 222 reflections, which instead remain quite narrow (see fig. 1(b)). Other cycles performed following different paths in the $`PT`$ diagram have demonstrated that a temperature of about 270 K or above is necessary to obtain the annealing of the crystal in the orthorhombic structure in a reasonable time.

Starting from the orthorhombic phase at low temperature, we have collected diffraction patterns during several isothermal decompression runs. To detect any possible monoclinic deformation we have fitted all the structures leaving the monoclinic angle $`\beta `$ as a free parameter. The diffraction patterns obtained during decompression from 6.5 to 5.3 GPa and on compression from 5.3 up to 6.2 GPa at 65 K are shown in fig. 1(c). Between 5.9 and 5.3 GPa the diffraction pattern changes drastically with a big shift to lower angle of the 200 line and with a splitting of the 111 line into 111 and $`11\overline{1}`$. This is an unambiguous evidence of the transition to the monoclinic phase. The diffraction pattern at 5.3 GPa can be satisfactorily indexed on the basis of a monoclinic structure with space group $`C2/m`$, the same of the $`\alpha `$–phase. Upon compression at low temperature, the monoclinic structure does not transform again to the orthorhombic one, as testified by the two diffraction patterns measured at the same pressure, 5.9 GPa, before and after the phase transition. To obtain the orthorhombic structure again the temperature had to be increased up to about 270 K, increasing accordingly the pressure in order to avoid the transformation to the $`\beta `$ phase. To our knowledge in the previous spectroscopic experiments the pressure was increased at low temperature, producing the metastability of the monoclinic phase. The same argument can be used to explain the reason for the missed observation of the phase transition in the diffraction experiment by Akahama et al. .

The enlarged monoclinic cell used in the fit of the diffraction patterns (the $`abc^{}`$ polyhedron in the inset of fig. 2(a)) contains 4 molecules and describes as a particular case the orthorhombic one. In this sense the evolution of the $`\beta ^{}`$ angle is of particular interest. This is reported in fig. 2(a) as a function of pressure along the isotherms that we have performed. At pressures above 6–7 GPa, depending on temperature, the values of the $`\beta ^{}`$ angle are all very close to 90 degrees. Decreasing pressure, the value of $`\beta ^{}`$ shows a sudden increase of about 6 degrees, that is the hallmark of the $`\delta `$$`\alpha `$ transition. The solid line in fig. 2 represents the values of $`\beta ^{}`$ resulting from lattice dynamics calculation performed at 0 K , plotted after having multiplied the calculated pressure by a factor, to make the calculated transition (2.4 GPa) pressure to coincide with the experimental one (about 5.5 GPa at 85 K). There is a quantitative agreement with the experiment both for the value of the change of $`\beta ^{}`$ at the transition and for the pressure slope in the $`\alpha `$–phase.

As mentioned before, the anti-ferromagnetic order of the $`\alpha `$–phase produces the distortion of the regular hexagon marking the molecular positions in the $`\beta `$–phase, and the departure of the ratio $`b/a`$ from the value $`b/a=1/\sqrt{3}`$ is a measure of the magnetic interaction between molecules. It is important therefore to monitor the evolution with pressure of the ratio $`b/a`$, whose values along the isotherms at 65, 85 and 145 K are reported in fig. 2(b,c,d), and compared to the lattice dynamics calculation data , plotted, also in this case, after having rescaled the theoretical pressure values. The ratio $`b/a`$ increases steadily with pressure, in both the $`\alpha `$ and $`\delta `$ phases, confirming the increase of the magnetic interaction between molecules with pressure.

The calculation is quantitatively correct in the $`\alpha `$–phase, and also the positive sudden increment (if not the size) at the $`\alpha `$$`\delta `$ transition is predicted by the theory.

The main conclusions concerning the low temperature phase diagram are summarized in fig. 3. Full and open circles mark the $`PT`$ values where diffraction has been measured, representing respectively the orthorhombic and monoclinic structures. The line represents a tentative phase boundary between $`\alpha `$ and $`\delta `$ phases. Due to the fact that no evidence of the phase boundary $`\delta `$’–$`\delta `$ was observed and that the structure is monoclinic with the same space group of the $`\alpha `$ phase, $`C2/m`$, we assign the whole region enclosed by the $`\beta `$ and $`\delta `$ phases to the $`\alpha `$ phase. The phase boundary identified by Yen and Nicol  on the basis of weak extra peaks in the low frequency region of the Raman spectrum is very close to the phase boundary determined in this work. For this reason we did not remove this phase line, as this might indicate the same $`\alpha `$$`\delta `$ phase transition.

To conclude, our results demonstrate the stability of the orthorhombic $`\delta `$–phase below room temperature, that can observed when an annealing of the crystal is performed at room temperature or slightly below, but not when a sample is decompressed from the $`ϵ`$–phase at low temperature across the $`ϵ`$$`\delta `$ transition. In this latter case a mixture of orthorhombic and monoclinic structure is obtained. The $`\alpha `$$`\delta `$ transition has been observed unambiguously, starting with a good orthorhombic sample and decreasing pressure along four isotherms, appearing with jump in the monoclinic angle $`\beta ^{}`$ from 90 degrees to about 96 degrees. The sign and size of this change, and the slope with which the angle $`\beta ^{}`$ changes with pressure in the $`\alpha `$–phase (but not the transition pressure) agree in a remarkable way with the theoretical prediction of Etters et al. .

The lattice parameters ratio $`b/a`$ also shows a discontinuous increase passing from $`\alpha `$–O2 to $`\delta `$–O2, also predicted by the calculations .

No other phase transitions are detected below 5 GPa. The phase diagram of solid oxygen therefore has been completely determined, deleting uncertain phase transition boundaries at low temperature. The sensitivity of the X-ray diffraction technique has been essential to clarify the weakness of previous spectroscopic work, where phase transitions can be detected only by subtle kinks in the pressure slope of mode frequencies. In addition, the results of this work can furnish a precise recipe to be followed to avoid metastability and obtain a good annealing of an solid oxygen sample.

This research has been supported by the EU contract HPRICT1999-00111.


text:
# Nontrivial Velocity Distributions in Inelastic Gases

## Abstract

We study freely evolving and forced inelastic gases using the Boltzmann equation. We consider uniform collision rates and obtain analytical results valid for arbitrary spatial dimension $`d`$ and arbitrary dissipation coefficient $`ϵ`$. In the freely evolving case, we find that the velocity distribution decays algebraically, $`P(v,t)v^\sigma `$ for sufficiently large velocities. We derive the exponent $`\sigma (d,ϵ)`$, which exhibits nontrivial dependence on both $`d`$ and $`ϵ`$, exactly. In the forced case, the velocity distribution approaches a steady-state with a Gaussian large velocity tail.

PACS: 47.70.Nd, 45.70.Mg, 05.40.-a, 05.20.Dd, 81.05.Rm

The dynamics of granular gases present novel challenges, previously not encountered in fluid dynamics . Specifically, the strong underlying energy dissipation leads to clustering instabilities and strong velocity correlations . A series of recent experimental and theoretical studies reveals a rich phenomenology. In particular, velocities are characterized by anomalous statistics, sensitive to the details of the driving conditions, the density, and the degree of dissipation .

Kinetic theory provides a systematic framework for deriving macroscopic properties from the microscopic collision dynamics . Yet, analysis of the corresponding Boltzmann equation often involves uncontrolled approximations or use of nearly Maxwellian distributions. Motivated by the latter issue, we examine both unforced and forced inelastic gases using a simplified Boltzmann equation. Specifically, we shall employ Maxwell’s collision rate which is proportional to the typical velocity rather than the relative velocity, and show that in general, it is analytically tractable.

Previous studies reveal that this kinetic theory leads to interesting behaviors in the freely evolving case. In one dimension, while moments of the velocity distribution exhibit multiscaling , the velocity distribution itself still approaches a scaling form with an algebraic large velocity tail . An algebraic tail was also found numerically in two dimensions . Here, we show analytically that in arbitrary spatial dimension the velocity distribution admits a scaling solution with an algebraic large velocity tail. The corresponding exponent, a root of a transcendental equation, depends on the spatial dimension and the restitution coefficient. In the driven case, we find that although the velocity distribution is non-Maxwellian, it does exhibit a Gaussian tail.

Our starting point is a homogeneous system of identical inelastic spherical particles in arbitrary spatial dimension $`d`$. The mass and the cross-section are set to unity without loss of generality. When two particles collide, the normal component of the relative velocity is reduced by the restitution coefficient $`r=12ϵ`$, while the tangential component remains the same. Denoting by $`𝐧`$ the unit vector connecting the centers of the colliding particles, the post-collision velocities $`𝐯_{1,2}`$ are given by a linear combination of the precollision velocities $`𝐮_{1,2}`$,

$$𝐯_{1,2}=𝐮_{1,2}(1ϵ)(𝐠𝐧)𝐧,$$
(1)
with the relative velocity $`𝐠=𝐮_1𝐮_2`$. The energy dissipated in a collision equals $`\mathrm{\Delta }E=ϵ(1ϵ)(𝐠𝐧)^2`$. When $`ϵ=0`$ collisions are elastic, while when $`ϵ=1/2`$ collisions are perfectly inelastic with maximal energy dissipation.

Let $`P(𝐯,t)`$ be the normalized density of particles with velocity $`𝐯`$ at time $`t`$. Assuming molecular chaos or perfect mixing, we arrive at the following Boltzmann equation for the velocity distribution

$$\frac{P(𝐯,t)}{t}=g𝑑𝐧𝑑𝐮_1𝑑𝐮_2P(𝐮_1,t)P(𝐮_2,t)\left\{\delta \left[𝐯𝐮_1+(1ϵ)(𝐠𝐧)𝐧\right]\delta (𝐯𝐮_1)\right\}.$$
(2)
The collision rate in the Maxwell approximation represent the typical velocity scale, $`g=\sqrt{T}`$, with $`T=\frac{1}{d}v^2`$ the “granular temperature” or the average kinetic energy per degree of freedom (for hard spheres, the collision rate equals the actual relative velocity $`𝐠𝐧`$). This evolution equation naturally describes a stochastic process where randomly chosen pairs of particles undergo inelastic collisions according to (1) with a randomly chosen impact direction $`𝐧`$. In writing Eq. (2) we tacitly ignored the restriction $`𝐠𝐧>0`$ on the integration range, because the integrand obeys the reflection symmetry $`𝐧𝐧`$. The integration measure should be normalized, $`𝑑𝐧=1`$.

In the absence of energy input the system “cools” indefinitely according to Haff’s law . Indeed, the time dependence of the temperature is found from the Boltzmann equation (2) to give $`\frac{d}{dt}T=\lambda T^{3/2}`$ with $`\lambda =2ϵ(1ϵ)𝑑𝐧n_1^2`$ (the first axis was conveniently chosen to be parallel to $`𝐠`$). Since $`n_1^2+\mathrm{}+n_d^2=1`$ and $`𝑑𝐧=1`$ one has $`\lambda =2ϵ(1ϵ)/d`$. Thus, the temperature decays according to $`T(t)=T_0(1+t/t_0)^2`$, with the initial temperature $`T_0`$ and the characteristic time scale $`t_0=d/\left[ϵ(1ϵ)\sqrt{T_0}\right]`$. The temperature quantifies velocity fluctuations. In the following, we focus on the natural case of isotropic velocity distributions. We will show that asymptotically, the temperature represents the only relevant velocity scale as the velocity distribution approaches the scaling form

$$P(𝐯,t)\frac{1}{T^{d/2}}𝒫\left(\frac{v}{\sqrt{T}}\right).$$
(3)
Given the convolution structure of the Boltzmann equation (2), we introduce $`F(𝐤,t)`$, the Fourier transform of the velocity distribution function, $`F(𝐤,t)=𝑑𝐯e^{i𝐤𝐯}P(𝐯,t)`$. Applying the Fourier transform to Eq. (2) and integrating over the velocities gives

$$\frac{1}{\sqrt{T}}\frac{}{t}F(𝐤,t)+F(𝐤,t)=𝑑𝐧F[𝐤𝐪,t]F[𝐪,t],$$
(4)
with $`𝐪=(1ϵ)(𝐤𝐧)𝐧`$. This rate equation reflects the momentum transferred between particles during collisions. We seek an isotropic scaling solution for the Fourier transform, the equivalent of (3),

$$F(𝐤,t)=\mathrm{\Phi }\left(k^2T\right),$$
(5)
with $`k|𝐤|`$. In the $`k0`$ limit, the Fourier transform, $`F(𝐤,t)1\frac{1}{2}k^2T`$, implies that the first two terms in the Taylor expansion of the corresponding scaling function $`\mathrm{\Phi }(x)`$ are universal, $`\mathrm{\Phi }(x)1\frac{1}{2}x`$.

Let us first consider the simpler 1D case. Equation (4) reduces to $`\frac{1}{\sqrt{T}}\frac{}{t}F(k,t)+F(k,t)=F[ϵk,t]F[kϵk,t]`$ and the scaling function (5) satisfies

$$\lambda x\mathrm{\Phi }^{}(x)+\mathrm{\Phi }(x)=\mathrm{\Phi }\left[ϵ^2x\right]\mathrm{\Phi }\left[(1ϵ)^2x\right].$$
(6)
This equation admits a very simple solution

$$\mathrm{\Phi }(x)=\left[1+\sqrt{x}\right]e^\sqrt{x}.$$
(7)
The inverse Fourier transform gives the scaling function of the velocity distribution as a squared Lorentzian

$$𝒫(w)=\frac{2}{\pi }\frac{1}{(1+w^2)^2}.$$
(8)
Therefore, the form of the scaling solution (8) is universal as it is independent of the dissipation degree. Another important feature is the algebraic tail of the velocity distribution, $`𝒫(w)w^4`$ as $`w\mathrm{}`$.

We now return to the general $`d`$-dimensional case. Substituting (5) into (4) and using the temperature cooling equation $`\frac{d}{dt}T=\lambda T^{3/2}`$ one finds that the scaling function $`\mathrm{\Phi }(x)`$ satisfies $`\mathrm{\Phi }(x)\lambda x\mathrm{\Phi }^{}(x)=𝑑𝐧\mathrm{\Phi }(\xi x)\mathrm{\Phi }(\eta x)`$, where $`\xi =1(1ϵ^2)n_1^2`$ and $`\eta =(1ϵ)^2n_1^2`$. To integrate over the impact direction $`𝐧`$ we use spherical coordinates and treat the first axis as the polar axis, $`n_1=\mathrm{cos}\theta `$. The $`\theta `$-dependent factor of the measure is $`d𝐧=N^1(\mathrm{sin}\theta )^{d2}d\theta `$ with the factor $`N=_0^\pi (\mathrm{sin}\theta )^{d2}𝑑\theta =B(\frac{1}{2},\frac{d1}{2})`$ ensuring proper normalization ($`B(a,b)`$ is the beta function). Using $`\mu =\mathrm{cos}^2\theta `$ as the integration variable, the above governing equation for the scaling function reads

$$\lambda x\mathrm{\Phi }^{}(x)+\mathrm{\Phi }(x)=_0^1𝒟\mu \mathrm{\Phi }(\xi x)\mathrm{\Phi }(\eta x),$$
(9)
where $`\xi (\mu )=1(1ϵ^2)\mu `$ and $`\eta (\mu )=(1ϵ)^2\mu `$. Additionally, the integration measure was re-written using $`𝒟\mu `$, defined via $`B(\frac{1}{2},\frac{d1}{2})𝒟\mu =\mu ^{\frac{1}{2}}(1\mu )^{\frac{d3}{2}}d\mu `$ (it remains normalized, $`_0^1𝒟\mu =1`$). In the elastic case, the velocity distribution is Maxwellian. Indeed, $`\xi +\eta =1`$ and $`\lambda =0`$ when $`ϵ=0`$, and thence $`\mathrm{\Phi }(x)=e^{x/2}`$.

Our primarily goal is to determine statistics of extremely fast particles, namely the tail of the velocity distribution. This can be accomplished by noting that the large-$`v`$ behavior of the velocity distribution is reflected by the small-$`k`$ behavior of its Fourier transform. This is seen from the small-$`x`$ expansion of the 1D solution (7) which contains both regular and singular terms: $`\mathrm{\Phi }(x)=1\frac{1}{2}x+\frac{1}{3}x^{3/2}+\mathrm{}`$. The dominant singular $`x^{3/2}`$ term reflects the $`w^4`$ tail of $`𝒫(w)`$. In general, an algebraic tail of the velocity distribution (3),

$$𝒫(w)w^\sigma \mathrm{when}w\mathrm{},$$
(10)
indicates the existence of a singular component in the Fourier transform,

$$\mathrm{\Phi }_{\mathrm{sing}}(x)x^{(\sigma d)/2}\mathrm{when}x0,$$
(11)
and vice versa. To see this, it is convenient to take the Fourier transform $`\mathrm{\Phi }(x)_0^{\mathrm{}}𝑑ww^{d1}𝒫(w)e^{iw\sqrt{x}}`$, and write $`x=s^2`$ thus recasting it into the Laplace transform $`I(s)_0^{\mathrm{}}𝑑ww^{d1}𝒫(w)e^{ws}`$. The small-$`s`$ expansion of the integral $`I(s)`$ contains regular and singular components. For example, when $`\sigma <d`$, the integral $`I(s)`$ diverges as $`s0`$ and integration over large-$`w`$ yields the dominant contribution $`I_{\mathrm{sing}}(s)s^{\sigma d}`$. When $`d<\sigma <d+1`$, $`I(0)`$ is finite, but the next term is the above singular term, so $`I(s)=I(0)+I_{\mathrm{sing}}(s)+\mathrm{}`$. In general, the singular contribution is $`I_{\mathrm{sing}}(s)s^{\sigma d}`$, thereby leading to the singular term of Eq. (11). This is of course consistent with the one-dimensional behavior where $`\sigma =4`$.

The exponent $`\sigma `$ can now be obtained by inserting $`\mathrm{\Phi }(x)=\mathrm{\Phi }_{\mathrm{reg}}(x)+\mathrm{\Phi }_{\mathrm{sing}}(x)`$ into Eq. (9) and equating the dominant singular terms. Combining $`\mathrm{\Phi }_{\mathrm{reg}}(0)=1`$, with the anticipated leading singular term of Eq. (11), we find that the exponent $`\sigma `$ is a root of the following integral equation

$$1\lambda \frac{\sigma d}{2}=_0^1𝒟\mu \left[\xi ^{(\sigma d)/2}+\eta ^{(\sigma d)/2}\right].$$
(12)
This equation has a trivial solution $`\sigma =d+2`$, following from the identity $`𝒟\mu (\xi +\eta )=1\lambda `$, where the singular and the regular components simply coincide, $`x^{(\sigma d)/2}=x`$. Since we seek the leading singular term, the solution of Eq. (12) must therefore satisfy $`\sigma >d+2`$.

The integral equation (12) can be written explicitly in terms of special functions

$`1ϵ(1ϵ){\displaystyle \frac{\sigma d}{d}}=`$ (13)
$`{}_{2}{}^{}F_{1}^{}[{\displaystyle \frac{d\sigma }{2}},{\displaystyle \frac{1}{2}};{\displaystyle \frac{d}{2}};1ϵ^2]+(1ϵ)^{\sigma d}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\sigma d+1}{2}\right)\mathrm{\Gamma }\left(\frac{d}{2}\right)}{\mathrm{\Gamma }\left(\frac{\sigma }{2}\right)\mathrm{\Gamma }\left(\frac{1}{2}\right)}},`$ (14)
with $`{}_{2}{}^{}F_{1}^{}(a,b;c;z)`$ the hypergeometric function . Interestingly, the exponent $`\sigma \sigma (d,ϵ)`$ is a root of the transcendental equation (13) and thence it depends in a nontrivial fashion on spatial dimension $`d`$ and the restitution coefficient $`ϵ`$.

We first consider the dependence on the restitution coefficient by considering the quasi-elastic limit $`ϵ0`$. As discussed above, in the elastic case the velocity distribution is Maxwellian and the Fourier transform is simply $`\mathrm{\Phi }(x)=e^{x/2}`$ . This implies a diverging exponent $`\sigma \mathrm{}`$ as $`ϵ0`$. Therefore, the right-hand side of Eq. (13) vanishes, and the leading behavior is

$$\sigma \frac{d}{ϵ}\mathrm{as}ϵ0.$$
(15)
One can further expand $`\sigma (d,ϵ)`$ in the $`ϵ0`$ limit to find: $`\sigma (d,ϵ)=dϵ^1+a_1(d)ϵ^{1/2}+a_2(d)+\mathrm{}`$. We merely quote the leading correction in the physically relevant spatial dimensions $`a_1(2)=2(e^2+1)/\sqrt{\pi }`$ and $`a_1(3)=\sqrt{3\pi /2}`$. Clearly, the quasi-elastic limit is singular. Dissipation, even if minute, seriously changes the nature of the system .

Next, we discuss the dependence on the dimension. First, one can verify that $`\sigma =4`$ when $`d=1`$ using the identity $`{}_{2}{}^{}F_{1}^{}(a,b;b;z)=(1z)^a`$. The case $`d=1`$ is unique in that the entire scaling function and in particular the exponent become independent of the restitution coefficient. The case of $`d\mathrm{}`$ is similar to the $`ϵ=0`$ case in that the inelastic nature of the collisions becomes irrelevant, the velocity distribution is Maxwellian, and the exponent diverges, $`\sigma \mathrm{}`$ as $`d\mathrm{}`$. In this limit, the second integral in Eq. (12) is negligible as it vanishes exponentially with the dimension. The first integral can be evaluated by taking the limits $`d\mathrm{}`$ and $`\mu 0`$, with $`z=\mu d/2`$ being finite. Then, the integration measure is transformed $`𝒟\mu (\pi z)^{1/2}e^zdz`$, and Eq. (12) becomes $`1\lambda u=_0^{\mathrm{}}𝑑z(\pi z)^{1/2}e^{[1+(1ϵ^2)u]z}`$ where $`u=\frac{\sigma }{d}1`$. Performing the integration yields $`1\lambda u=[1+(1ϵ^2)u]^{1/2}`$, from which we find $`u`$ and

$$\sigma d\frac{1+\frac{3}{2}ϵϵ^3ϵ^{1/2}\left(1+\frac{5}{4}ϵ\right)^{1/2}}{ϵ(1ϵ^2)},$$
(16)
as $`d\mathrm{}`$. In general, $`\sigma d`$, and therefore, the algebraic decay becomes sharper as the dimension increases. The exponent $`\sigma (d,ϵ)`$ increases monotonically with increasing dimension, and additionally, it increases monotonically with decreasing $`ϵ`$ (see Fig. 1). Both features are intuitive as they mirror the monotonic dependence of the energy dissipation rate $`\lambda =2ϵ(1ϵ)/d`$ on $`d`$ and $`ϵ`$. Hence, the completely inelastic case provides a lower bound for the exponent, $`\sigma (d,ϵ)\sigma (d,ϵ=1/2)`$ with $`\sigma (d,1/2)=6.28753`$, $`8.32937`$, for $`d=2`$, $`3`$, respectively. The former value should be compared with $`\sigma 5`$, obtained from numerical simulations . The algebraic tails are characterized by unusually large exponents which may be difficult to measure accurately in practice (for typical granular particles $`ϵ0.1`$ yielding $`\sigma 30`$). Figure 1 also shows that the quantity $`\sigma /d`$ weakly depends upon the dimension, and the large-$`d`$ limit (16) provides a good approximation even for small dimensions.

Thus far, we discussed only freely cooling systems where the energy decreases indefinitely. In typical experimental situations, however, the system is supplied with energy to balance the energy dissipation . Theoretically, it is natural to consider white noise forcing , i.e., coupling to a thermal heat bath which leads to a nonequilibrium steady state. Specifically, we assume that in addition to changes due to collisions, velocities may also change due to an external forcing: $`\frac{dv_j}{dt}|_{\mathrm{heat}}=\xi _j`$ with $`j=1,\mathrm{},d`$. We use standard uncorrelated white noise $`\xi _i(t)\xi _j(t^{})=2D\delta _{ij}\delta (tt^{})`$ with a zero average $`\xi _j=0`$. The temperature rate equation is modified by the additional source term $`\frac{d}{dt}T+\lambda T^{3/2}=D`$, and the system approaches a steady state, $`T_{\mathrm{}}=(D/\lambda )^{2/3}`$. The relaxation toward this state is exponential, $`|T_{\mathrm{}}T|e^{t/\tau }`$.

Uncorrelated white noise forcing amounts to diffusion in velocity space, and Eq. (4) is modified as follows $`\frac{1}{\sqrt{T}}\frac{}{t}\frac{1}{\sqrt{T}}\frac{}{t}+Dk^2`$. In the steady state, the Fourier transform, $`F(𝐤,t=\mathrm{})\mathrm{\Psi }(y)`$ with $`y=Dk^2`$, obeys

$$(1+y)\mathrm{\Psi }(y)=\mathrm{\Psi }(\xi y)\mathrm{\Psi }(\eta y),$$
(17)
where integration with respect to the measure $`𝒟\mu `$ is denoted by $`f=_0^1𝒟\mu f(\mu )`$.

Equation (17) is solved recursively by employing the cumulant expansion $`\mathrm{\Psi }(y)=\mathrm{exp}\left[_{n1}(y)^nF_n\right]`$. Writing $`1+y=\mathrm{exp}\left[_{n1}(y)^n/n\right]`$, we recast Eq. (17) into

$`\mathrm{exp}\left[{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(y)^n}{n}}\right]=\mathrm{exp}\left[{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(y)^nG_n\right],`$ (18)
where $`G_n=F_n(1\xi ^n\eta ^n)`$. The desired cumulants $`F_n`$ are obtained by solving for $`G_n`$ recursively and then using $`F_n=G_n/1\xi ^n\eta ^n`$. In one dimension $`\mu ^n=1`$ and one immediately obtains $`G_n=n^1`$ . In higher dimensions, the averages acquire non-trivial dependence on $`n`$ though the qualitative nature of the solution remains the same since $`G_n`$ and hence the cumulants $`F_n`$ (as well as the moments of the distribution) are positive. This implies the following small-$`k`$ behavior of the Fourier transform at the steady state: $`\mathrm{ln}F_{\mathrm{}}(𝐤)(k/k_0)^2`$. Consequently, the large velocity tail of the velocity distribution is Gaussian

$$P_{\mathrm{}}(v)e^{(v/v_0)^2}\mathrm{when}vv_0.$$
(19)
This behavior agrees with the prediction of Ref.  derived via a small-$`ϵ`$ expansion, but differs from the stretched exponential behavior $`\mathrm{exp}(v^{3/2})`$ found for the driven inelastic hard sphere gas.

We have studied inelastic gases within the framework of the Boltzmann equation with a uniform collision kernel. In the freely evolving case, we have shown analytically that the density of high-energy particles is suppressed algebraically. The algebraic tails are characterized by remarkably large exponents, and may be hard to distinguish from (stretched) exponential tails. Our results, combined with previous kinetic theory studies which find exponential, stretched exponential, and Gaussian tails, indicate that the extremal characteristics can be very sensitive to parameters such as the restitution coefficient, and the dimension. These findings are all the more intriguing when compared with molecular dynamics simulations and scaling studies which suggest universal velocity distributions, at least in some cases . On the other hand, our results in the forced case support the near-Maxwellian assumptions typically used to obtain macroscopic transport coefficients from kinetic theory .

We thank A. Baldassari for communicating the results of Ref.  prior to publication, and for fruitful discussions. We also thank G. D. Doolen and S. Redner for their useful comments. This research was supported by DOE (W-7405-ENG-36) and NSF(DMR9978902).