text:
# Compton reflection and iron fluorescence in BeppoSAX observations of Seyfert type 1 galaxies

## 1 Introduction

Observations with the Ginga satellite in the 2–18 keV band of a sample of bright Active Galactic Nuclei (Nandra & Pounds 1994), following the discoveries with the same satellite by Piro et al. (1990), Matsuoka et al. (1990) and Pounds et al. (1990), definitely confirmed the presence of reprocessed, in addition to the direct, power law emission, as predicted by e.g. Guilbert & Rees (1988), Lightman & White (1988), George & Fabian (1991), Matt et al. (1991). The reprocessed radiation consists of two main components. The first is seen as a broad excess rising above about 10 keV, due to Compton reflection of the power law by optically thick, low ionization matter, the second is Kα fluorescence of mainly “cold” iron at 6.4 keV. Both appear to be common features in Seyfert galaxies of type 1 and intermediate between type 1 and 2 (objects of the latter type are often referred to as Narrow Emission Line Galaxies, NELG, their X-ray spectra are typically affected by a substantial, order 1022-23 H atoms cm-2, neutral absorbing column).

Once the inclination angle $`i`$ between line of sight and material is adopted, the amplitude of the reflection relative to the power law can be converted into a measure $`R=\mathrm{\Omega }/2\pi `$ of the solid angle subtended by the reflector. Through a reanalysis of Ginga data for 23 radio–quiet Seyfert 1 and NELG, Zdziarski et al. (1999) found evidence of a strong correlation between $`R`$ (estimated for a fiducial value of $`i=30^{}`$) and the photon slope $`\mathrm{\Gamma }`$ of the power law, which they tentatively describe with a function where the average $`R\mathrm{\Gamma }^v`$, $`v`$ about 12: $`R`$ rises from about 0.1 at $`\mathrm{\Gamma }`$=1.7 to about 1 at $`\mathrm{\Gamma }`$=2.1. They argue that the correlation is not an orientation effect, and, in the framework of the models where the production of the power law is associated with a very hot medium that comptonizes soft photons, they tentatively attribute it to a feedback of the additional photons from the reflection on the physical state of that medium. At present, however, its real significance, at least in its quantitative terms, is in doubt, because the existence of a strong correlation between these two variables in the fit of single observations is likely to bias substantially the result. Vaughan & Edelson (2001), for instance, on the basis of extensive simulations of data obtained with the PCA instrument onboard the RossiXTE satellite in an energy band similar to that of the Ginga data, have shown that a spurious correlation can be produced. Done & Nayakshin (2001) proposed yet an alternative possibility, namely that the correlation could be an artifact of the single–zone ionization (neutral in practice) model adopted for the reflector: however the simulations of Ginga spectra with their sophisticated model show that the correlation should be much shallower than the analytical description quoted above.

The equivalent width $`W_\alpha `$ of the iron fluorescence line, which inevitably arises in the reflecting material, is expected to be proportional to $`R`$. However, a line with $`W_\alpha `$ up to about 100 eV, that is a substantial fraction of the observed values, in an optimal geometrical arrangement can be produced also by marginally thin gas (e.g. Woźniak et al. 1998). Thanks to the superior energy resolution, the observations with the ASCA satellite of a representative sample of objects (Nandra et al. 1997) have shown that in most cases the line is resolved and broader than the velocity spread of the optical Broad Line Region (BLR). Moreover the first ASCA observation of MCG–6-30-15 (Tanaka et al. 1995) revealed a profile consistent with the special and general relativistic effects predicted for a line originating in the innermost regions of a flat accretion disk (Fabian et al. 1989). Although the line width in this object turned out to be rather extreme compared to the others in the Nandra et al. (1997) sample, it is widely believed (e.g. Fabian et al. 2000) that the fluorescent gas is most likely associated with the optically thick accretion disk, and therefore the broad line $`W_\alpha `$ should correlate with $`R`$. Unfortunately the ASCA band, being restricted to less than 10 keV, is inadequate to estimate $`R`$, and because of the variability which typically affects the continuum emission of these objects (in addition to reported cases of variations in the line, e.g. Weaver et al. 2001 and references in Fabian et al. 2000), the comparison of $`W_\alpha `$, as measured with ASCA, and of $`R`$ as measured with Ginga, would lack significance. An alternative approach chosen by Lubiński & Zdziarski (2001) aims at looking for a correlation between the broad line $`W_\alpha `$ and $`\mathrm{\Gamma }`$ in ASCA observations of 25 objects, which should correspond to the one between $`R`$ and $`\mathrm{\Gamma }`$. Their Fig. 1 shows indeed a trend, on average, in $`W_\alpha `$ (single gaussian fit) to increase with $`\mathrm{\Gamma }`$. This trend is quantified by fitting a relativistic diskline profile and an additional narrow gaussian profile to three average spectra, each one grouped according to the slope: the hardest group ($`\mathrm{\Gamma }_{\mathrm{av}}`$=1.62) has broad $`W_\alpha `$=72$`{}_{18}{}^{}{}_{}{}^{+17}`$ eV, the softest group ($`\mathrm{\Gamma }_{\mathrm{av}}`$=1.94) has broad $`W_\alpha `$=124$`\pm `$26, with an average inclination angle (a product of the diskline fit) practically the same for the three groups, $`i45^{}`$. The additional gaussian has an average narrow $`W_\alpha `$ of about 50 eV in all three groups. The authors interpret this component as originating in marginally thin gas: this is a possibility, as mentioned above, but the alternative one that it might instead be associated with highly reflecting material well beyond the accretion disk looks at present at least equally viable. For instance, the case of NGC 4051 observed with BeppoSAX in a state of exceptionally low intensity of its power law (Guainazzi et al. 1998) shows that the narrow line can be accompanied by the right amount of reflection, likely associated with a thick molecular torus well beyond the accretion disk, as predicted (e.g. Ghisellini et al. 1994) in Seyfert type 1 and 2 unified models (Antonucci 1993). The results in Lubiński & Zdziarski (2001) have been recently criticized by Yaqoob et al. (2002), and should be regarded with some caution. Besides this criticism, we note that the trend in their Fig. 1, mentioned above and which is based on individual measurements of $`W_\alpha `$ and $`\mathrm{\Gamma }`$, appears to be shallower than the one between $`R`$ and $`\mathrm{\Gamma }`$ claimed by Zdziarski et al. (1999).

In this context, the goal of this paper is to provide an assessement of the observational contribution of the italian–dutch satellite BeppoSAX (Boella et al. 1997a). This satellite is characterized by: a very broad energy band, 0.1–200 keV, covered by its Narrow Field Instruments (NFI), namely the LECS (0.1–10 keV, Parmar et al. 1997), the MECS (1.8–10.5 keV, Boella et al. 1997b) and the PDS (13–200 keV, Frontera et al. 1997); the energy resolution of the imaging ones (LECS and MECS) $`\mathrm{\Delta }E/E`$ 0.08($`E`$/6 keV)$`^{\frac{1}{2}}`$, which is intermediate between that of ASCA and of Ginga at the iron line; the very effective control of the background in the PDS. These properties have shown their value for the estimate in several Seyfert 1 galaxies of the spectral parameters of interest for this paper, provided their flux in the 2–10 keV band is at least (2–3)$`\times 10^{11}`$ erg cm-2 s-1 (e.g. Guainazzi et al 1999a, Perola et al. 1999).

Thanks to the signal detected in the PDS band, the strength of the reflection is basically measured in the fit by the excess counts in the 20–50 keV, where the broad maximum of this component is located, over the power law whose slope is constrained by the spectral counts in the LECS and MECS. Determining the slope of the power law is therefore the main issue, and the confidence beyond the purely statistical errors is related to the confidence that can be associated, object by object, to the treatment of additional spectral components, such as the cold and warm (ionized gas) absorbers and a generic soft excess. In addition, the power law itself may be affected by a high energy “break” (as first seen in NGC 4151, Jourdain et al. 1992), which we shall, as it is often done in the literature, empirically describe with an exponential factor, $`\mathrm{exp}(E/E_\mathrm{f})`$. Since the relative strength of the reflection below and at its broad maximum is very sensitive to the Compton recoil in the scattering of higher energy photons, hence to $`E_\mathrm{f}`$, the ability to estimate both $`R`$ and $`E_\mathrm{f}`$ thanks to the wide band and the efficient background control of the PDS (collimator used in rocking mode) is an asset of BeppoSAX. Nonetheless caution must be retained on single observation measurements of the three parameters $`\mathrm{\Gamma }`$, $`R`$ and $`E_\mathrm{f}`$, because systematic errors on $`\mathrm{\Gamma }`$, arising from an imperfect description of the spectrum at low energies, could be substantial compared to the purely statistical errors.

When the iron line is resolved, it might seem more appropriate to use models, both for the line and the associated reflection, which incorporates Doppler and relativistic effects. This approach would be physically more self-consistent, and in principle provide in addition the inclination $`i`$. However, as confirmed by observations with the Chandra and XMM-Newton satellites (detailed in Sect. 3 for our objects), the presence of a substantial narrow contribution is a rather common property, which is difficult to separate from the broad contribution with the BeppoSAX resolution, thus making this approach very tricky, particularly for the objects in our sample, where the line width is modest, $`\sigma 0.4`$ keV. We have therefore chosen to follow a conservative approach, where the model adopted for the reflection component is static, and the line is described as a single gaussian. We cared, though, to verify the outcome of adding to the model a narrow line at 6.4 keV: in general no improvement in the fit is achieved, moreover, when both a narrow and a broad component are found, their combined $`W_\alpha `$ does not differ significantly from that obtained with a single gaussian.

Despite their BeppoSAX observations fulfilled the selection criteria given in Sect. 2, in order to limit the incidence in our sample of the systematic errors mentioned above, we prudently excluded objects with a too complex and variable absorber, like NGC 4151 (Piro et al. 2001), NGC 3516 (Costantini et al. 2000, Guainazzi et al. 2001) and NGC 2992 (Gilli et al. 2000). The observation of MCG–6-30-15 (Guainazzi et al. 1999b) confirmed the ASCA finding (Tanaka et al. 1995) of a very broad and asimmetrical iron line, and was also excluded because it could not be treated homogeneously with the rest of the sample.

The paper is organized in the following manner. Sect. 2 describes the sample of the observations used. Sect. 3 describes the spectral fitting procedures and the results. Sect. 4 illustrates the distribution of the parameters in connection with the reprocessing. Sect. 5 is devoted to discussion and conclusions.

## 2 The sample data

The sample contains observations of objects with flux in the 2–10 keV band greater than 3$`\times 10^{11}`$ erg cm-2 s-1, duration longer than 50000 s in the MECS. For some objects, with separate observations shorter than 50000 s, after verifying with a hardness ratio analysis, as well as independent fits, that no spectral changes had occurred of such entity as to impair the reliability of the outcome relevant to the present paper, we merged two (NGC 5506, Mrk 509) or three observations (IC4329A(2)). Several of these observations are the subject of papers already published, whose reference will be given in Sect. 3.

Table 1 lists for each object: the observation start-dates and the good exposure times in the LECS, MECS and PDS, the redshift and the galactic hydrogen column, $`N_\mathrm{g}`$, adopted, the flux in the 2–10 keV band as observed, $`F_\mathrm{o}`$, corrected for any absorption as from the model fit, $`F_{\mathrm{na}}`$, and its corresponding luminosity ($`H_o`$=75 km s-1 Mpc-1).

The reduction procedure and screening criteria adopted are standard (see e.g. Guainazzi et al. 1999a). Concerning the PDS data, the Variable Rise Time (VRT) option was generally adopted, except for the observation of NGC 4593, where the Fixed Rise Time (FRT) was used. The PDS net spectra were obtained from direct subtraction of the off- from the on–source products. We verified that for none of the sample objects there was a significant contamination by hard sources in both the on- and the off–source field of view. Concerning the data from the LECS and the MECS images, the background was subtracted using spectra from blank sky files at the source position in the detectors.

The exposure times in the LECS are much shorter than in the MECS, because the former instrument is operated only during the night–time fraction of each orbit. Since the incident flux is generally recorded as variable during each observation, to take care of this inhomogeneity in the time coverage, in the fits the normalization of the LECS to the MECS was left as a free parameter. The normalization of the PDS to the MECS was instead held fixed at 0.8 (VRT selection) or 0.86 (FRT selection, Fiore et al. 1999). For the LECS data we used the version Sept. 1997 of the transfer matrix, and limited the range to 0.1–4.0 keV, where this matrix is more reliable. The energy bins chosen correspond to about one third of the instrument resolution, which itself is energy dependent.

## 3 Spectral analysis

The analysis of the spectral counts was performed using the software package XSPEC (Arnaud 1996, version 11.0.1). The Baseline Model Spectrum (BMS) includes: a photon power law (PL) with exponential cut–off, $`AE^\mathrm{\Gamma }\mathrm{exp}(E/E_\mathrm{f})`$, together with a reflection component (RC) from a cold slab isotropically illuminated by the PL photons, subtending the solid angle $`R=\mathrm{\Omega }/2\pi `$ and with inclination cos $`i`$ to the line of sight (Magdziarz & Zdziarski 1995, module pexrav); a uniform cold absorber column at the source, $`N_\mathrm{s}`$, in addition to $`N_\mathrm{g}`$ given in Table 1; a uniform warm absorber at the source, as described in the module absori, with column $`N_\mathrm{w}`$ and ionization parameter $`X_\mathrm{i}=L_\mathrm{i}/nD^2`$ erg cm s-1, where $`n`$ is the gas density and $`D`$ its distance from the source of ionizing luminosity $`L_\mathrm{i}`$ in the interval 5 eV to 20 keV (for $`L_\mathrm{i}`$, assumed with a power law spectrum, we set the slope equal to $`\mathrm{\Gamma }`$; the temperature of the absorbing gas was fixed as a rule at 3$`\times 10^4`$ K); a gaussian line to represent the iron Kα fluorescence, with energy $`E_\alpha `$, width $`\sigma _\alpha `$, intensity $`I_\alpha `$, accompanied by a gaussian of identical width, to represent the iron Kβ (with $`E_\beta `$ and $`I_\beta `$ set equal to 1.1$`E_\alpha `$ and 0.11$`I_\alpha `$, George & Fabian 1991). For both reflection and absorption we adopted the element abundances in Anders & Grevesse (1993). All fits were performed with cos $`i`$ fixed equal to 0.9 ($`i=26^{}`$). Alternatively cos $`i`$ was left free with $`R`$ fixed equal to 1. Because the dependence on cos $`i`$ of the RC shape is modest, all other parameters turned out identical in practice to those obtained with the angle fixed, and the differences in $`\chi ^2`$ are insignificant.

The results on the parameters of main interest for this paper are collected in Table 2. The line parameters are given in the source frame, and the line intensity $`I_\alpha `$ is also expressed as an equivalent width, $`W_\alpha `$. The errors represent the 90% confidence interval for two parameters of interest ($`\mathrm{\Delta }\chi ^2`$=4.61, Lampton et al. 1976), a realistic choice dictated by the statistical correlation which exists between couples such as ($`\mathrm{\Gamma }`$, $`R`$), ($`\mathrm{\Gamma }`$, $`E_\mathrm{f}`$), ($`\sigma _\alpha `$, $`I_\alpha `$).

To be noted that in general the value of $`E_\alpha `$ is perfectly consistent with 6.4 keV, the value for “cold” iron (in agreement with previous findings, e.g. Nandra & Pounds 1994, Nandra et al. 1997, Lubiński & Zdziarski 2001). In two objects (NGC 5506 and Mrk 509) the (per se marginal) evidence of a larger value finds explanation in recent XMM-Newton results, reported below.

Next we briefly describe the case of each source, including the best fit values of the other BMS parameters and additional spectral components adopted to improve the quality of the fit when the BMS gave an unacceptable $`\chi ^2`$. Generally $`X_\mathrm{i}`$ is very poorly constrained, hence we give only its best fit value. Reference is also given to some recent results on the iron line.

MCG 8-11-11. This observation was published by Perola et al. (2000), to which we refer for the choice of neglecting in the fit the LECS counts below 0.4 keV. $`N_\mathrm{s}`$ is negligible compared to $`N_\mathrm{g}`$, the warm absorber has $`N_\mathrm{w}`$=(6$`{}_{6}{}^{}{}_{}{}^{+9}`$)$`\times 10^{20}`$ cm-2 with $`X_\mathrm{i}`$7 erg cm s-1. The iron line is not resolved, the upper limit on its width is comparable to that of the lines resolved in our sample.

MCG–5-23-16. The cold absorber in this NELG is substantial, $`N_\mathrm{s}`$=(1.6$`\pm `$0.9)$`\times 10^{22}`$ cm-2; no signature is detected of a warm absorber. In ASCA (Weaver et al. 1997) and RossiXTE (Weaver et al. 1998) observations (comparable in flux to this one) a strong ($`W_\alpha `$250 eV) and broad iron line was consistently found. In the BeppoSAX observation the line is much narrower than $`\sigma _\alpha `$ 0.4 keV, the value obtained from a single gaussian description in both the ASCA and the RossiXTE data sets (Weaver et al. 1998); the line intensity is about twice that of the narrow gaussian, but about one third of the sum of the narrow and broad gaussians, as measured by Weaver et al. (1997) in their fit to the ASCA data. While this is probably yet another case of variations in the iron line (see also Smith & Done 1996 for a formally significant difference in $`W_\alpha `$ between two Ginga observations made a few days apart), we note that, contrary to expectation, the amount of reflection measured with RossiXTE (for a choice of $`E_\mathrm{f}`$=200 keV, the closest to our best fit value of the three adopted; $`\mathrm{\Gamma }`$ is identical to our best fit value; Weaver et al. 1998) is definitely not larger than in the BeppoSAX observation. Finally we note that a Chandra HETGS observation, briefly described in Weaver (2001), shows a narrow (FWHM less than 3000 Km s-1) line with $`W_\alpha `$ 90 eV, similar to our value.

NGC 3783. The BMS fit is poor ($`\chi ^2`$=208/150), with an excess in the residuals at low energies. With the addition of a line the $`\chi ^2`$ reduces to 160/147, with $`E_\mathrm{l}`$=0.59 keV and $`W_\mathrm{l}`$=140 eV: its strength is however almost ten times larger than that of the OVII emission lines detected by Kaspi et al. (2001) in a Chandra HETGS observation (when the source was fainter by only about 15%). A formally better improvement is achieved if, instead of a line, a black body is added, provided the temperature of the warm absorber is set equal to 105 K: the parameters given in Table 2 refer to this fit. $`N_\mathrm{s}`$ is negligible compared to $`N_\mathrm{g}`$, $`N_\mathrm{w}`$=(9.6$`{}_{1.4}{}^{}{}_{}{}^{+3.2}`$)$`\times 10^{21}`$ cm-2, with $`X_\mathrm{i}`$15 erg cm s-1, the black body $`kT`$=0.18$`{}_{0.03}{}^{}{}_{}{}^{+0.02}`$ keV. The iron line is definitely resolved. The Chandra HETGS observation (Kaspi et al. 2001) has revealed a narrow line at 6.4 keV, with an upper limit on the width ($`\sigma `$ less than 0.030 keV) that places its origin beyond the BLR: its intensity of (6.6$`\pm `$2.1)$`\times 10^5`$ cm-2 s-1 is about 55% of that measured by us. If we include in the fit a narrow, in addition to the broad line, no improvement is attained in $`\chi ^2`$ and the intensity of the narrow component is not well constrained, such that a value identical to that from the Chandra observation is acceptable. In a recent publication (De Rosa et al. 2002) on this BeppoSAX observation, the same model adopted by us (their Model E) yields somewhat different values, in particular $`\mathrm{\Gamma }`$=1.86$`\pm `$0.03, $`R`$=0.71$`{}_{0.28}{}^{}{}_{}{}^{+0.20}`$, $`W_\alpha `$=210$`\pm `$45, together with $`N_\mathrm{w}`$ about twice larger. The discrepancy is due to the different code adopted to describe the warm absorber, and it is not such as to influence the conclusions of the analysis following in Sect. 4. We note on the other hand the good agreement in our values of $`\mathrm{\Gamma }`$ and $`N_\mathrm{w}`$ with those obtained by Kaspi et al. (2001) in fitting the “line free zones” of the HETGS spectral data with a power law and a warm absorber.

NGC 4593. This observation was published by Guainazzi et al. (1999a). $`N_\mathrm{s}`$ is negligible compared to $`N_\mathrm{g}`$, and $`N_\mathrm{w}`$=(2.3$`{}_{0.9}{}^{}{}_{}{}^{+1.1}`$)$`\times 10^{21}`$ cm-2 with $`X_\mathrm{i}`$9 erg cm s-1. The iron line is resolved, $`W_\alpha `$ together with $`R`$ are the largest in our sample.

IC 4329A. Observation (1) was published by Perola et al. (1999). The average flux of observation (2) is 20% lower, yet there are no statistically significant differences between the two sets of parameters (the same holds for the three pointings merged, when fitted individually, despite a 25% difference between the lowest and the highest of their flux levels). The mean values of the absorbers are: $`N_\mathrm{s}`$=(2.7$`\pm `$1.0)$`\times 10^{21}`$ cm-2, $`N_\mathrm{w}`$=(2.8$`{}_{1.0}{}^{}{}_{}{}^{+1.2}`$)$`\times 10^{21}`$ cm-2 with $`X_\mathrm{i}`$5 erg cm s-1. The iron line is resolved. Given the strength of the source (the brightest of the sample) and the correspondingly good statistics, Perola et al. (1999) attempted a fit with a narrow gaussian plus a relativistic disk profile, and obtained marginal evidence for the existence of the narrow component. Another marginal detection of the narrow component in a simultaneous ASCA and RossiXTE observation is presented by Done et al. (2000). In a XMM-Newton observation with the EPIC instrument (Gondoin et al. 2001), at a flux level 25% higher than in observation (1), the narrow ($`\sigma `$ less than 0.06 keV) line is detected, with $`W_\alpha `$=43$`\pm `$1 eV; a broad component is seemingly absent, but unfortunately the authors do not provide an upper limit on its equivalent width. A simultaneous and relatively short (34000 s in the MECS) BeppoSAX observation is used, together with the EPIC data, by Gondoin et al. (2001) to estimate $`\mathrm{\Gamma }`$=1.93$`\pm `$0.03 and $`R`$=1.1$`\pm `$0.3: we believe that the estimate of $`R`$ is incorrect, because obtained with a PDS to MECS normalization factor equal to an improbable value of 0.7; by fitting the BeppoSAX data, as retrieved from the public archive, with this factor equal to 0.86 (for FRT selection, see Sect. 2), we obtain best fit values of $`R`$=0.54 and $`\mathrm{\Gamma }`$=1.91, fully consistent with those in Table 2.

NGC 5506. The two observations here merged have an almost identical average flux level, and when fitted individually yield similar values of the spectral parameters. The BMS fit with all three instruments gives an unacceptable $`\chi ^2`$=180/117, due to an excess in the residuals well below the strong photoelectric cut–off. Since the parameters of interest do not crucially depend on how this excess is modeled, we have chosen to present in Table 2 the results obtained excluding the LECS. This NELG has a substantial $`N_\mathrm{s}`$=(3.7$`\pm `$0.2)$`\times 10^{22}`$ cm-2, and no signatures of a warm absorber. The iron line is resolved and $`E_\alpha `$ is marginally inconsistent with 6.4 keV. This line is likely a blend of the two components recently revealed in a simultaneous XMM-Newton (with EPIC) and BeppoSAX observation (flux level similar to ours, Matt et al. 2001), one unresolved at 6.4 keV, the other resolved ($`\sigma `$ about 0.25 keV) at 6.75 keV, the two with comparable intensities whose sum is consistent with the one in our observation. Matt et al. (2001) present arguments, centred on the interpretation of the iron edge seen at 7.1 keV, that a substantial fraction of the RC in this object is very likely associated with the narrow line component. Evidence supporting this interpretation can be found in a spectral variability investigation conducted with RossiXTE by Lamer et al. (2000).

NGC 5548. This observation was published by Nicastro et al. (2000). They find internally evidence of a correlation between $`\mathrm{\Gamma }`$ and the flux (previously discovered in Ginga observations by Magdziarz et al. 1998, confirmed with RossiXTE observations by Chiang et al. 2000), which has some impact on the quality of our fit. In a Chandra LETGS observation Kaastra et al. (2002) find evidence of a very complex warm absorber, of a soft excess and of low energy emission lines. Following Kaastra & Barr (1989), who first reported a soft excess in EXOSAT observations, they model it as a “modified black body” with temperature $`kT`$0.1 keV. They suggest that the NVII emission line at 0.5 keV could be identified with the low energy line detected by Nicastro et al. (2000) in the BeppoSAX observation. We therefore added to the BMS a soft component, as described by the module DISKBB in XSPEC, with $`kT`$ fixed at 0.1 keV, and a line with $`E_\mathrm{l}`$ fixed at 0.5 keV. The fit yields $`W_\mathrm{l}`$=76$`{}_{54}{}^{}{}_{}{}^{+61}`$ eV, compatible with the Chandra result. $`N_\mathrm{s}`$ turns out negligible compared to $`N_\mathrm{g}`$, while $`N_\mathrm{w}`$=(2.4$`\pm `$0.7)$`\times 10^{21}`$ cm-2 (with $`X_\mathrm{i}`$14 erg cm s-1) is marginally consistent with the column density of the dominant component in the warm absorber as estimated by Kaastra et al. (2002). The iron line is unresolved and significantly narrower (less than 0.13 keV) than in any of the other sample objects. A Chandra HETGS observation (Yaqoob et al. 2001), with the source at a flux level about 30% lower, has revealed a narrow line at 6.4 keV with $`\sigma _\alpha `$=41$`{}_{24}{}^{}{}_{}{}^{+32}`$ eV, corresponding to a velocity spread consistent with an origin in the outer BLR (a substantial contribution to the line from thick matter beyond that region is not excluded). The intensity of this line matches the one measured by us. Using the same prescriptions for a relativistic disk line used by Yaqoob et al. (2001) to obtain an upper limit of 7$`\times 10^5`$ cm-2 s-1, we find a more stringent limit of 4$`\times 10^5`$ cm-2 s-1: this is a factor about 2 less than in earlier ASCA records of a substantially broader line, as reported in Yaqoob et al. (2001) or in Chiang et al. (2000). In a spectral variabilty study with RossiXTE, Chiang et al. (2000) found the intriguingly paradoxical (for the reprocessing paradigm) result that $`W_\alpha `$ is anticorrelated to $`R`$.

Mrk 509. The two observations and the results from their merging were published by Perola et al. (2000): in that paper it is shown that the BMS fit leaves the line width practically unconstrained (about 3 keV) and that in this respect the situation changes radically if a soft component is added in the form of a power law (preferred to a black body). This component has $`\mathrm{\Gamma }_\mathrm{s}`$=2.5$`{}_{0.2}{}^{}{}_{}{}^{+0.9}`$ and normalization 1.1$`\times 10^2`$ cm-2 s-1 keV-1 at 1 keV. $`N_\mathrm{s}`$ is negligible compared to $`N_\mathrm{g}`$, and $`N_\mathrm{w}`$=(4$`{}_{2}{}^{}{}_{}{}^{+8}`$)$`\times 10^{20}`$ cm-2 with $`X_\mathrm{i}`$16 erg cm s-1. The line width is poorly constrained and the line energy is marginally inconsistent with 6.4 keV. A XMM-Newton observation (with EPIC, Pounds et al. 2001) has revealed the presence of a resolved cold and of a broader ionized (6.9 keV) components. The $`W`$ of the former is about twice that of the latter, their sum is consistent with the value measured by us. The source, however, was almost 3 times fainter than in our observation, thus the blend must have experienced a change in intensity by a comparable factor.

NGC 7469. The BMS fit gives a marginally acceptable $`\chi ^2`$=191/151, but, like in the case of Mrk 509 as illustrated in Perola et al. (2000), the line width is practically unconstrained (about 2 keV). The addition of a soft excess in the form of a power law changes the situation radically. The slope of the soft power law is $`\mathrm{\Gamma }_\mathrm{s}`$=3.3$`\pm `$0.5, its normalization 5.5$`\times 10^3`$ cm-2 s-1 keV-1 at 1 keV. $`N_\mathrm{s}`$ is negligible compared to $`N_\mathrm{g}`$, and $`N_\mathrm{w}`$=(6.7$`{}_{3.7}{}^{}{}_{}{}^{+4.3}`$)$`\times 10^{20}`$ cm-2 with $`X_\mathrm{i}`$0.13 erg cm s-1. The line is not resolved, the upper limit on its width is comparable to that of the lines resolved in our sample. In a spectral variabilty study with RossiXTE, Nandra et al. (2000) detected changes in the line intensity which are not correlated with the absolute normalization of the reflection component, but they themselves doubt the reliability of their estimate of the second of the two parameters. We note that George et al. (1998) did not find evidence of a soft excess in their analysis (above 0.6 keV) of a ASCA observation with a flux 25% lower than in our data: thus the additional soft component, that we suggest here to provide a better description of the continuum of this observation, should be variable. We further note that the value of $`\mathrm{\Gamma }`$ in Table 2 is in good agreement with $`\mathrm{\Gamma }`$=1.92$`\pm `$0.02 obtained by Nandra et al. (2000) from the integrated fit of the data in the study mentioned above, with an average flux only 10% lower than in our observation.

## 4 An analysis of the parameters

This Section is devoted to illustrate the distribution of the parameters in the sample, in relation with the reprocessing of the PL.

### 4.1 $`W_\alpha `$ and R

Irrespective of model details, under the hypothesis that the line emission is entirely associated with optically thick material, to a good approximation $`W_\alpha `$ should correlate linearly with the value of $`R`$ obtained for an arbitrarily fixed inclination angle. Apart from the one introduced by the measurements errors, an intrinsic scatter in $`W_\alpha `$ about the correlation can be caused by a variance in the iron abundance (for differences within a factor 2 of the cosmic value, the dependence of $`W_\alpha `$ on the iron abundance is approximately linear, see Matt et al. 1997). Furthermore, for a given $`R`$, $`W_\alpha `$ may differ according to the value of $`\mathrm{\Gamma }`$, because it is proportional to the ratio of the photon number above 7.1 keV responsible for the fluorescence to the photon flux density at the line energy. The plot in Fig. 1 appears scatter dominated, and because of the rather large errors on both parameters little can be said about individual deviations from the hypothetical linear trend. The Spearman rank–order (Press et al. 1992) test, which does not take the errors into account, yields a correlation coefficient which is positive, $`r_\mathrm{s}`$=0.43, but with a large, 0.2, significance level for this quantity deviating from zero. Following another approach, where the error bars are also ignored, but which seems more useful in this context, we can compute for each object the ratio of the best fit value of $`W_\alpha `$ to the value expected, given the $`\mathrm{\Gamma }`$ of the PL and the best fit value of $`R`$ for cos $`i`$=0.9. We did that by using the $`\mathrm{\Gamma }`$ and angular dependence of $`W_\alpha `$ given in George & Fabian (1991), rescaled to the Anders & Grevesse (1993) abundances (see Matt et al. 1997). We obtain 0.97 for the mean value, and 0.29 for the standard deviation of the sample. A similar outcome (mean value 0.88, standard deviation 0.27) obtains when we proceed instead with R fixed equal to 1 and the best fit values of cos $`i`$ to estimate the expected $`W_\alpha `$. In our sample, then, on average the observed $`W_\alpha `$ is of the right amount expected from the simultaneously measured strength of the RC relative to the PL for the cosmic value of the element abundances. The standard deviation, being inevitably affected by the large errors on both $`W_\alpha `$ and $`R`$ can be considered a limit on the variance of the iron abundance in our sample.

Concerning $`R`$, we note the concentration (6 out of 9 objects) of best fit values in the interval 0.50–0.73, or the corresponding concentration of cos $`i`$ (for $`R`$=1) in the interval 0.35–0.58. The absence of objects with cos $`i`$ less than 0.35 (inclined more than about 70) can be understood in the frame of the Seyfert type 1 and type 2 unified models (Antonucci 1993), where the first type takes the appearence of the second at large inclinations. The paucity of objects with small, relative to intermediate inclinations could well be a statistical fluctuation, given the modest sample size, but it may reflect intrinsic differences in $`R`$, with objects where the intrinsic $`R`$ is indeed substantially less than unity.

### 4.2 $`R`$ and $`W_\alpha `$ versus $`\mathrm{\Gamma }`$

In a diagram of $`R`$ versus $`\mathrm{\Gamma }`$, scatter is expected due to differences in orientation. Fig. 2 shows scatter in $`R`$ only above $`\mathrm{\Gamma }`$=1.8, where the largest values are confined together with values comparable to those found below. While it could be due to a statistical fluctuation in the sampling of the angles, this mere fact provides some evidence of a trend in $`R`$ to increase on average with increasing $`\mathrm{\Gamma }`$. This is reflected in the Spearman rank–order correlation coefficient $`r_\mathrm{s}`$, which is positive and equal to 0.67, with a significance level for $`r_\mathrm{s}`$ deviating from zero equal to 0.03. The marginal significance of this evidence is further weakened by the consideration (see Sect. 1) that the two variables $`R`$ and $`\mathrm{\Gamma }`$ remain far from independent also in the fits to our very broad band data, as can be appreciated with a simple exercise: by fitting the data of each observation with $`\mathrm{\Gamma }`$ fixed at the two extremes of its error bar, the values of $`R`$ obtained differ typically by as much as 60 to 80% of their 90% confidence intervals. In any case, the trend in our sample (admittedly much smaller, 9 objects as opposed to 23) is definitely much shallower than the correlation found by Zdziarski et al. (1999; see also Sect. 4.3 for a supplementary cause likely to contribute to the difference).

In addition to the one due to differences in the orientation, a scatter due to differences in abundance is expected in a graph of $`W_\alpha `$ versus $`\mathrm{\Gamma }`$. A weak anticorrelation is predicted, with $`W_\alpha `$ decreasing with increasing $`\mathrm{\Gamma }`$, as mentioned above. The line (arbitrary normalization) describing this effect (deduced from George & Fabian 1991) is plotted in the graph of $`W_\alpha `$ versus $`\mathrm{\Gamma }`$ in Fig. 3. Notably $`W_\alpha `$ shows instead a weak tendency to increase, on average, with $`\mathrm{\Gamma }`$, which resembles that obtained by Lubiński & Zdziarski (2001, their Fig. 1). This is reflected in the sign of $`r_\mathrm{s}`$=0.57, which, contrary to the prediction, is positive. The significance level for $`r_\mathrm{s}`$ deviating from zero is however 0.08, three times larger than in the previous case. We note, though, that in the fits the variable $`W_\alpha `$ is less dependent on $`\mathrm{\Gamma }`$ than the variable $`R`$: the same exercise performed for $`R`$, and just mentioned, yields differences in $`W_\alpha `$ which are at most, and typically less than 35% of their 90% confidence interval. Thus we are inclined to regard the evidence of a shallow trend in $`W_\alpha `$ increasing on average with $`\mathrm{\Gamma }`$, albeit very marginal, no more so than it appears to be for $`R`$.

### 4.3 $`E_\mathrm{f}`$ versus $`\mathrm{\Gamma }`$

The plot of $`E_\mathrm{f}`$ versus $`\mathrm{\Gamma }`$ in Fig. 4 appears to confirm the existence of a correlation, with $`E_\mathrm{f}`$ increasing on average with $`\mathrm{\Gamma }`$, noted first in BeppoSAX data by Piro (1999) and then by Petrucci et al. (2001, their Fig. 6). Although the variable $`E_\mathrm{f}`$ in the fits, very much like $`R`$, is rather strongly dependent on $`\mathrm{\Gamma }`$, and therefore one should not disregard the possibility that our result might be biased by this effect, we note that in this case the correlation coefficient $`r_\mathrm{s}`$ is equal to 0.88, and the significance level for $`r_\mathrm{s}`$ deviating from zero is comfortably lower and equal to 0.0007. Discussing this correlation in terms of the mechanisms that give rise to the main spectral component goes beyond the scope of this paper (but see Petrucci et al. 2001). Our purpose here is to point out that, so long as one adopts the empirical representation of the direct continuum as a power law with an exponential cut–off, the existence of this type of correlation could, if ignored, artificially contribute to a correlation between $`R`$ and $`\mathrm{\Gamma }`$: in fact, as Ginga data do not allow to evaluate $`E_\mathrm{f}`$, by either neglecting the exponential cut–off or keeping it at a fixed value irrespective of $`\mathrm{\Gamma }`$ (Zdziarski et al. 1999 choose $`E_\mathrm{f}`$ = 400 keV fixed), one would systematically obtain lower values of $`R`$ for lower values of $`\mathrm{\Gamma }`$, for the reason already mentioned in Sect. 1. As an example of the interplay between $`E_\mathrm{f}`$ and $`R`$, we show the case of NGC 3783, whose $`\mathrm{\Gamma }`$ has an intermediate value in our sample: after excluding the cut–off we obtain $`\chi ^2`$=212.1/149 ($`\mathrm{\Delta }\chi ^2`$=61.7 demonstrates the high significance of the cut–off in the fit given in Table 2) and a best fit value of $`R`$=0.38 as opposed to $`R`$=0.63. Although the effect is sizeable, it alone is probably insufficient to explain the difference in steepness of the correlation found in the Ginga sample by Zdziarski et al. (1999) and in our sample, but it does represent an explanation of the offset between Ginga and BeppoSAX results noted, and tentatively attributed to calibration differences, by the same authors.

## 5 Discussion and conclusions

In confronting theory of reprocessing in AGN and observations, first comes a test of the predicted correlation between the strength of the reflection and that of the iron line fluorescence relative to the primary photons in the PL, from material mainly al low ionization, as indicated by the observed energy of the line. Our test, while improving on it, does confirm the result originally obtained by Nandra & Pounds (1994), who compared the mean sample values of $`R`$ and $`W_\alpha `$ with expectations from a face–on slab. We take properly into account the dependence of $`W_\alpha `$ on $`\mathrm{\Gamma }`$ in each observation and show that, be $`R`$ intrinsically equal to one or otherwise, the amount of fluorescence is close to what expected from the strength of the reflection, and that a modest variance in the iron abundance is likely present. In doing that, we assume that the fairly substantial “narrow” contribution to the line, which from Chandra and XMM-Newton observations appears to be a rather common feature in this type of AGN, is associated with optically thick gas. This assumption needs to be checked in the future, at present, to our knowledge, it is observationally confirmed in only one of the objects in our sample, NGC 5506.

We have then tested the existence of a correlation in $`R`$ and in $`W_\alpha `$ with respect to $`\mathrm{\Gamma }`$. In both cases we find a marginal evidence of a very gradual increase, on average, with increasing spectral slope. The trend in $`W_\alpha `$ is qualitatively similar to that found in the larger sample of ASCA observations (gaussian line fit) by Lubiński & Zdziarski (2001, their Fig. 1), while the one that we find in $`R`$ is much shallower than that found by Zdziarski et al. (1999) in Ginga data. Partly this discrepancy could be due to the existence of a correlation, of which we have found evidence in our sample, between the e–folding energy $`E_\mathrm{f}`$ (spanning from about 80 keV to more than 300 keV) and $`\mathrm{\Gamma }`$, which we took automatically into account in our fits, while Zdziarski et al. (1999) adopted a constant value equal to 400 keV. Thus our assessment, based on an admittedly much smaller sample of objects, is that the same phenomenological description of the spectral components, which are relevant to this issue, does not confirm the existence of a strong correlation between the relative strength of the Compton reflection, and the slope of the power law, while it provides some evidence that a very gradual correlation, and only on the average, for this quantity as well as for the iron line equivalent width, cannot be altogether excluded, but need to be confirmed with a sample of objects much larger than could be studied with BeppoSAX.

We conclude by commenting on three further limitations of this work. The first has to do with the narrow component of the line. On this matter further observational steps, which could not be pursued with BeppoSAX alone, are required, to disentangle the broad (accretion disk) component from the “narrow” component of the line, and simultaneously the fraction of the reflection associated with either of them. As recently shown by Malzac & Petrucci (2002) the nearly constant contribution by distant and thick matter, combined with a power law spectrum pivoting, as it varies in intensity, around energies in the range 1 to 10 keV, can introduce correlations between the reprocessing parameters and $`\mathrm{\Gamma }`$. Their attempt to explain in this scenario the paradoxical result obtained on NGC 5548 by Chiang et al. (2000), mentioned in Section 3, is particularly suggestive.

The second limitation concerns the ionization structure and radiative transfer properties of the reprocessing material. Possibly with the exception of a fraction of line emission in NGC 5506 and Mrk 509, in our sample the adoption of a single zone, low ionization structure of the reflector seems fully justified. However, Done & Nayakshin (2001) have developed a disk model with a sharp transition from outer regions of low ionization to inner regions with a thick skin, highly or fully ionized, where the reflected continuum would be hard to distinguish from the direct continuum. In this model obviously the reflected fraction pertaining to the outer regions would present a value of $`R`$ substantially less than unity. Moreover, their simulations of this model with the Ginga response, when fitted with the same model used by us, show that both $`R`$ and $`W_\alpha `$ turn out to correlate with $`\mathrm{\Gamma }`$ in a rather gentle manner, qualitatively similar to our results. It would be interesting to see this exercise repeated with the wider band BeppoSAX response.

The third limitation has to do with the shape of the spectrum of the primary component. The conventional exponential factor adopted is useful to reproduce the steepening of the PL at high energies, as required by the data, but clearly the value of $`R`$ which obtains is very sensitive to this shape. When using instead models of the continuum based on anisotropic comptonization in a disk corona (Haardt & Maraschi 1991, Haardt 1993) in fitting a sample of BeppoSAX observations, Petrucci et al. (2001) have shown that the value of $`R`$ turns out different, typically larger than that provided by a description like the one used by us.

Thus on the whole there is still a very ample margin of progress to be made, before the confrontation between theory and observations could be accomplished in a more satisfactory way.

###### Acknowledgements.

The BeppoSAX satellite is a joint Italian–Dutch program. We are grateful to the BeppoSAX Scientific Data Centre (now ASI–SDC) for assistance. GCP thanks Fabio La Franca for help on a problem in statistics. This work was partially supported by the Italian Space Agency (ASI), and by the Ministry for University and Research (MIUR) under grant COFIN–00–02–36. POP acknowledges a grant of the European Commission under contract ERBFMRX-CT98-0195.


format: text

text:
# Quantum Dynamics of Ultracold Na + Na2 collisions

## Abstract

Ultracold collisions between spin-polarized Na atoms and vibrationally excited Na2 molecules are investigated theoretically, using both an inelastic formalism (neglecting atom exchange channels) and a reactive formalism (including atom exchange). Calculations are carried out on both pairwise additive and non-additive potential energy surfaces for the quartet electronic state. In both inelastic and reactive calculations, the Wigner threshold laws are followed for energies below $`10^6`$ K. It is found that vibrational relaxation processes dominate elastic processes for temperatures below $`10^310^4`$ K. For temperatures below $`10^5`$ K, the rate coefficients for vibrational relaxation ($`v=10`$) from the full calculation are $`4.8\times 10^{11}`$ and $`5.2\times 10^{10}`$ cm3 s-1 for the additive and non-additive potentials respectively.

Methods for creating diatomic molecules in atomic Bose-Einstein condensates (BECs) are now starting to be realised experimentally. They include photoassociation spectroscopy  and magnetic tuning through a Feshbach resonance . An important long-range goal of such experiments is the production of a stable molecular BEC. However, in most cases the molecules are produced initially in vibrationally excited states, and their lifetime is limited by collisions with other atoms or molecules. Since the magnetic trap depth is typically 1 mK, any vibrationally or rotationally inelastic collision will release enough kinetic energy for both collision partners to be ejected from the trap.

Very little is known about vibrational relaxation rates for ultracold alkali dimers. Wynar et al.  produced 87Rb2 molecules in very high vibrational levels of the ground electronic state in an atomic BEC by stimulated Raman scattering. They analyzed their line shapes to obtain an upper bound on the inelastic rate coefficient $`k^{\mathrm{inel}}<8\times 10^{11}`$ cm3/s. In subsequent experiments  on ultracold (but not condensed) Rb2 molecules in a different set of vibrational levels, they measured $`k^{\mathrm{inel}}=3\times 10^{11}`$ cm3/s. These rates are too high for the production of long-lived molecules, but there is a hope that lower-lying states will relax more slowly  and that methods for stabilizing the molecules can be found.

Several calculations have been carried out on the vibrational relaxation of molecules such as H2  and CO  at ultralow energies. Spin-changing collisions of O2 have also been investigated.  However, collisions of alkali dimers present new theoretical challenges that are not present for collisions of stabler molecules. In particular, the potential energy surfaces are such that barrierless atom-exchange reactions can occur; even if the products are indistinguishable from the reactants, the reactive channels must be taken into account in a full treatment of the collision dynamics. Barrierless reactions are significantly different from reactions such as F + H2,  which have substantial barriers, and have not yet been investigated at ultralow energies. This Letter presents an initial investigation of alkali + alkali dimer collisions, for the case of Na + Na2 collisions occurring on the lowest quartet surface for Na3. This corresponds physically to collisions of atoms in their “stretched” spin states, with $`F=F_{\mathrm{max}}=I+S`$ and $`|M_F|=F`$.

In this work we use the potential energy surfaces of Higgins et al.  for the $`1^4A_2^{}`$ state of Na3. Their full (non-additive) potential was obtained from a grid of ab initio calculations with a large basis set, interpolated using the reproducing-kernel Hilbert space scheme of Ho and Rabitz . The potential has a global minimum at $`1222.1`$ K, with the atoms in an equilateral triangle configuration 4.41 Å apart. The symmetric linear geometry is a saddle point at $`554.4`$ K (667.7 K above the minimum), with $`r=5.10`$ Å. Comparison with the corresponding triplet Na2 pair potential , which has $`r_e=5.192`$ Å and $`D_e=255.7`$ K, shows that large non-additive effects are present. In the present work, we test the sensitivity of our results to the potential surface by performing calculations using both the full non-additive potential and a potential obtained by neglecting the non-additive part.

We initially carried out inelastic scattering calculations of vibrational relaxation from the $`v=1`$ state of triplet Na2, using the MOLSCAT program . Such calculations are carried out in Jacobi coordinates $`(R,r,\theta )`$ for a single arrangement of the atoms, and neglect the reactive channels. The wavefunction in the interaction region is expanded in a basis set of Na2 vibrational functions that are eigenfunctions of the Hamiltonian for free Na2, supplemented by a wall at large $`R`$ to provide a representation of the Na2 continuum.

For full close-coupling or coupled states calculations, it proved impossible to converge the basis set of Na2 vibrational functions. This may readily be understood in terms of the potential energy surface. Consider Na approaching Na2 for a T-shaped geometry, $`\theta =90^{}`$, allowing the Na2 bond length $`r`$ to relax to minimize the energy at each intermolecular distance $`R`$. The system passes through the minimum-energy equilateral geometry, and then the two atoms of the Na2 move apart to allow the Na atom to insert between them. However, even at $`R=0`$ (the linear geometry), the energy is below that for separated Na + Na2. At this configuration, the optimum value of $`r`$ is 10.4 Å, and the basis set of Na2 vibrational functions centered around $`r_e=5.1`$ Å is completely inadequate for representing such dramatically expanded geometries.

This problem does not arise for collinear scattering, for which the vibrational basis set is adequately converged for $`v_{\mathrm{max}}=20`$. We therefore carried out collinear calculations of the elastic and vibrational relaxation cross sections. The coupled equations were integrated from $`R=2.3`$ Å to 2500 Å. The resulting elastic and inelastic cross sections for collisions of Na2 initially in $`v=1`$ are shown in Figs. 1 and 2 for both the pairwise-additive and non-additive potentials.

The cross sections all follow the expected Wigner threshold laws  at energies below $`10^6`$ K. For elastic collisions, the cross sections are proportional to $`E^0`$, $`E^2`$ and $`E^3`$ for s, p and d waves respectively. The $`E^3`$ dependence for the d wave arises because the long-range dispersion ($`R^6`$) term in the atom-diatom potential modifies the threshold behavior for $`l>1`$ . For vibrational relaxation, the cross sections are proportional to $`E^{1/2}`$, $`E^{1/2}`$ and $`E^{3/2}`$ for s, p and d waves respectively. The oscillations that occur for $`E>1`$ mK arise simply from zeroes in $`\mathrm{sin}^2\eta `$, where $`\eta `$ is the scattering phase shift.

The elastic cross sections are dramatically different for the additive and non-additive potentials (by a factor of almost $`10^3`$), while the inelastic cross sections are much more similar (differing by a factor of only 10). The large difference in the elastic cross sections presumably occurs because the zero-energy elastic cross section may be written as $`\sigma _0^{\mathrm{el}}=4\pi |a|^2`$, where $`a`$ is the atom-diatom complex scattering length . This is accidentally close to zero for the non-additive potential: $`a=(0.901.19i)`$ nm, compared to $`(33.830.46i)`$ nm for the additive potential. This has the effect that $`\sigma ^{\mathrm{el}}>\sigma ^{\mathrm{inel}}`$ for $`T>10^{10}`$ K on the additive surface, but only for $`T>10^4`$ K on the non-additive surface.

An intriguing feature of the results is the appearance of a second linear region (in the log-log plots) for inelastic s-wave scattering between about $`10^5`$ and 10-2 K. The cross section in this region is proportional to $`E^1`$ for the non-additive potential but to about $`E^{4/3}`$ for the additive potential. A possible reason for this is that the vibrational coupling at long range decays as $`R^6`$ for the full potential but as $`R^8`$ for the pairwise-additive potential.

As described above, the MOLSCAT calculations are restricted to collinear geometries and a single arrangement channel. To go beyond this requires a reactive scattering formalism in which all three arrangement channels are included. Such calculations are more expensive, and have not previously been attempted for systems containing three non-hydrogen atoms.

As a first step towards a rigorous treatment, we have performed three-dimensional quantum reactive scattering calculations for total nuclear orbital angular momentum $`J=0`$. The configuration space is divided into an inner and an outer region, depending on the atom-diatom distance. In the inner region, we use a formalism based on body-frame democratic hyperspherical coordinates . This has already proved successful in describing atom-diatom insertion reactions such as N(2D) + H2$``$ NH + H  and O(1D) + H2$``$ OH + H . A related approach has also been used in studies of three-body recombination of ultracold atoms . The scattering wave function is expanded on a set of hyperspherical adiabatic states. This yields a set of close-coupling equations, which in our method are solved using the Johnson-Manolopoulos log-derivative propagator . In the outer region, we use the Arthurs-Dalgarno formalism , which is based on Jacobi coordinates. Matching of the wavefunctions in the inner and outer regions is performed on a boundary which is an hypersphere of radius 25 Å. This yields the reactance $`K`$-matrix and the scattering $`S`$-matrix.

The inner region starts at a hyperradius of 4 Å and is split into 297 sectors. The adiabatic states in each sector are obtained by a variational expansion on a basis of hyperspherical harmonics with $`A_1`$ symmetry. They are fully symmetric with respect to particle permutations to account for the indistinguishability of atoms. At large hyperradius, the adiabatic states concentrate into the arrangement channels and describe Na2 molecules in even $`j`$ states. At small hyperradius, they span a large fraction of configuration space and allow for atom exchange. The hyperspherical harmonic basis is truncated at $`\mathrm{\Lambda }_{\mathrm{max}}`$, the maximum value of the grand angular momentum. $`\mathrm{\Lambda }_{\mathrm{max}}`$ varies from 198 (867 harmonics) at small hyperradius to 398 (3400 harmonics) at large hyperradius. A fixed number of 135 adiabatic states is used in the close-coupling expansion in each sector. At the boundary between the inner and outer regions, the adiabatic states are projected onto a set of Na2 rovibrational states, with $`j_{\mathrm{max}}=48`$, 44, 40, 36, 30, 26, 20, 10 for vibrational levels $`v=0,1,\mathrm{},7`$. The boundary between the inner and outer regions was placed at a distance such that couplings due to the atom-diatom residual interaction can be neglected outside the boundary. In the outer region, regular and irregular solutions of a radial Schrödinger equation which includes the isotropic ($`R^6`$) part of the interaction were integrated inwards from very large distances (5000 Å).

Fig. 3 shows the cross sections as a function of collision energy for the additive and non-additive potentials. All collisions that are energetically elastic, with or without atom exchange, are included in the elastic cross section. All other processes (which produce Na2 $`(v=0,j)`$) contribute to the quenching cross section. Na2 rotational levels up to $`j=20`$ are energetically accessible at the energy of the $`v=1`$ state (23.5 cm-1), and all accessible levels are populated in the products.

It may be seen that the Wigner threshold laws for the elastic and quenching cross sections hold below $`10^5`$ K for both potentials. The cross sections are larger for the non-additive than for the additive potential, by a factor of about 10 for both elastic and quenching cross sections. The $`E^{1/2}`$ dependence of the quenching cross sections corresponds to a constant rate coefficient below $`10^5`$ K, which is $`k^{\mathrm{inel}}=4.8\times 10^{11}`$ cm3 s-1 for the additive potential and $`5.2\times 10^{10}`$ cm3 s-1 for the non-additive potential. The corresponding scattering lengths are $`a=(2.260.80i)`$ and $`(2.388.54i)`$ nm, respectively.

Outside the Wigner region, the cross sections have a more complicated energy dependence. The quenching probability increases with increasing energy and approaches unity at the limit of the Wigner region. The quenching cross sections thus vary approximately as $`E^1`$ above $`10^4`$ K, because of the $`k^2`$ factor in the expression for the cross section. Above $`10^3`$ K, the elastic cross sections show oscillations which are similar to but less pronounced than the ones in Fig. 1.

Finally, the ratio of quenching to elastic cross sections is larger than 1 at energies below $`10^4`$ K for the additive potential and below $`10^3`$ K for the non-additive potential. It increases up to 500 in the nK range for both potentials.

In future work, we intend to investigate the dependence of quenching rates on the initial vibrational quantum number and to investigate the effects of magnetic fields and nuclear spin coupling.

The three-dimensional quantum dynamical calculations were performed on a NEC-SX5 vector supercomputer, through a grant from the “Institut du Développement des Ressources en Informatique Scientifique” (IDRIS) in Orsay (France). PS is grateful to the EPSRC for a Research Associateship under grant no. GR/R17522/01.


text:
# 1 Introduction

## 1 Introduction

There are two main motivations to study a deformed many-body system. The first one comes from fundamental physics. It is well-known that the noncommutative geometry plays a very important role in unravelling the properties of the Planck scale physics. It has for a long time been suspected that the noncommutative spacetime might be a realistic picture of how spacetime behaves near the Planck scale. In fact, the noncommutative geometry naturally enters the theory of open string in a background $`B`$-field. In particular, the noncommutative geometry makes the holography (e.g. the AdS/CFT correspondence) of a higher dimensional quantum system of gravity and lower dimensional theory possible. A kind of special regularization with exponentially increasing spacetime cutoff has been introduced based on the noncommutative geometry. A very small displacement of the noncommutative deformation parameter from its classical value reduces sharply the entropy of quantum system of gravity. An adequately adopted noncommutative deformation of geometry makes the holography of a higher dimensional quantum system of gravity and lower dimensional theory possible. It was also discovered that simple limits of $`M`$ theory and superstring theory lead directly to the noncommutative gauge field theory. The noncommutative field theory has been intensively studied in the past two decads. The most realistic laboratory for testing the Planck scale physics is the early universe and the black hole physics. From the standard picture of cosmology, the universe started from a very hot, dense phase about 15 billion years ago in the big bang. The very early universe was opaque due to the constant interchange of energy between matter and radiation. Several problems of standard big bang cosmology, such as the horizon, homogeneity, entropy and flatness problems, can be explained under the assumption that the very early universe underwent a period of extremely rapid inflation, driven by the potential of some assumed inflation fields. In particular, the inflation scenario is also able to explain the observed perturbations in the cosmic microwave background radiation, namely as originating ultimately from quantum fluctuations of the inflation field. However, in most of the current models of inflation, the physical length of perturbations was much smaller than the Planck length at the beginning of inflation. It is then of interest to investigate whether the predictions of the spectrum of cosmological perturbations are sensitive to the unknown super-Planck-scale physics-. Recent studies- by physicists and astrophysicists show that without strange inflation fields but just the noncommutative geometry and the deformed radiation dominated universe subject to the noncommutative spacetime quantization indeed give rise to inflation as the radiation temperature exceeds the Planck temperature. Based on the noncommutative geometry and the deformed radiation, one can even get an alternative explanation for the horizon, homogeneity, entropy and flatness problems in the standard cosmology.

Another motivation comes from quantum optics as well as molecular and atomic spectroscopy. A fascinating development of quantum optics in last two decades has taken place within two specific areas, one concerning the study of nonclassical states of light, and the other studying the interaction of light with matter. The nonclassical states of light offer the possibility of the so-called reduction of quantum noise below the standard limit, determined by the Heisenberg’s uncertainty relation. Study on interaction of light with matter shows that the old Bohr quantum jumps are real. Dynamics of deformed photonic field interacting with matter presented many interesting properties both in the two important areas-. Great success has also been made by making use of a deformed oscillator to study molecular spectrum-.

In this paper, we study thermodynamics of a deformed Bose gas, and in particular, of a deformed radiation field. On the topic, many papers have been published in more than ten years-. However, we still have not a general form of statistical distribution function of the deformed Bose gas, which should reduce to the Bose-Einstein distribution function for an ideal Bose gas in the $`q=1`$ limit. In this situation, it is difficult to investigate systematically the influence of the deformation on thermodynamic characteristics of a Bose system, which is interesting in fundamental physics and other subjects as discussed above. For this purpose, by making use of the double-time Green function technique, we obtain an algebra equation, which should be satisfied by the distribution function of a deformed Bose gas. The equation for distribution function can not be solved exactly in general. The situation is similar to the case of magnetic lattices. From this point of view, a deformed Bose gas is really a strongly correlated system. In the case of $`q1`$, we give a useful expression of the distribution function for a deformed Bose gas. Thermodynamic properties, such as internal energy, specific heat, free energy and equation of state, are discussed in detail.

The paper is organized as follows. In section 2, a deformed Bose gas is discussed. We get the distribution function for a deformed Bose gas by using the double-time Green function technique in section 3. In section 4, we calculate some thermodynamic quantities for a deformed photon gas and compare them with those for the ideal Bose gas. Concluding remarks are given in section 5.

## 2 A deformed Bose gas system

The Lagrange density of a scalar field $`\varphi (x)`$ with mass $`m`$ reads

$$(x)=\frac{\mathrm{}^2}{2}\frac{\varphi }{x_\mu }\frac{\varphi }{x^\mu }\frac{1}{2}m^2c^2\varphi ^2.$$
(1)
The Euler-Lagrange equation immediately leads to the Klein-Gordon equation

$$\left(\mathrm{}^2\frac{}{x_\mu }\frac{}{x^\mu }+m^2c^2\right)\varphi (x)=0.$$
(2)
The canonically conjugate field is

$$\pi (x)=\frac{}{\dot{\varphi }(x)},$$
(3)
which leads to the Hamilton density

$$(x)=\frac{1}{2}\left(\frac{1}{\mathrm{}^2}\pi ^2(x)+(\mathrm{}\varphi (x))^2+m^2c^2\varphi ^2(x)\right).$$
(4)
The Fourier decomposition of the scalar field is of the form,

$$\begin{array}{c}\varphi (x)=\frac{d^3kc}{\sqrt{(2\pi )^32\mathrm{}\omega _k}}\left(a(k)e^{ikx}+a^{}(k)e^{ikx}\right),\hfill \\ \pi (x)=\frac{d^3ki\mathrm{}^2\omega _k}{\sqrt{(2\pi )^32\mathrm{}\omega _k}}\left(a(k)e^{ikx}+a^{}(k)e^{ikx}\right),\hfill \end{array}$$
(5)
where we have used the notation $`\omega _k=\frac{c}{\mathrm{}}\sqrt{(\mathrm{}𝐤)^2+m^2c^2}`$. The canonical quantization of the scalar field reads

$$[\widehat{\varphi }(𝐱,t),\widehat{\pi }(𝐲,t)]=i\mathrm{}c\delta ^3(𝐱𝐲).$$
(6)
In terms of the Fourier modes $`\widehat{a}(k)`$ and $`\widehat{a}^{}(k)`$, the commutation relations can be rewritten into the form

$$[\widehat{a}(k),\widehat{a}^{}(k^{})]=\delta ^3(𝐤𝐤^{}).$$
(7)
It is not difficult to express the Hamiltonian in terms of these Fourier modes,

$$\widehat{H}=\frac{1}{2}d^3k\mathrm{}\omega _k\left(\widehat{a}(k)\widehat{a}^{}(k)+\widehat{a}^{}(k)\widehat{a}(k)\right).$$
(8)
In the quantum field theory, the vacuum state is defined as follows

$$\widehat{a}(k)|0=0.$$
(9)
And the $`N`$-particle Fock space can be written down directly as

$$|n(k_1),n(k_2),\mathrm{},n(k_m)=\frac{\left(\widehat{a}^{}(k_1)\right)^{n(k_1)}}{\sqrt{n(k_1)!}}\frac{\left(\widehat{a}^{}(k_2)\right)^{n(k_2)}}{\sqrt{n(k_2)!}}\mathrm{}\frac{\left(\widehat{a}^{}(k_m)\right)^{n(k_m)}}{\sqrt{n(k_m)!}}|0.$$
(10)
To measure how many particles of a certain momentum exist in the field system, one can introduce the particle number operator $`\widehat{N}(k)`$,

$$\widehat{N}(k)=\widehat{a}^{}(k)\widehat{a}(k).$$
(11)
It is easy to check that

$$\widehat{N}(k_i)|n(k_1),n(k_2),\mathrm{},n(k_m)=n(k_i)|n(k_1),n(k_2),\mathrm{},n(k_m).$$
(12)
Thus, a field can be seen as a many-body system that consists of free harmonic oscillators with different momentum. Properties of the field system can be described exactly by the Weyl-Heisenberg algebra,

$$\begin{array}{c}[\widehat{a},\widehat{a}^{}]=1,\hfill \\ [\widehat{N},\widehat{a}]=\widehat{a},[\widehat{N},\widehat{a}^{}]=\widehat{a}^{},\hfill \end{array}$$
(13)
with the Fock space $`|n=\frac{(\widehat{a}^{})^n}{\sqrt{n!}}|0`$ as its representation. All properties , in particular, the thermodynamic properties of the many-body system can be investigated systematically, by making use of the Hamiltonian (8) and the Weyl-Heisenberg algebra (13).

The natural quantum of the quantum field theory on noncommutative geometry is the so-called deformed harmonic oscillator. On the topic of deformed oscillators, a lot of papers have been published-. Similar to the case of an ordinary quantum field system, the many-body properties of a noncommutative quantum field system can be described by a collection of the deformed harmonic oscillators,

$$\widehat{H}_q=\frac{1}{2}\underset{𝐤}{}\mathrm{}\omega _k\left(\widehat{a}_q(k)\widehat{a}_q^{}(k)+\widehat{a}_q^{}(k)\widehat{a}_q(k)\right).$$
(14)
The deformed harmonic oscillator is related with the simple harmonic oscillator as follows,

$$\widehat{a}_q=\sqrt{\frac{[\widehat{N}]_q}{\widehat{N}}}\widehat{a},\widehat{a}_q^{}=\widehat{a}^{}\sqrt{\frac{[\widehat{N}]_q}{\widehat{N}}},$$
(15)
where we have used the notation $`[x]_q={\displaystyle \frac{q^xq^x}{qq^1}}`$.
The deformed Weyl-Heisenberg algebra reads

$$\begin{array}{c}\widehat{a}_q\widehat{a}_q^{}q\widehat{a}_q^{}\widehat{a}_q=q^{\widehat{N}},\hfill \\ [\widehat{N},\widehat{a}_q^{}]=\widehat{a}_q^{},[\widehat{N},\widehat{a}_q]=\widehat{a}_q.\hfill \end{array}$$
(16)
The representation of the deformed Weyl-Heisenberg algebra is obtained by constructing the Fock space,

$$\begin{array}{c}|n_q=\frac{(\widehat{a}_q^{})^n}{\sqrt{[n]_q!}}|0,\hfill \\ \widehat{a}_q|n_q=\sqrt{[n]_q}|n1_q,\widehat{a}_q^{}|n_q=\sqrt{[n+1]_q}|n+1_q,\hfill \end{array}$$
(17)
where $`[n]_q![n]_q[n1]_q\mathrm{}[2]_q[1]_q`$. It is easy to check that one still can measure how many particles of a certain momentum exist in the deformed field system by making use of the particle number operator $`\widehat{N}(k)`$,

$$\widehat{N}(k_i)|n(k_1),n(k_2),\mathrm{},n(k_m)_q=n(k_i)|n(k_1),n(k_2),\mathrm{},n(k_m)_q.$$
(18)
The Hamiltonian of a deformed harmonic oscillator- can be written as

$$\widehat{H}_q=\frac{\mathrm{}\omega }{2}([\widehat{N}]_q+[\widehat{N}+1]_q).$$
(19)
Below we will calculate the distribution function of a deformed harmonic oscillator system by making use of the double-time Green function technique.

## 3 Distribution function

The double-time Green function (retarded or advanced) for a Bose system is generally defined as

$$G_\rho (tt^{})=\frac{i}{2}\left((\rho +1)\theta (tt^{})+(\rho 1)\theta (t^{}t)\right)[\widehat{A}(t),\widehat{B}(t^{})],$$
(20)
where we have used the notations,

$$\begin{array}{c}\widehat{A}(t)\mathrm{exp}(\frac{iHt}{\mathrm{}})\widehat{A}\mathrm{exp}(\frac{iHt}{\mathrm{}}),\hfill \\ \widehat{O}=Z^1\mathrm{Tr}(e^{\beta H}\widehat{O}),Z=\mathrm{Tr}(e^{\beta H}),\hfill \end{array}$$

here $`\rho =1`$ for the retarded one and $`\rho =1`$ for the advanced one.

In Eq.(20), setting

$$\widehat{A}(t)\widehat{a}(t),\widehat{B}(t^{})\widehat{a}^{}(t^{}),\rho =1,$$

(the retarded Green function is assumed), we get

$$G_R(tt^{})=i\theta (tt^{})\widehat{a}(t)\widehat{a}^{}(t^{})\widehat{a}^{}(t^{})\widehat{a}(t).$$
(21)
The equation of motion of the retarded Green function reads

$$i\mathrm{}\frac{}{t}G_R(tt^{})=\mathrm{}\left(\delta (tt^{})\widehat{a}(t)\widehat{a}^{}(t^{})\widehat{a}^{}(t^{})\widehat{a}(t)+\theta (tt^{})\frac{d\widehat{a}(t)}{dt}\widehat{a}^{}(t^{})\widehat{a}^{}(t^{})\frac{d\widehat{a}(t)}{dt}\right).$$
(22)
The Heisenberg equation for the operator $`\widehat{a}(t)`$

$$i\mathrm{}\frac{d\widehat{a}(t)}{dt}=[\widehat{a}(t),\widehat{H}(t)]$$

gives,

$$\frac{d\widehat{a}(t)}{dt}=i\frac{\omega }{2}\left(q^{\widehat{N}+1}+q^{(\widehat{N}+1)}\right)\widehat{a}(t).$$
(23)
Then, we have

$$i\mathrm{}\frac{}{t}G_R(tt^{})=\mathrm{}\delta (tt^{})+\frac{\mathrm{}\omega }{2}\left(q^{\widehat{N}+1}+q^{\widehat{N}+1}\right)G_R(tt^{}).$$
(24)
Making use of the Fourier transformation

$$G_R(tt^{})=_{\mathrm{}}^{\mathrm{}}\frac{dE}{2\pi \mathrm{}}G_R(E)e^{iE(tt^{})/\mathrm{}},$$
(25)
we get

$$G_R(E)=\frac{\mathrm{}}{E\frac{\mathrm{}\omega (q^{\widehat{N}+1}+q^{\widehat{N}+1})}{2}}.$$
(26)
The spectrum theorem of Green functions tells us

$$\widehat{a}^{}(t^{})\widehat{a}(t)=i_{\mathrm{}}^{\mathrm{}}𝑑E\frac{(G_R(E+i0^{})G_R(Ei0^{}))e^{iE(tt^{})/\mathrm{}}}{2\pi \mathrm{}(e^{\beta E}1)}.$$
(27)
Then, we can write formally the distribution function of a deformed Bose gas as

$$n(\epsilon )=\widehat{N}=\frac{1}{\mathrm{exp}(\frac{\mathrm{}\omega \beta (q^{\widehat{N}+1}+q^{\widehat{N}+1})}{2})1}.$$
(28)
Evidently, in the case of $`q1`$, $`n(\epsilon )`$ recovers the familiar Bose-Einstein distribution function. In principle, all interesting properties of a deformed Bose gas can be obtained by solving this equation for distribution function. This is to be one of the most advantageous features of the Green function approach.

Now we try to solve Eq.(28) to gain an insight into the influence of the deformation on the distribution function and to compare the deformed distribution function with the Bose-Einstein distribution function.

Expanding the right hand side of Eq.(28) in the neighborhood of $`q=1`$ to the second order, we get an exactly solvable equation for the deformed distribution function. The solution is of the form,

$$\begin{array}{ccc}\hfill n(\epsilon )& =& \frac{1}{\delta ^2\beta \epsilon e^{\beta \epsilon }}\left(\delta ^2\beta \epsilon e^{\beta \epsilon }+\left(e^{\beta \epsilon }1\right)^2\right)\hfill \\ & & +\frac{1}{\delta ^2\beta \epsilon e^{\beta \epsilon }}\sqrt{\left(\delta ^2\beta \epsilon e^{\beta \epsilon }+\left(e^{\beta \epsilon }1\right)^2\right)^2\delta ^2\beta \epsilon e^{\beta \epsilon }\left(\delta ^2\beta \epsilon e^{\beta \epsilon }2(e^{\beta \epsilon }1)\right)},\hfill \end{array}$$
(29)
where we have used the notation $`\delta q1`$ for short.

We note that it is difficult to compare the distribution function of the deformed Bose gas, Eq.(29), with the Bose-Einstein distribution function directly. For the aim, we rewrite Eq.(29) as follows,

$$n(\epsilon )=\frac{(e^{\beta \epsilon }1)2^1\delta ^2\beta \epsilon e^{\beta \epsilon }}{(e^{\beta \epsilon }1)^2+\delta ^2\beta \epsilon e^{\beta \epsilon }}\frac{\delta ^2\beta \epsilon e^{\beta \epsilon }\left(\delta ^2\beta \epsilon e^{\beta \epsilon }2\left(e^{\beta \epsilon }1\right)\right)^2}{\left(2\delta ^2\beta \epsilon e^{\beta \epsilon }+2\left(e^{\beta \epsilon }1\right)^2\right)^3}.$$
(30)
One can further expand the distribution function into a series over deformation parameter $`\delta `$,

$$n(\epsilon )=\frac{1}{e^{\beta \epsilon }1}\delta ^2\left(\frac{\beta \epsilon e^{\beta \epsilon }}{2(e^{\beta \epsilon }1)^2}+\frac{\beta \epsilon e^{\beta \epsilon }}{(e^{\beta \epsilon }1)^3}+\frac{\beta \epsilon e^{\beta \epsilon }}{2(e^{\beta \epsilon }1)^4}\right).$$
(31)
From the above form of the distribution function of the deformed Bose gas, we see clearly that, at the case of $`q=1`$, the Bose-Einstein distribution is recovered. To further investigate the effect of deformation on the distribution function, we present a plot of the distribution function in Fig.1.

FIG.1 The distribution function of a deformed Bose gas

( start from up $`q=1.000,1.001,1.002,1.003,1.004,1.005,1.006,1.010`$, here we have used the notation $`\beta \epsilon =E`$.)

To describe behaviors of the distribution function of the deformed Bose gas near the ground state more accurately, we give another plot of the deformed distribution function in Fig.2.

FIG. 2 The behaviors of distribution function

of a deformed Bose gas near the ground state

(start from up $`q=1.000,1.001,1.002,1.003,1.004,1.005,1.006,1.010`$, here we have used the notation $`\beta \epsilon =E`$.)

## 4 Thermodynamic properties

In this section, we concentrate on the study of thermodynamic properties of a deformed photonic field. As being well-known, the dispersion relation of a photon is

$$ϵ=cp.$$
(32)
We know that a photonic field can be described by a collection of simple harmonic oscillators and have relation

$$\widehat{H}=\mathrm{}\omega \left(\widehat{N}+\frac{1}{2}\right).$$
(33)
One can then read out that

$$p=\frac{\mathrm{}\omega }{c}\left(n+\frac{1}{2}\right).$$
(34)
This relation and the Hamiltonian (19) give the dispersion relation of a deformed Bose gas,

$$\epsilon =\frac{1}{2}\mathrm{}\omega \left(\left[\frac{c}{\mathrm{}\omega }p+\frac{1}{2}\right]_q+\left[\frac{c}{\mathrm{}\omega }p\frac{1}{2}\right]_q\right).$$
(35)
Expanding the right hand side of (35) in the neighborhood of $`q=1`$ to the second order, we get

$$\epsilon =cp\frac{\delta ^2}{24}cp+\frac{\delta ^2}{6\mathrm{}^2\omega ^2}c^3p^3.$$
(36)
The density of energy levels is as follows,

$$\begin{array}{ccc}\hfill a(\epsilon )& =& \frac{8\pi V}{3h^3c^3}\frac{3\sqrt{3B}\epsilon +\sqrt{4A^3+27B\epsilon ^2}}{4B\sqrt{4A^3+27B\epsilon ^2}}\frac{\left(2\times 6^{1/3}A+\left(18\sqrt{B}\epsilon +2\sqrt{3(4A^3+27B\epsilon ^2)}\right)^{2/3}\right)^2}{\left(18\sqrt{B}\epsilon +2\sqrt{3(4A^3+27B\epsilon ^2)}\right)^{1/3}}\hfill \\ & & \times \frac{9\sqrt{B}\epsilon +\sqrt{3(4A^3+27B\epsilon ^2)}+6^{1/3}A\left(18\sqrt{B}\epsilon +2\sqrt{3(4A^3+27B\epsilon ^2)}\right)^{1/3}}{\left(9\sqrt{B}\epsilon +\sqrt{3(4A^3+27B\epsilon ^2)}\right)^2},\hfill \end{array}$$
(37)
where we have used the notations $`A(1{\displaystyle \frac{\delta ^2}{24}})c`$ and $`B{\displaystyle \frac{\delta ^2}{6\mathrm{}^2\omega ^2}}c^3`$.

The mean value of the total energy of the system has the form,

$$\overline{E}=_0^{\mathrm{}}\epsilon n(\epsilon )a(\epsilon )𝑑\epsilon .$$
(38)
To calculate $`\overline{E}`$, we expand the right hand side of (37) to the third order in the neighborhood of $`q=1`$ and rewrite the density of energy levels as the form

$$a(\epsilon )=\frac{8\pi V}{3h^3c^3}\left(3\frac{17}{8}\delta ^2\right)\epsilon ^2.$$
(39)
Then we get

$$\overline{E}=\frac{8\pi ^5V}{15}\frac{(kT)^4}{(hc)^3}\pi V\delta ^2\left(\frac{8\pi ^2}{3}+\frac{11\pi ^4}{15}+48\zeta (3)\right)\frac{(kT)^4}{(hc)^3},$$
(40)
where $`\zeta (n)`$ is the Riemann zeta function.

We can further calculate the grand partition function $`\mathrm{ln}\mathrm{\Xi }`$ through the thermodynamic relation $`\overline{E}=\frac{}{\beta }\mathrm{ln}\mathrm{\Xi }`$ ,

$$\mathrm{ln}\mathrm{\Xi }=\frac{8\pi ^5V}{45h^3c^3}(kT)^3\frac{\pi V\delta ^2}{h^3c^3}\left(\frac{8\pi ^2}{9}+\frac{11\pi ^4}{45}+16\zeta (3)\right)(kT)^3.$$
(41)
It is a straightforward calculation to obtain thermodynamic functions by making use of the grand partition function. The pressure of the deformed Bose gas is

$$P=\frac{1}{\beta }\frac{}{V}\mathrm{ln}\mathrm{\Xi }=\frac{8\pi ^5}{45h^3c^3}(kT)^4\frac{\pi \delta ^2}{h^3c^3}\left(\frac{8\pi ^2}{9}+\frac{11\pi ^4}{45}+16\zeta (3)\right)(kT)^4.$$
(42)
The free energy of the deformed Bose gas is of the form,

$$F=\frac{8\pi ^5V}{45}\frac{(kT)^4}{(hc)^3}+\pi V\delta ^2\left(\frac{8\pi ^2}{9}+\frac{11\pi ^4}{45}+16\zeta (3)\right)\frac{(kT)^4}{(hc)^3}.$$
(43)
The entropy of the deformed Bose gas is as

$$S=\frac{UF}{T}=\frac{32\pi ^5}{45}\frac{k^4T^3}{h^3c^3}V\pi \delta ^2\left(\frac{32\pi ^2}{9}+\frac{44\pi ^4}{45}+64\zeta (3)\right)\frac{k^4T^3}{h^3c^3}V.$$
(44)
The first term in the above equation is just the entropy of an ordinary Bose gas. And the terms following it denote the effect of noncommutative deformation of geometry on the entropy. It is clear that the effect can not be ignored at least in the high temperature region, where the early universe lived. We give a plot of the entropy vs temperature for different deformation parameters in Fig. 3.

FIG.3 The entropy per unit volume of a deformed Bose gas

( start from up $`q=1.0,1.1,1.2,1.3,1.4`$).

Finally, the specific heat of the deformed Bose gas is

$$C_V=T(\frac{S}{T})_V=\frac{32\pi ^5V}{15}\frac{k^4T^3}{h^3c^3}\pi V\delta ^2\left(\frac{32\pi ^2}{3}+\frac{44\pi ^4}{15}+192\zeta (3)\right)\frac{k^4T^3}{h^3c^3}.$$
(45)
At the case of $`q=1`$, we recover the result for an ordinary Bose gas. To show the effect of deformation on the specific heat, we give a plot of the specific heat per unit volume of a deformed Bose gas vs temperature for different deformation parameters in Fig. 4.

From the above discusses, we can also obtain the Planck’s formula for the deformed Bose gas,

FIG.4 The specific heat per unit volume of a deformed Bose gas

( start from up $`q=1.0,1.1,1.2,1.3,1.4`$).

$$\begin{array}{ccc}\hfill u(\omega ,T)& =& \frac{\mathrm{}\omega ^3}{c^3\pi ^2(e^{\beta \mathrm{}\omega }1)}\frac{17\delta ^2\mathrm{}\omega ^3}{24c^3\pi ^2(e^{\beta \mathrm{}\omega }1)}\frac{\delta ^2\beta \mathrm{}^2\omega ^4e^{\beta \mathrm{}\omega }}{2c^3\pi ^2(e^{\beta \mathrm{}\omega }1)^2}\hfill \\ & & \frac{\delta ^2\beta \mathrm{}^2\omega ^4e^{\beta \mathrm{}\omega }}{c^3\pi ^2(e^{\beta \mathrm{}\omega }1)^3}\frac{\delta ^2\beta \mathrm{}^2\omega ^4e^{\beta \mathrm{}\omega }}{2c^3\pi ^2(e^{\beta \mathrm{}\omega }1)^4}.\hfill \end{array}$$
(46)
FIG.5 The Planck’s formula of a deformed Bose gas

( start from up $`q=1.0,1.1,1.2,1.3,1.4`$, here we have used the notation $`f=\beta \mathrm{}\omega `$.)

In Fig.5, we give a plot of the Planck’s formula for the deformed Bose gas. The first term in equation(46) describes an ideal gas consisting of photons. The linearity of the equations of electrodynamics on commutative spacetime guarantees that photons do not interact with one another. However, the other terms in this equation denote interactions among photons or nonlinearity of equations of electrodynamics caused by noncommutative deformation of spacetime. At high frequencies ($`\beta \mathrm{}\omega 1`$), we have

$$\begin{array}{ccc}\hfill u(\omega ,T)& =& \frac{\mathrm{}\omega ^3e^{\beta \mathrm{}\omega }}{c^3\pi ^2}\frac{17\delta ^2\mathrm{}\omega ^3e^{\beta \mathrm{}\omega }}{24c^3\pi ^2}\frac{\delta ^2\beta \mathrm{}^2\omega ^4e^{\beta \mathrm{}\omega }}{2c^3\pi ^2}\hfill \\ & & \frac{\delta ^2\beta \mathrm{}^2\omega ^4e^{2\beta \mathrm{}\omega }}{c^3\pi ^2}\frac{\delta ^2\beta \mathrm{}^2\omega ^4e^{3\beta \mathrm{}\omega }}{2c^3\pi ^2}.\hfill \end{array}$$
(47)
This is the Wien’s formula of the deformed Bose gas. We give a plot of the Wien’s formula for different deformation parameters in Fig.6. In spite of a complex form of the Planck’s formula and the Wien’s formula, we note that the density of the spectral frequency distribution of the energy of the deformed Bose gas has a maximum at a frequency $`\omega _m`$ ($`\beta \mathrm{}\omega _m2.8`$), which is roughly independent of deformation parameter. Thus the displacement law is still remained for the deformed Bose gas.

FIG.6 The Wien’s formula of a deformed Bose gas

( start from up $`q=1.0,1.1,1.2,1.3,1.4`$, here we have used the notation $`f=\beta \mathrm{}\omega `$).

In the opposite limiting case of low frequencies ($`\beta \mathrm{}\omega 1`$), we get the Rayleigh-Jeans formula for the deformed Bose gas,

$$\begin{array}{ccc}\hfill u(\omega ,T)& =& \frac{kT\omega ^2}{c^3\pi ^2}\frac{17\delta ^2kT\omega ^2}{24c^3\pi ^2}\frac{\delta ^2kT\omega ^2}{2c^3\pi ^2}(\frac{\mathrm{}\omega }{kT}+1)\hfill \\ & & \frac{\delta ^2(kT)^2\omega }{c^3\pi ^2\mathrm{}}(\frac{\mathrm{}\omega }{kT}+1)\frac{\delta ^2(kT)^3}{2c^3\pi ^2\mathrm{}^2}\left(\frac{\mathrm{}\omega }{kT}+1\right).\hfill \end{array}$$
(48)
Fig.7 shows a graph corresponding to the Rayleigh-Jeans formula with different deformation parameters.

FIG.7 The Rayleigh-Jeans formula of a deformed Bose gas

( start from up $`q=1.00,1.06,1.12,1.18`$).

We obtain the Stefan-Boltzmann’s law for the deformed Bose gas of the form,

$$J=\frac{\pi ^2k^4}{60\mathrm{}^3c^2}T^4\frac{\delta ^2k^4}{32\pi ^2\mathrm{}^3c^2}\left(\frac{8\pi ^2}{3}+\frac{11\pi ^4}{15}+48\zeta (3)\right)T^4.$$
(49)
FIG.8 The Stefan-Boltzmann law of a deformed Bose gas

( start from up $`q=1.0,1.1,1.2,1.3,1.4`$).

Similar to the case of an ideal Bose gas, the total energy of the deformed Bose gas is proportional to the fourth power of the temperature. But the Stefan-Boltzmann’s constant is effectively reduced by the noncommutative deformation of geometry. The fact gives also a signal that the deformed Bose gas is a strongly correlated system. We give a plot of the Stefan-Boltzmann’s law for different deformation parameters in Fig.8.

In an adiabatic expansion (or compression) process, the entropy of the system is a constant and the valume and temperature are related by $`VT^3=`$constant. From (42), we can get the equation of state,

$$PV^{4/3}=c_q.$$
(50)
## 5 Concluding Remarks

The noncommutative geometry plays a very important role in the study of the Planck scale physics. Quantum field theory and superstring theory (or M theory) on noncommutative spacetime have been studied for a long time. It is well-known that the most realistic laboratory for testing the Planck scale physics is the early universe. The very early universe was opaque due to constant interchange of energy between matter and radiation. It is clear that the radiation field at the early universe can not be described well by an ideal Bose gas system. There should be strong correlations among photons at this stage of evolution of the universe. In fact, attention has been paid recently to the investigation on the deformed radiation dominated universe. Preliminary results show that spacetime quantization gives rise to inflation as the radiation temperature exceeds the Planck temperature. Of course, the study in this direction is only on its very beginning. Many efforts should be made to set up a self-contained inflation theory on noncommutative spacetime. We still have little knowledge on thermodynamics of the deformed radiation field.

From the linearity of the equations of electrodynamics on commutative spacetime, we know that photons do not interact with one another. The thermodynamics of the photonic field can be described well by an ideal Bose gas system. However, the photons in the deformed radiation field have strong correlations among themselves. The strong correlations among the photons are caused by the nonlinearity of equations of electrodynamics on noncommutative spacetime. It is a difficult task to deal with the interactions among the photons. Numerous papers have been published- on the thermodynamics of the deformed photonic field. However, we still have not got a satisfactory distribution function of the deformed Bose gas.

In this paper, we have treated the deformed photonic field as a strongly correlated system by making use of the Green’s function theory. Equations of motion of the Green’s function for the deformed photonic field was set up. The Tyablikov procedure was used to decouple these equations of motion. A general form of statistical distribution function for the deformed photonic field was got systematically. Based on this distribution function, we calculated the free energy, the entropy and the specific heat of the deformed photonic field in detail. The equation of state was obtained. We found that the Planck’s formula for the deformed Bose gas is really different from the ordinary one. We got terms in Planck’s formula, which correspond to the interactions among photons. Similar to the case of an ideal Bose gas, the total energy of the deformed Bose gas is proportional to the fourth power of the temperature. But the Stefan-Boltzmann’s constant is effectively reduced by the noncommutative deformation of geometry. We noted that the density of the spectral frequency distribution of the energy of the deformed Bose gas has a maximum at a frequency $`\omega _m`$ ($`\beta \mathrm{}\omega _m2.8`$), which is roughly independent of deformation parameter. Thus the displacement law is still remained for the deformed Bose gas. The results obtained in the paper can be used to investigate properties of the early universe (e.g. inflation model on noncommutative spacetime). All of these works are in progress.

Acknowledgements

The work was supported in part by the Natural Science Foundation of China.