text:
# Radial Temperature Profiles of 11 Clusters of Galaxies Observed With BeppoSAX

## 1. Introduction

Knowledge of the radial temperature profile of the hot gas contained within galaxy clusters is a crucial element in determining the total gravitational mass of clusters. Through the equation of hydrostatic equilibrium, the total mass of the cluster can be derived if the density gradient, temperature, and temperature gradient of the gas is known. The latter quantity is the most difficult quantity to obtain, and it has generally been assumed that the gas is isothermal in most previous approaches, since collimated X-ray instruments such as Ginga and EXOSAT could not determine the temperature structure of clusters.

The assumption of isothermality of the hot gas has been called into question in recent years, mainly as a result of studies done with ASCA. The ability of ASCA to perform spatially-resolved spectroscopy over the 1–10 keV energy range made it the first X-ray instrument capable of addressing the issue of temperature structure in hot clusters. Many ASCA studies have found that the gas within clusters is not isothermal, but decreases with increasing radius, in some cases up to a factor of two (e.g., Markevitch 1996; Markevitch et al. 1998; Markevitch et al. 1999). However, other studies of clusters using ASCA data have come to the conclusion that the gas is largely isothermal (e.g., White 1999; Fujita et al. 1996; Ohashi et al. 1997; Kikuchi et al. 1999), at least outside of the cooling radius of cooling flow clusters. A likely cause of this discrepancy is the handling of the large, energy-dependent point spread function (PSF) of ASCA that preferentially scatters hard X-rays. This creates an artificial increase in the temperature profile with radius if not dealt with properly. The PSF-correction method applied by Markevitch et al. consistently leads to significantly decreasing temperature profiles, while other methods (most notably the method of White 1999) lead to isothermal profiles.

The discrepancy among the different PSF-correction methods prompted an analysis of ROSAT PSPC data by Irwin, Bregman, & Evrard (1999). Although ROSAT was only sensitive to photon energies up to 2.4 keV and was therefore not the most ideal instrument with which to study hot clusters, large (factor of two) differences in temperature should have been detected, but were not. The composite X-ray “color” profiles for 26 clusters in the Irwin et al. (1999) survey indicated isothermality outside of the cooling radius. In fact, a 20% temperature drop within 35% of the virial radius was ruled out at the 99% confidence level.

In this paper, we attempt to resolve the temperature profile discrepancy using BeppoSAX data. BeppoSAX is sensitive to photon energies up to 10.5 keV, and has a half-power radius that is one-half that of the ASCA GIS instrument. In addition, the PSF of BeppoSAX is only weakly dependent on energy. Thus, BeppoSAX is better-suited for determining temperature profiles for clusters of galaxies than previous X-ray telescopes. Using a sample of 11 clusters found in the BeppoSAX archive, we derive radial temperature profiles for each cluster. In a future paper, we will discuss the abundance profiles of the 11 clusters. Throughout this paper, we assume $`H_0=50`$ km s-1 Mpc-1 and $`q_0=0.5`$.

## 2. Sample and Data Reduction

From the BeppoSAX Science Data Center (SDC) archive (available at http://www.sdc.asi.it/sax\_main.html) we have obtained data for seven cooling flow clusters (A85, A496, A1795, A2029, A2142, A2199, and 2A0335+096) and four non-cooling flow clusters (A2163, A2256, A2319, and A3266). All the clusters have redshifts in the range $`z=0.030.09`$ except A2163 ($`z=0.203`$). We analyze data from the Medium Energy Concentrator Spectrometer (MECS), which consists of three identical gas scintillation proportional counters (two detectors after 1997 May 9) sensitive in the 1.3–10.5 keV energy range. A detailed description of the MECS is given in Boella et al. (1997). The event files were subjected to the standard screening criteria of the BeppoSAX SDC.

Since our goal is spatially-resolved spectroscopy, accounting for scattering from the PSF of the MECS is important. As stated above, scattering from the PSF has an enormous impact on the temperature profiles derived from ASCA data. Fortunately, the detector + telescope PSF of BeppoSAX is nearly independent of energy, unlike the ASCA GIS. This is because the Gaussian PSF of the MECS detector improves with increasing energy, while the PSF of the grazing incidence Mirror Unit degrades with increasing energy (D’Acri, De Grandi, & Molendi 1998), leading to a partial cancellation when these two effects are combined. Still, it is important to account for the PSF accurately when deriving temperature profiles.

To correct for the PSF we have used the routine effarea, available as part of the SAXDAS 2.0 suite of BeppoSAX data reduction programs. This routine is described in detail in Molendi (1998) and D’Acri et al. (1998). Briefly, effarea creates an appropriate effective area file that corrects for vignetting and scattering effects for an azimuthally-symmetric circular or annular region. It does so by creating correction vectors that are a function of energy, and which when multiplied by the observed spectrum yields the corrected spectrum. The surface brightness profile of the cluster (determined from the analysis of ROSAT data by Mohr, Mathiesen, & Evrard 1999 and Ettori & Fabian 1999) is convolved with the PSF of the MECS in order to determine the extent to which scattering from other regions of the cluster have contaminated the emission from the extraction region in question. This information is incorporated into the auxiliary response file (the .arf file), which is subsequently used in the spectral fitting. The correction to the observed spectrum is modest; D’Acri et al. (1998) and Kaastra, Bleeker, & Mewe (1998) found only small changes between the uncorrected and corrected temperature profiles for Virgo and A2199, respectively. The mismatch in energy bandpasses between the ROSAT (0.2–2.4 keV) surface brightness profile and BeppoSAX (1.65–10.5 keV) does not appear to have significantly affected the results.

For each cluster, spectra were extracted from concentric annular regions centered on the peak of emission of the cluster with inner and outer radii of $`0^{}2^{}`$, $`2^{}4^{}`$, $`4^{}6^{}`$, and $`6^{}9^{}`$. We also extracted one global spectrum from $`0^{}9^{}`$. At $`9^{}`$ the telescope entrance window support structure (the strongback) becomes a factor. In addition, for off-axis angles greater than $`10^{}`$, the departure of the PSF from radial symmetry becomes noticeable (Boella et al. 1997). This coupled with the fact that some of the clusters have poor photon statistics outside of $`10^{}`$ prompted us to end our profiles at $`9^{}`$. At this radius, our temperature profiles extend out to 55% of the virial radius, $`r_{virial}`$, for A2163 and 17%–33% for the other clusters, where $`r_{virial}=3.9(T/10\mathrm{keV})^{1/2}`$ Mpc. Background was obtained from the deep blank sky data provided by the SDC. We used the same region filter to extract the background as we did the data, so that both background and data were affected by the detector response in the same manner. The energy channels were rebinned to contain at least 25 counts.

The procedure outlined above does not fit the spectrum of the various regions within the cluster simultaneously. Instead, it assumes a uniform spectrum throughout when correcting for contamination from other regions. Whereas this is not important for the innermost bin (since very few photons are scattered in from larger radii compared to the number of photons truly belonging in the innermost bin), this might affect the outer bins if the temperature profile is varying strongly. This does not seem to be the case though. The spectrum correction vectors presented in D’Acri et al. (1998) for A2199 were quite modest. Other than their innermost bin (which loses some flux via scattering but does not receive much scattered flux from exterior bins), the corrections amounted to 5% or less for energies above 3 keV. In addition, if the exterior bins were significantly contaminated by emission from the interior of the cluster that was at a considerably different temperature, it is likely that no single-component thermal model would give an adequate fit to the data. However, all the spectral fits of the third and fourth spatial bins of the clusters in our study were adequate, with nearly all fits having $`\chi _\nu <1.05`$. This coupled with the fact that the PSF-corrected temperatures were not significantly different from the uncorrected temperatures (see § 3.2) indicates that our results are not strongly affected by not fitting the spectra from different regions simultaneously.

For one of the lower temperature clusters (A2199) a long pointed ROSAT PSPC observation was available in the HEASARC archive for which no ROSAT-determined temperature profile had been published. The observation (RP800644N00) was filtered such that all time intervals with a Master Veto Rate above 170 counts s-1 were excluded, in order to discard periods of high background. This resulted in a net exposure of 34,232 seconds. Background was taken from an annulus with inner and outer radii of $`30^{}`$ and $`40^{}`$. The background was scaled to and subtracted from the source spectra, which were subsequently binned such that each energy channel contained at least 25 counts. Energy channels below 0.2 keV were ignored in the fit.

## 3. Temperature of the Hot Gas

### 3.1. Global Temperatures

Using XSPEC we used a MEKAL model with an absorption component fixed at the Galactic value in all spectral fits, except the ROSAT PSPC observation of A2199. The temperature, metallicity, and redshift were allowed to vary. Spectral fits of BeppoSAX MECS data of the Perseus cluster with the redshift fixed showed significant residuals in the iron line complex region around 6.7 keV (R. Dupke 1999, private communication). These residuals disappeared when the redshift was allowed to vary. The cause of this feature is a systematic shift of 45–50 eV in the MECS channel–to–energy conversion (F. Fiore 1999, private communication). This systematic shift was evident in our sample; when the redshift was allowed to vary, the measured redshift was less than the optically-determined redshift in all 11 clusters, and inconsistent with the optically-determined redshift at the 90% confidence level for eight of them. A modest decrease in the reduced $`\chi ^2`$ also occurred for most of the clusters when the redshift was allowed to vary. However, freeing the redshift did not affect the values obtained for the temperature and metallicity significantly (less than a 5% change in either quantity).

With the model described above, we fit the data from MECS2 and MECS3 (and MECS1 when available) separately, but with the same normalization. In accordance with the Cookbook for BeppoSAX NFI Spectral Analysis energy channels below 1.65 keV and above 10.5 keV were ignored in the fit. On the whole, the fits to the global spectra were rather poor, ranging from $`\chi _\nu ^2=1.131.8`$ for 170–556 degrees of freedom. Inspection of the residuals revealed that the poor fits resulted from excess emission in the 1.65–3.0 keV range. The best-fit spectrum from the non-cooling flow cluster A2256 is shown in Figure 1 and illustrates the positive residuals below 3.0 keV. This effect was even more pronounced in the cooling flow clusters. We performed the fits again, this time only using data in the 3.0–10.5 keV range. The fits were much better, with the fits to all clusters having $`\chi _\nu ^21.20`$. In addition, the global temperature values for the fits performed in the 3.0–10.5 keV range were much closer to the global values determined from ASCA data (see Table 1). Ten of the 11 temperatures derived from the 1.65–10.5 keV fit were below the ASCA value, and in eight cases the 90% error bars did not overlap. Conversely, nine of the 11 temperatures derived from the 3.0–10.5 keV fit have 90% errors bars that overlap with the error bars from the ASCA temperatures. All quoted errors are 90% confidence levels unless otherwise noted.

Since the improvement in $`\chi _\nu ^2`$ when channels in the 1.65–3.0 energy range were excluded was more pronounced for the cooling flow clusters than in the non-cooling flow clusters, we investigated the possibility that the excess emission below 3.0 keV was a result of a cooling flow component. Indeed, the addition of a cooling flow model to the MEKAL model provided a substantially better fit in the 1.65–10.5 keV case, and in some cases the fit became formally acceptable. However, the inferred cooling rates were several hundred solar masses per year higher than previous published values (e.g., Peres et al. 1998). In fact, cooling rates of several hundred solar masses per year were found for the non-cooling flow clusters A2163, A2256, A2319, and A3266. In addition, significant cooling rates were found for spectra extracted from regions of the cluster far from the cluster center, where no cooling gas should be found. This clear contradiction illustrates how a physically implausible model can still yield good spectral fits, and the danger in interpreting such a result. We conclude that although some of the excess emission in the 1.65–3.0 keV range is from cooling gas, there is a clear excess of soft emission beyond what is expected from cooling gas, possibly due to uncertainties in the calibration of the MECS instruments. Since including this energy range leads to global temperatures significantly below the ASCA value, we only fit the data in the 3.0–10.5 keV range for the remainder of the paper. Given the good agreement in the 3.0–10.5 keV fit and the ASCA-determined temperatures, we are confident of the calibration of BeppoSAX above 3.0 keV.

### 3.2. Radial Temperature Profiles

The PSF-corrected and -uncorrected radial temperature profiles for each of the 11 clusters are shown in the left panels of Figure 2. As was found in previous studies, correction for the MECS PSF does not seriously affect the temperature profile. In the right panels are our BeppoSAX temperature profiles along with temperature profiles derived from other ASCA and BeppoSAX studies for comparison. We note that for A85, A496, A1795, A2029, A2142, and A2199 the innermost bin has been fit with a cooling flow component in addition to a thermal model for the profiles of Markevitch et al. (1998), Markevitch et al. (1999), and Sarazin, Wise, & Markevitch (1998), which accounts for the large discrepancy in this bin compared to other studies that fit the spectra with only single-component models.

A85: We have included the subclump to the south of the cluster center in our analysis. Outside of the cooling radius, the temperature profile increases moderately from 6.4 keV to 7.9 keV at a significance level of $`2.2\sigma `$. This is in agreement with the ASCA analysis by White (1999), but contrasts with Markevitch et al. (1998) who found a profile that decreased from 8.0 keV in a $`1\stackrel{}{\mathrm{.}}56^{}`$ annular bin to 6.3 keV in a $`6^{}12^{}`$ annular bin with ASCA data. Pislar et al. (1997) and Kneer et al. (1995) analyzed ROSAT PSPC data for this cluster and found a roughly isothermal profile outside the cooling region. However, the ROSAT analysis found a significantly lower temperature, with most of the cluster below 5 keV. This is possibly due to a gain calibration problem of the ROSAT PSPC instrument.

A496: This nearby cluster ($`z=0.0326`$) has a cooling rate of $`\dot{M}=95M_{}`$ yr-1 (Peres et al. 1998). The temperature is lowest in the center (typical of a cooling flow), and is consistent with a constant value of 4.5–5.0 keV at larger radii. This is consistent with the ASCA analysis of Dupke & White (1999), who did not perform a PSF correction. For low temperature clusters, the PSF does not introduce a significant spurious positive gradient to the temperature profile (Takahashi et al. 1995). The profile of White (1999) also agreed for the most part with our profile, although one of their bins deviated by $`2\sigma `$ from ours. Markevitch et al. (1999) found a continuous drop in temperature of 5.6 keV to 3.5 keV from a $`2^{}5^{}`$ annular bin to a $`10^{}17^{}`$ annular bin.

A1795: We find a temperature profile consistent with a temperature of 6–7 keV outside of the cooling flow region. A significant decrease in temperature was not found in the center for this cooling flow cluster. This result is similar to the ASCA results of Mushotzky et al. (1995; uncorrected for the PSF), White (1999), and Ohashi et al. (1997), although the Ohashi et al. (1997) result found a somewhat lower overall temperature for this cluster. Markevitch et al. (1998) detected a decline in temperature with radius, but not at a high significance with ASCA. ROSAT found a low temperature in the inner cooling flow region, and a moderately increasing profile at larger radii (Briel & Henry 1996). The discrepancy in the innermost bin is likely the result of the low energy bandpass (0.1–2.4 keV) of ROSAT sampling lower temperature gas than the higher energy bandpass (3.0–10.5 keV) of BeppoSAX in the cooling region.

A2029: We find a basically flat temperature profile out to $`6^{}`$ and a marginally significant ($`1.7\sigma `$) rise from $`6^{}9^{}`$, consistent with the ASCA analysis of White (1999). This is the opposite trend found by Sarazin et al. (1998) with ASCA data, who found a decline in temperature from $``$9 keV to 6 keV outside of $`5^{}`$ (but with large errors in the outer regions). Still, they found that an isothermal profile was rejected at the $`>96\%`$ confidence level. This was one of the few clusters for which Irwin et al. (1999) found evidence for a statistically significant temperature decline with ROSAT data, although the drop did not occur until outside of 10. Molendi & De Grandi (1999) analyzed the same BeppoSAX data and found a temperature profile consistent with 8 keV out to 8 and dropping to 5 keV from $`8^{}12^{}`$ (see § 4.3 for a detailed comparison of the two analyses).

A2142: We find a very flat profile with a temperature of 8–9 keV. A similar result was found by White (1999), although the errors were large. Markevitch et al. (1998) found evidence for a temperature decline, but not at a high significance level. Henry & Briel (1996) analyzed ROSAT data and found temperatures of 10 keV or higher outside the cooling region, with a peak in the $`2\stackrel{}{\mathrm{.}}55^{}`$ bin.

A2163: This cluster has the highest redshift in our sample ($`z=0.203`$), and is also the hottest. The cluster does not possess a cooling flow and probably underwent a merger in the recent past (e.g., Elbaz, Arnaud, & Böhringer 1995). Since the centroid of the cluster lies over $`5^{}`$ from the detector center, we have excluded data from behind the strongback support structure and beyond. We have included data in the last annular bin out to $`12^{}`$ (as long as it fell inside the strongback) to improve the statistics in the last bin. Our profile extends out to 73% of the virial radius, considerably farther than any other cluster in our sample. We find that the cluster has a temperature of 10–11 keV out to 4, before experiencing a marginally significant ($`<2\sigma `$) rise in temperature at larger radii. The errors are quite large at large radii, although the temperature is greater than 8 keV at the 90% level. This agrees with the ASCA result of White (1999), but disagrees strongly with the result from ASCA and ROSAT by Markevitch et al. (1996), who found temperatures of 12.2$`{}_{1.2}{}^{}{}_{}{}^{+1.9}`$, 11.5$`{}_{2.9}{}^{}{}_{}{}^{+2.7}`$, and 3.8$`{}_{0.9}{}^{}{}_{}{}^{+1.1}`$ for annular bins of $`0^{}3^{}`$, $`3^{}6^{}`$, and $`6^{}10^{}`$ in extent, respectively. For these regions, we find temperatures of 10.1$`{}_{0.8}{}^{}{}_{}{}^{+0.9}`$, 11.7$`{}_{1.8}{}^{}{}_{}{}^{+2.6}`$, and 13.2$`{}_{4.6}{}^{}{}_{}{}^{+17.9}`$, respectively, using only data inside the strongback (90% errors). Thus, the outermost bin of the BeppoSAX data differs from that of Markevitch et al. (1998) at the 3.3$`\sigma `$ confidence level.

A2199: Outside of the cooling flow region, we find a constant temperature of about 4.5 keV. Analysis of the same BeppoSAX data by Kaastra et al. (1998) and of ASCA data by White (1999) found a similar result. However, Markevitch et al. (1999) found a steadily decreasing profile from 5.2 keV to under 5 keV at $`7\stackrel{}{\mathrm{.}}5`$, with a further decline at larger radii with the same ASCA data.

We have analyzed the ROSAT PSPC data for this cluster, and derived a temperature profile out to $`18^{}`$ (0.9 Mpc). The angular extent of this cluster is large, so to confirm that our background region was not significantly contaminated by the cluster emission, we derived a temperature

profile using a background annulus of $`40^{}50^{}`$ and found the same results as we did with a background annulus of $`30^{}40^{}`$. The temperature profile is shown in Figure 3. The profile shows a drop in the center indicative of a cooling flow. At large radii the profile is flat. However, the ROSAT-derived temperature is significantly lower than the value from ASCA or BeppoSAX. Whereas this might be expected in the central cooling flow region where ROSAT is sampling cooler gas than the other instruments because of its low energy bandpass, this tendency persists outside the cooling region. From $`2^{}9^{}`$, the temperature is $`3.2\pm 0.2`$ keV, whereas it is $`4.5\pm 0.1`$ for BeppoSAX. Data from ROSAT appears sometimes to have a tendency to measure lower temperatures than other instruments for clusters, such as A3558 (Markevitch & Viklinin 1997) and A85 (see above).

To compensate for this effect, we have adjusted the gain of the observation such that the temperature in the $`2^{}9^{}`$ region matched that of BeppoSAX. We excluded the inner $`2^{}`$ to avoid complications from the cooling flow region. An adjustment of 1.5% in the gain was necessary to bring the global temperature determined by the two instruments into agreement. A gain adjustment of this magnitude was within the range of values found by Henry & Briel (1996) when they analyzed five different pointings of A2142 with ROSAT. With this new gain value the temperature profile remains flat outside of the cooling region out to $`18^{}`$ (35% of the virial radius), albeit at a higher value than before. The main conclusion drawn from the ROSAT result of A2199 is that the temperature profile appears flat outside of the cooling flow region regardless of whether or not the gain was adjusted.

A2256: We find a flat temperature profile consistent with a temperature of 7 keV, in excellent agreement with the ASCA-determined profile of White (1999). This contrasts with the Markevitch (1996) result from ASCA which found the temperature to decrease from 8.7 keV to 7.3 keV from the $`0^{}6^{}`$ to $`6^{}11^{}`$, with a steep drop to 4 keV at larger radii. Markevitch (1996) claims that the ROSAT data confirm this result, albeit with larger errors, contrary to the claim of Briel, & Henry (1993) who analyzed the same ROSAT data and found a roughly isothermal profile.

A2319: We find a flat temperature profile out to $`6^{}`$ and a marginally significant ($`2\sigma `$) increase from $`6^{}9^{}`$. This is consistent with the White (1999) result. Markevitch (1996) also found an isothermal profile out to $`10^{}`$ with ASCA, and a decreasing profile at larger radii. The ROSAT data suggested isothermality (Irwin et al. 1999). Molendi et al. 1999 analyzed the same BeppoSAX data and found a flat temperature profile out to 16.

A3266: This probable merging cluster exhibits a slight increase in temperature in the center and levels off at radii out to 9. With ASCA data, Markevitch et al. (1998) found a steady decrease in temperature from almost 10 keV in the central $`2\stackrel{}{\mathrm{.}}5`$ to 5 keV outside of $`10^{}`$, and White (1999) also found a modestly decreasing profile (from 9 keV to 7.5 keV). The ROSAT data suggested isothermality although a modest decrease in temperature could not be ruled out (Irwin et al. 1999). De Grandi & Molendi (1999) analyzed the same BeppoSAX data and found a more substantial decrease in temperature, with a drop from 10 keV in the center to 4.5 keV out to 20.

2A0335+096: The coolest cluster in our sample, 2A0335+096 possesses a rather strong cooling flow ($`400M_{}`$ yr-1; Irwin & Sarazin 1995). The temperature is lowest in the innermost bin and levels off to a value of 3.4 keV out to 9. Kikuchi et al. (1999) analyzed the ASCA data for this cluster and found a flat temperature profile out to 10, while White (1999) found a flat profile out to 4 and a jump to 4.6 keV at larger radii.

In conclusion, the temperature profiles of the 11 clusters in our sample are roughly constant, and in general agreement with those derived from ASCA data by White (1999). This is in in disagreement with the profiles derived from the same ASCA data with the method of Markevitch et al. (1996), which find a significant decrease in temperature in many of the clusters. A more detailed comparison is given in § 4.

### 3.3. Normalized Temperature Profiles

In Figure 4a we plot the temperature profiles normalized to the global temperature for all 11 clusters versus radius in units of the virial radius, $`r_{virial}`$. The error bars represent the 1 $`\sigma `$ uncertainties. At small radii, the normalized profiles are typically less than one, owing to the presence of cooling flows in seven of the 11 clusters. As a result of the low temperature in the center, the outer regions are necessarily normalized to a value greater than one. To compensate for this effect, we have calculated global temperatures for the seven cooling flow clusters excluding the inner $`2^{}`$, and then normalized the $`2^{}4^{}`$, $`4^{}6^{}`$, and $`6^{}9^{}`$ bins by this temperature, while excluding the innermost bin. The result is shown in Figure 4b. Out to 20% of the virial radius, the temperature profiles appear flat. From 20% to 30% of the virial radius, the profiles rise somewhat, although the temperatures are not well constrained in this region. Most of the values are consistent with unity at the 1$`\sigma `$ confidence level.

A constant temperature model ($`T/T_{mean}=1`$) provided a good fit to the data ($`\chi ^2=43.7`$ for 37 degrees of freedom). A linear model of the form $`(T/T_{mean})=a+b(r/r_{virial})`$ provided a somewhat better fit ($`\chi ^2=36.0`$ for 35 degrees of freedom), with values of $`a=0.942`$ and $`b=0.440`$. The 90% confidence range on the slope $`b`$ was 0.123–0.752. This best-fit line is shown in Figure 4b, with the dashed lines representing the 90% confidence levels on the slope of the line. We find that a temperature drop of 14% from the center out to 30% of the virial radius can be ruled out at the 99% confidence level.

The data indicate that the gas is isothermal or mildly increasing in temperature out to 30% of the virial radius (1.0$`h_{50}^1`$ Mpc for a 7 keV cluster). Beyond this radius, the temperature profile may well decline. Most hydrodynamical cluster simulations predict a drop in temperature past 50% of the virial radius (see, e.g., Frenk et al. 1999). XMM will be the ideal instrument for determining the temperature profiles of clusters at large radii.

## 4. Comparison to Other Results

The presence of temperature gradients in clusters of galaxies has been a topic of hot debate in recent years, and there have been many investigations as to the temperature structure of clusters. Here, we summarize previous results grouped by X-ray telescope and compare them to our results.

### 4.1. Temperature Profiles Determined With ASCA

The claim of the existence of large temperature declines with radius was first suggested following the results from ASCA data using the PSF-correction technique of Markevitch et al. (1996), including Markevitch (1996), Markevitch et al. (1998), Sarazin et al. (1998), and Markevitch et al. (1999). Markevitch et al. (1998) observed 30 clusters with ASCA and found in general that the temperature profiles declined with radius (up to a factor of two), and that the decline was described well by a polytropic equation of state, i.e.,

$$T(r)\left(1+\frac{r^2}{a_x^2}\right)^{3\beta (\gamma 1)/2},$$
(1)
where $`r`$ is the projected distance, $`a_x`$ is the core radius, $`\beta `$ has its usual meaning in the context of isothermal $`\beta `$-models, and $`\gamma `$ is the polytropic index. The best-fit polytropic index was $`\gamma =1.24_{0.12}^{+0.20}`$ (90% confidence levels). This model clearly does not fit the BeppoSAX data. The polytropic fit led to $`\chi ^2=577`$ for 37 degrees of freedom.

It was pointed out by Irwin et al. (1999) that the trend in the temperature profile found by the method of Markevitch et al. (1996) differed from that found by other PSF-correction techniques from various authors (e.g., Ikebe 1995; Fujita et al. 1996; Kikuchi et al. 1999). Of the 28 clusters analyzed using the method of Markevitch et al. (1996) that did not show very strong evidence for a merger, 22 (14) of them were inconsistent with a constant temperature profile at the 70% (90%) probability level. Conversely, clusters analyzed using the method of Ikebe (1995), Fujita et al. (1996), or Kikuchi et al. (1999) were found to have basically flat temperature profiles, with all 11 clusters consistent with a constant temperature.

Of particular interest are clusters analyzed by more than one method, or clusters whose temperature is low enough that the PSF of ASCA does not influence the temperature profile significantly (for clusters below 5 keV). A399, A401, MKW3S, A1795, and A496 are examples where the method of Markevitch et al. (1996) leads to different trends in the radial temperature profiles than with other methods (Fujita et al. 1996; Kikuchi et al. 1996; Ohashi et al. 1997; Dupke & White 1999). The latter two are also in our sample, and show no evidence for a decline in temperature outside of the cooling region (Figure 2).

Recently, White (1999) has determined the temperature profiles of a large number of clusters observed with ASCA using the PSF-correction technique outlined in White & Buote (1999). He found that 90% of the 98 clusters in his sample were consistent with isothermality at the 3$`\sigma `$ confidence level. On a cluster by cluster comparison of the 11 clusters common in both samples, we find our temperature profiles are in excellent agreement with those of White (1999), with only two exceptions. Our temperature value of the outermost bin of A3266 is $``$2 keV higher than the value of White (1999). However, the 1$`\sigma `$ error bars just touch so this discrepancy is not at a high significance. Also, our temperature profile of A496 differs somewhat from his. This difference occurs around 4 and is significant at $``$$`2\sigma `$. However, given that we have 44 temperature values for the 11 clusters, it is not unreasonable that two of our values differ by $`2\sigma `$. In fact, this is the number of $`2\sigma `$ discrepancies one would expect from a statistical point of view.

Given the conflicting results regarding this topic, it is worthwhile to examine the PSF-correction techniques of the two largest studies, namely those using the techniques of Markevitch et al. (1996; hereafter method M) and White & Buote (1999; hereafter method WB). Both techniques assume a spatial distribution for the cluster emission; method M assumes the emissivity profile obtained from ROSAT PSPC data for the cluster in question, while method WB assumes the spatial distribution derived from a maximum-likelihood deconvolution of the ASCA image. Method M creates simulated events drawn from this spatial distribution and from an initial guess for the cluster spectrum, and convolves the events with the spatially-variable PSF. The input cluster spectrum is varied until the simulated data matches the actual data. Method WB ray-traces events drawn from the spatial distribution through the telescope optics and attempts to find the most likely association between events in the deconvolved and convolved planes on a PI energy bin by PI energy bin basis. Once an energy has been assigned to each event in the deconvolved plane, the event list can be analyzed via standard spectral-fitting procedures.

Both methods have their advantages and disadvantages. Whereas method WB employs only a spatially-invariant PSF, method M uses a PSF that varies with position. On the other hand, method M is reliant on ROSAT data to determine the emissivity profile, while method WB is not. It is not clear if the emissivity profile derived from the 0.2–2.0 keV ROSAT band is appropriate for use over the ASCA bandpass. The only way to determine which of these disadvantages are leading to incorrect results is to create simulated data with a known temperature and surface brightness profile and convolve it with the response of the instrument, and then see if the PSF-correction technique can recover the input temperature and surface brightness profiles. White & Buote (1999) have tested their code on simulated cluster data with a giant cooling flow model, a medium cooling flow model, an isothermal profile, and a profile that decreases by a factor of two from the center out to 20, as well as data of varying signal-to-noise ratios. In all cases, the original temperature and surface brightness profiles are recovered. Conversely, method M has not been as extensively tested on a variety of simulated temperature distributions and signal-to-noise ratios.

As a final note, it should be noted that Markevitch et al. (1998) used a cooling flow model to fit the innermost bin of cooling flow clusters, whereas White (1999) used a single component thermal model. This will naturally lead to a lower value for White (1999) for the innermost bin (typically $`5\%`$ of the virial radius). However, White (1999) claims that correcting for the cooling gas leads to an ambient core temperature within the innermost bin consistent with the outer regions of the cluster. Again, this is inconsistent with the results of Markevitch et al. (1998) who find a value for the ambient core temperature greater than the rest of the cluster.

### 4.2. Temperature Profiles Determined With ROSAT

Although the limited bandpass of the ROSAT PSPC (0.2-2.4 keV) precludes tight constraints to be put on the temperatures of hot clusters, large (factor of two) temperature changes should be detectable with large enough signal-to-noise ratios. Unfortunately, most clusters observed with ROSAT did not have good enough statistics to accomplish this on an individual basis. To circumvent this problem, Irwin et al. (1999) averaged together the radial color profiles (ratio of counts in various bands covering the ROSAT PSPC bandpass) of 26 clusters observed with the PSPC. If large-scale deviations from isothermality were common in clusters, such a feature lost in the noise for an individual cluster would become apparent when the clusters were added together. Although a drop in temperature was found in the center of cooling flow clusters (indicating that the method could indeed detect changes in temperatures even for hot clusters), the temperature profiles were flat outside of the cooling region out to 35% of the virial radius. It was found that a 20% temperature drop within 35% of the virial radius was ruled out at the 99% confidence level. This is in agreement with the BeppoSAX data presented here, where a decline in temperature of 14% out to 30% of the virial radius is ruled out at the 99% confidence level.

### 4.3. Temperature Profiles Determined With BeppoSAX

Several of the clusters in our sample have been analyzed by other authors: A2199 (Kaastra et al. 1998), A2029 (Molendi & De Grandi 1999), A2319 (Molendi et al. 1999), and A3266 (De Grandi & Molendi 1999). The profile of A2199 is fully consistent with ours. For the other three clusters, somewhat different steps were taken in the data reduction process than what we did. Whereas we excluded data below 3.0 keV, other authors have included data down to 2.0 keV. In addition, they have frozen the redshift whereas we have let it be a free parameter. We have re-analyzed the BeppoSAX data for these three clusters, this time including data down to 2.0 keV and freezing the redshift. Although inclusion of channels below 3.0 keV lowered the derived temperatures somewhat (see § 3.1), it did not change the overall shape of the profiles.

Our new profile of A2319 was consistent with that of Molendi et al. (1999). The profiles of A2029 and A3266 differed somewhat from Molendi & De Grandi (1999) and De Grandi & Molendi (1999). However, the differences occurred only in the 2$`{}_{}{}^{}4^{}`$ bins, with our 1$`\sigma `$ error bars overlapping their 1$`\sigma `$ error bars for the other three spatial bins. For both A2029 and A3266 our 2$`{}_{}{}^{}4^{}`$ temperature was $``$1.5 keV lower than those obtained by Molendi & De Grandi (1999) and De Grandi & Molendi (1999), and the difference was significant at the 3.3$`\sigma `$ and 1.9$`\sigma `$ confidence levels for A2029 and A3266, respectively. The cause of this discrepancy may be the choice of the assumed surface brightness profile, especially for the case of A2029, which possesses a large cooling flow. As mentioned in § 2, we have used the double-beta model profile obtained by Mohr et al. (1999) from ROSAT PSPC data. It is not stated where the other authors obtained their assumed surface brightness profiles. However, the agreement between the other three spatial bins (especially the first and fourth spatial bins) is encouraging, indicating a flat temperature profile out to $`9^{}`$ for A2029 and A2319, and a modestly decreasing profile for A3266. It should be noted that Molendi & De Grandi (1999) and De Grandi & Molendi (1999) find significant temperature drops in A2029 and A3266, respectively, at very large radii, where we have truncated our profiles because of the presence of the strongback, and the increasing asymmetry of the PSF at larger radii.

## 5. A Summary of Cluster Temperature Profiles

Given the effort put forth in determining temperature profiles for a large sample of clusters in recent years it is worthwhile to summarize the major contributions to this subject: Markevitch et al. (1998; and references within) and White (1999) with ASCA data, Irwin et al. (1999) with ROSAT PSPC data, and this current study with BeppoSAX data. The combined results are shown in Figure 5. The 98 temperature profiles of White (1999) (shown as crosses) have been normalized to the global temperature for each individual cluster and scaled in units of the virial radius. The long dashed lines enclose the middle 68% of the data points. The error bars are typically rather large for each data point, and have been excluded for clarity. The temperature profile derived by Markevitch et al. (1998) is shown as a dotted line, and represents a fit with polytropic index of $`\gamma =1.24`$. The results presented in this paper is shown as a solid line. The ROSAT PSPC result is not shown but is very similar to best-fit line derived here, only somewhat less constrained.

The results of White (1999), Irwin et al. (1999), and this study all point to the same general conclusion: outside of the cooling region but inside 30% of the virial radius there appears to be no decline in the temperature profile. White (1999) extends this result out to 45% of the virial radius, although it should be noted that large drops in the temperature have been found in A2029 and A3266 with BeppoSAX data (Molendi & De Grandi 1999; De Grandi & Molendi 1999) in this region. Outside of 45% of the virial radius White (1999) finds evidence for a decline in temperature although the statistics are sparse in this region. This drop is not surprising though since nearly all cluster simulations show a decline in temperature at large radii.

JAI thanks R. Dupke for many useful comments and conversations. We thank D. White for kindly providing us with his ASCA temperature profiles, and also the anonymous referee for many insightful comments and suggestions to improve the paper. This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center, and also the BeppoSAX Science Data Center. This work has been supported by Chandra Fellowship grant PF9-10009, awarded through the Chandra Science Center. The Chandra Science Center is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-39073.


format: text

text:
# 1 Introduction

## 1 Introduction

The inclusion of heavy quark effects in deeply inelastic scattering (DIS) is an interesting theoretical problem involving two hard scales in a perturbative analysis. This issue is also of phenomenological importance. The charm contribution to the total structure function $`F_2`$ at small $`x`$, at HERA, is sizeable, up to $`25\%`$. Through the charm contribution to scaling violations, the treatment of charm also has a significant impact on the interpretation of fixed target deeply inelastic scattering data. Thus, a proper description of charm contributions to deeply inelastic scattering is required for a global analysis of structure function data and a precise extraction of the parton densities in the proton.

At scales $`Q\mathrm{\Gamma }<M_H`$, the contribution to deeply inelastic scattering of a heavy quark of mass $`M_H`$ can be calculated in the so-called fixed-flavor-number (FFN) prescription from hard processes initiated by light quarks ($`u,d,s,\mathrm{}`$) and gluons, where all effects of the heavy quark ($`H`$) are contained in the perturbative coefficient functions. This prescription incorporates the correct threshold behavior, but for large scales, $`QM_H`$, the coefficient functions at higher orders in $`\alpha _s`$ contain potentially large logarithms $`\mathrm{ln}^i(Q^2/M_H^2)`$, which may need to be summed . Such a summation can be achieved by including the heavy quark as an active parton in the proton. The simplest approach incorporating this idea is the so-called zero-mass variable-flavor-number (ZM-VFN) prescription, where heavy quarks are omitted entirely below some scale $`Q_0M_H`$ and included as massless partons above this threshold. This prescription has been used in global analyses of parton distributions for many years, but it has an error of order $`M_H^2/Q^2`$ and is not suited for quantitative analyses unless $`QM_H`$.

Considerable effort has been devoted to including heavy quark effects in deeply inelastic scattering in such a way that the calculated structure functions match those of the FFN prescription in the region $`QM_H`$ while they match those of the ZM-VFN prescription for $`QM_H`$. Two prescriptions of this sort, the Aivazis–Collins–Olness–Tung (ACOT)  and the Thorne–Roberts  prescriptions have been used in recent global analyses of parton distributions . More recently, additional variable-flavor-number prescriptions with non-zero mass have been defined in the literature . If one could sum perturbation theory, the calculated structure functions should be identical for any prescription that does not neglect the heavy quark mass. However, the way of ordering the perturbative expansion is not unique, so that the results generally differ at any finite order in perturbation theory.

In this paper, we will investigate a modification of the ACOT prescription advocated by Collins . It has the advantage of being simple to state and of allowing relatively simple calculations. This simplicity should be convenient for phenomenological analyses at the Born level. In addition, it could be crucial for implementing a variable-flavor-number prescription with non-zero mass at next-to-leading order in global analyses of parton distributions.

## 2 Schemes and partons

We consider the structure function $`F_2(x,Q)`$. The observable structure function can be written in terms of parton distribution functions $`f`$ and a calculable partonic structure function $`\widehat{F}_2`$ as

$$F_2(x,Q)=_x^1𝑑\xi \underset{a}{}f_{a/p}(\xi ,\mu )\widehat{F}_2(a,x/\xi ,\mu /Q,\alpha _s(\mu ))+𝒪(\mathrm{\Lambda }^2/Q^2).$$
(1)
This factorization formula has corrections of order $`\mathrm{\Lambda }^2/Q^2`$, where $`\mathrm{\Lambda }`$ is a typical scale of hadronic physics, perhaps the mass of the rho meson. Here $`\mu `$ is the renormalization/factorization scale. (For simplicity, we do not distinguish these two scales.) The partonic $`\widehat{F}_2`$ has a perturbative expansion in powers of $`\alpha _s(\mu )`$. There is a sum over parton types $`a=g,u,\overline{u},d,\overline{d},\mathrm{}`$, and there is a function $`f_{a/p}`$ for each parton.

To begin, we must pick a scheme for the definition of $`\alpha _s`$ and for the definition of the parton distribution functions. In fact, we will use multiple schemes, following the prescription of . Each scheme is designated by the number $`N`$ of “active” quark flavors. Scheme $`N`$ is designed to be most useful for physical scales $`Q`$ in the range $`M_N\mathrm{\Gamma }<Q\mathrm{\Gamma }<M_{N+1}`$, where $`M_N`$ is the mass of the $`N`$th quark flavor.

We define what we mean by the running coupling $`\alpha _s^N(\mu )`$ in scheme $`N`$ by defining how we perform renormalization. We renormalize using the CWZ prescription . Briefly, divergences involving active parton loops are removed with an   MS  subtraction, but when active parton external lines couple to a loop containing “non-partonic” lines (quarks $`N+1,N+2,\mathrm{}`$) the renormalization is by subtraction at zero external momentum. The running of $`\alpha _s^N(\mu )`$ is controlled by the usual   MS  renormalization group equation with the contributions from quarks $`1,\mathrm{},N`$ in the beta function.

In the scheme $`N`$ there are parton distributions for gluons and for quark flavors $`1,\mathrm{},N`$, but not for quark flavors $`N+1,\mathrm{}`$. The parton distribution functions are defined to be proton matrix elements of certain operators . The operator products are ultraviolet divergent and are renormalized according to the CWZ prescription. They obey a renormalization group equation – the usual DGLAP equation – in which contributions from quark flavors $`N+1,\mathrm{}`$ do not appear in the kernel.

Note that the parton distribution functions are non-perturbative objects. There is no question of neglecting any masses in the definition. The quark masses do not appear in the kernel of the evolution equation, but this is because renormalization counter terms are mass independent, not because the parton distributions themselves are mass independent.

Note also that by using CWZ renormalization throughout, quark loops for the non-partonic flavors decouple from calculations when the momentum scale is small compared to $`M_{N+1}`$. Thus for $`QM_{N+1}`$ it is a good approximation to leave the non-partonic flavors out of calculations altogether. (But leaving heavy flavors out of calculations is not part of the definition of scheme $`N`$ as used here; it is a separate approximation.) When one does leave the non-partonic flavors out of calculations, the remaining renormalizations are via the   MS  prescription. Thus one commonly refers to parton distributions in the prescription described here as   MS  parton distributions.

There is a perturbative connection between schemes $`N`$ and $`N+1`$ which we can represent by:

$$f_{a/p}^{N+1}(x,\mu )=f_{a/p}^N(x,\mu )+\underset{b}{}_x^1\frac{d\xi }{\xi }A_{ab}(x/\xi ,\mu /M_{N+1},\alpha _s(\mu ))f_{b/p}^N(\xi ,\mu ).$$
(2)
Here the index $`a`$ runs over $`g`$, the $`N`$ light quarks and their anti-quarks, and the heavy quark $`H`$ and anti-quark $`\overline{H}`$. For $`a=H,\overline{H}`$ we define $`f_{a/p}^N(x,\mu )=0`$ in the first term on the right hand side. The index $`b`$ runs over only the gluon and light quarks. At order $`\alpha _s`$, the only nonvanishing $`A_{Hb}`$ is $`A_{Hg}`$. Thus we may say that the $`H`$ distribution arises perturbatively from $`gH\overline{H}`$.

The $`\alpha _s^1`$ terms in the perturbative expansion of the $`A_{ab}`$ vanish at $`\mu =M_{N+1}`$. Thus one derives the simple matching conditon : $`f_{a/p}^{N+1}(x,M_{N+1})=f_{a/p}^N(x,M_{N+1})+𝒪(\alpha _s^2)`$. The analogous connection is known at order $`\alpha _s^2`$ but we do not repeat it here.111At order $`\alpha _s^2`$, the matching at $`\mu =M_{N+1}`$ is no longer continous in general .

## 3 Why is there an ambiguity?

When we construct the hard scattering cross section in the presence of heavy quarks, there is a factorization ambiguity that is not present in the case of light quarks. Let us see why. Consider two enormously simplified examples that illustrate the principle.

First, suppose that there are only gluons and one light quark $`L`$ with mass $`M_L\mathrm{\Gamma }<\mathrm{\Lambda }`$. Suppose, in addition, that the light quark is its own anti-particle, so that a quark $`L`$ and the anti-quark $`\overline{L}`$ are the same particle. Both the gluon and the light quark are considered to be active partons. In order to simplify the notation, let us take a moment $`𝑑xx^nF_2(x,Q)`$ of $`F_2`$, so that we get a factorization formula involving the corresponding moments of $`\widehat{F}_2`$ and of the parton distributions. The dependence on the moment number will not be indicated. To further simplify the notation, let us set the factorization and renormalization scale $`\mu `$ to $`Q`$. Then

$$F_2(Q)\widehat{F}_2(L,\alpha _s(Q))f_{L/p}(Q)+\widehat{F}_2(g,\alpha _s(Q))f_{g/p}(Q).$$
(3)
Now $`F_2`$ is an observable, so its definition is fixed. We have defined the parton distribution functions, and the two parton distribution functions are independent. Thus this equation all but fixes the definition of $`\widehat{F}_2(L,\alpha _s(Q))`$ and $`\widehat{F}_2(g,\alpha _s(Q))`$. The only possible modification would be to add terms proportional to powers of $`M_L^2/Q^2`$ to $`\widehat{F}_2(a,\alpha _s(Q))`$, at the cost of subtracting the same terms from the power suppressed remainder, $`𝒪(\mathrm{\Lambda }^2/Q^2)`$, in the factorization formula (Eq. (1)). The simplest solution, which is uniformly adopted, is not to allow a $`M_L^2/Q^2`$ dependence in $`\widehat{F}_2(a,\alpha _s(Q))`$.

Now suppose that there are only gluons and one (self-conjugate) heavy quark $`H`$ with mass $`M_H\mathrm{\Lambda }`$. Then in the scheme in which both the gluon and the heavy quark are considered to be active partons we have

$$F_2(Q)\widehat{F}_2(H,M_H/Q,\alpha _s(Q))f_{H/p}(Q)+\widehat{F}_2(g,M_H/Q,\alpha _s(Q))f_{g/p}(Q).$$
(4)
Here the $`\widehat{F}_2(a,M_H/Q,\alpha _s(Q))`$ (a=g,H) depend on the heavy quark mass and we cannot move terms of order $`M_H^2/Q^2`$ into the power suppressed corrections because only terms of order $`\mathrm{\Lambda }^2/Q^2`$ are allowed there. Thus it seems that we have no freedom. But we do: $`f_{H/p}(Q)`$ and $`f_{g/p}(Q)`$ are not independent. Since, according to Eq. (2), heavy quarks evolve from gluons, we have a relation of the form

$$f_{H/p}(Q)=V_{H/g}(\mathrm{ln}(Q/M_H),\alpha _s(Q))f_{g/p}(Q).$$
(5)
Here $`V`$ has a perturbative expansion that is obtained by solving the evolution equation. The first term has the form $`V\alpha _s\gamma \mathrm{ln}(Q/M)`$ where $`\gamma `$ is a constant. Using Eq. (5) in Eq. (4) we obtain

$$F_2(Q)\left\{\widehat{F}_2(H,M_H/Q,\alpha _s(Q))V_{H/g}(\mathrm{ln}(Q/M_H),\alpha _s(Q))+\widehat{F}_2(g,M_H/Q,\alpha _s(Q))\right\}f_{g/p}(Q).$$
(6)
Evidently, there is some freedom to move pieces from the first term in braces to the second. There is a constraint. For $`M_H/Q0`$, it is possible to neglect $`M_H`$ in the calculation of the $`\widehat{F}_2`$. On the other hand, $`V(\mathrm{ln}(Q/M_H),\alpha _s(Q))`$ does not have a smooth $`M_H/Q0`$ limit. This is not a problem in applications because by solving the evolution equation one sums the leading logarithms $`[\alpha _s\mathrm{ln}(Q/M_H)]^n`$ in $`V`$. It is important not to undo this summation. Thus the factorization scheme that we adopt should have the property that the functions $`\widehat{F}_2(a,M_H/Q,\alpha _s(Q))`$ have a finite limit as $`M_H/Q0`$. This still leaves us the option of adding a term like $`c\times (M_H^2/Q^2)`$ to $`\widehat{F}_2(H,M_H/Q,\alpha _s(Q))`$ and subtracting $`c\times (M_H^2/Q^2)V`$ from $`\widehat{F}_2(g,M_H/Q,\alpha _s(Q))`$.

Suppose now that we have calculated $`\widehat{F}_2(g,M_H/Q,\alpha _s(Q))`$ and $`\widehat{F}_2(H,M_H/Q,\alpha _s(Q))`$ in some convenient prescription – for example the prescription analyzed in  based on “on-shell” heavy quarks. Then we could define a new prescription with

$$\widehat{F}_2^{\mathrm{new}}(H,M_H/Q,\alpha _s(Q))=\widehat{F}_2^{\mathrm{old}}(H,0,\alpha _s(Q))$$
(7)
and

$`\widehat{F}_2^{\mathrm{new}}(g,M_H/Q,\alpha _s(Q))=\widehat{F}_2^{\mathrm{old}}(g,M_H/Q,\alpha _s(Q))`$
$`+\left\{\widehat{F}_2^{\mathrm{old}}(H,M_H/Q,\alpha _s(Q))\widehat{F}_2^{\mathrm{old}}(H,0,\alpha _s(Q))\right\}\times V_{H/g}(\mathrm{ln}(Q/M_H),\alpha _s(Q)).`$ (8)
In the following subsection, we shall give a prescription  based on this observation.

## 4 A prescription for resolving the ambiguity

Having seen the main idea, let us put the parton indices and the momentum fraction variables back. Suppose that there are $`N+1`$ quark flavors that we consider to be active. Let the heaviest active quark, quark $`N+1`$, be labelled $`H`$. Suppose that quark $`H`$ has a mass that is large compared to the hadronic mass scale $`\mathrm{\Lambda }`$. In the $`N+1`$ flavor scheme the factorization equation is

$$F_2(x,Q)_x^1𝑑\xi \underset{a}{}f_{a/p}(\xi ,\mu )\widehat{F}_2(a,x/\xi ,M_H/Q,\mu /Q,\alpha _s(\mu )),$$
(9)
where the sum over $`a`$ includes $`a=H`$ and $`a=\overline{H}`$.

The parton distributions for the gluon and $`N+1`$ quark flavors are obtained from the distributions in the scheme with only $`N`$ quark flavors. First we use the perturbative matching relation (Eq. (2)) at a scale near $`\mu =M_H`$, then we use the evolution equations to give the distribution functions at scale $`\mu `$. Since there are more output functions than input, we obtain a perturbative relation giving the heavy quark distribution functions at scale $`\mu `$ in terms of the light quark and gluon distribution functions at the same scale. This relation has the form

$$f_{H/p}(\xi ,\mu )=_\xi ^1\frac{d\tau }{\tau }\underset{aH\overline{H}}{}f_{a/p}(\tau ,\mu )V_{H/a}(\xi /\tau ;\mu /M_H,\alpha _s(\mu )),$$
(10)
with an analogous equation for the heavy anti-quark $`\overline{H}`$. Inserting this relation into Eq. (9), we have

$$F_2(x,Q)_x^1𝑑\xi \underset{aH,\overline{H}}{}f_{a/p}(\xi ,\mu )T_a(x/\xi ,M_H/Q,\mu /Q,\alpha _s(\mu )),$$
(11)
where

$`T_a(z,M_H/Q,\mu /Q,\alpha _s(\mu ))=\widehat{F}_2(a,z,M_H/Q,\mu /Q,\alpha _s(\mu ))`$
$`+{\displaystyle _z^1}𝑑\lambda \widehat{F}_2(H,z/\lambda ,M_H/Q,\mu /Q,\alpha _s(\mu ))V_{H/a}(\lambda ,\mu /M_H,\alpha _s(\mu ))`$
$`+{\displaystyle _z^1}𝑑\lambda \widehat{F}_2(\overline{H},z/\lambda ,M_H/Q,\mu /Q,\alpha _s(\mu ))V_{\overline{H}/a}(\lambda ,\mu /M_H,\alpha _s(\mu )).`$ (12)
We see that we have the same situation as in the simple example given earlier. Taking $`\mu `$ to be of order $`Q`$, one can shift contributions of order $`M_H^2/Q^2`$ between the hard scattering functions $`\widehat{F}_2`$ for $`H`$ and $`\overline{H}`$ and the corresponding functions for the light quarks and the gluon while keeping the functions $`T`$ unchanged for each light quark and gluon flavor $`a`$ and without ruining the property that all of the functions $`\widehat{F}_2`$ have finite $`M_H/Q0`$ limits.

This freedom can be exploited to make the calculation of the functions $`\widehat{F}_2`$ simpler. In particular, we can adopt a prescription proposed by Collins :

> Simplified ACOT (S-ACOT) prescription. Set $`M_H`$ to zero in the calculation of the hard scattering functions $`\widehat{F}_2`$ for incoming heavy quarks.

This observation tremendously simplifies the calculation of $`\widehat{F}_2`$ for $`a=H`$ as it reduces to that of the light-quark result. For example, at order $`\alpha _s^1`$, the $`\widehat{F}_2`$ functions for heavy quarks and light quarks are independent of $`M_H`$, and the calculation reduces to that of the simple massless result. The $`\alpha _s^1`$ gluon contribution to $`\widehat{F}_2`$ acquires an $`M_H`$ dependence when the gluon couples to a heavy quark loop, which is probed by the virtual photon ($`g\gamma H\overline{H}`$).

Note that the hard scattering functions $`\widehat{F}_2`$ obey a renormalization group equation that is different from the standard renormalization group equation that applies in the case of light flavors only. To see why this is so, imagine that the $`\widehat{F}_2`$ functions obeyed the usual renormalization equation and that $`\widehat{F}_2(H,z/\lambda ,M_H/Q,\mu /Q,\alpha _s(\mu ))`$ were independent of $`M_H`$ at some fixed value of $`\mu `$. Then $`\widehat{F}_2`$ for $`a=H,\overline{H}`$ would depend on $`M_H`$ at other values of $`\mu `$ because the standard renormalization group equation mixes the heavy and light flavors and would mix $`M_H`$ dependence into $`\widehat{F}_2(H,z/\lambda ,M_H/Q,\mu /Q,\alpha _s(\mu ))`$.

This Simplified ACOT prescription has the advantage of being simple to state. In addition, its calculational simplicity could be crucial for applying variable-flavor-number prescriptions with mass in analyses that go beyond first order in $`\alpha _s`$.

## 5 The ACOT and S-ACOT prescriptions at first order

In this section, we analyze the ACOT prescription and its simplified version, the S-ACOT prescription, at order $`\alpha _s^1`$. We consider neutral current deeply inelastic scattering from a proton target in a $`Q^2`$ regime in which quarks $`1,\mathrm{},N`$ can be considered as light while a single quark, $`H`$, is considered to be heavy or light depending on the value of $`Q^2`$. In order to keep the analysis simple, contributions from all heavier quarks are ignored. We further simplify the problem by supposing that the vector boson current that probes the proton couples only to the heavy quark $`H`$ and its anti-quark $`\overline{H}`$, but not to the lighter quarks or gluons. (In a crude approximation, this is like taking ordinary deeply inelastic scattering but demanding that $`H`$ appear in the final state.) Let $`F`$ be one of the structure functions $`F_1`$, $`F_2/x`$ or $`F_3`$ for our special vector boson. We choose the factorization and renormalization scales $`\mu `$ equal to $`Q`$. We write the factorization formula for $`F`$ in the shorthand notation

$$F=\underset{a}{}\widehat{F}_af_{a/p}+𝒪(\mathrm{\Lambda }^2/Q^2),$$
(13)
where the sum runs over all partons $`a`$ including $`H`$ and $`\overline{H}`$. The $``$ denotes a convolution, so that Eq. (13) means

$$F(x,Q)=\underset{a}{}\frac{d\xi }{\xi }\widehat{F}_a(x/\xi ,Q;\alpha _s(Q))f_{a/p}(\xi ,\mu )+𝒪(\mathrm{\Lambda }^2/Q^2).$$
(14)
The functions $`\widehat{F}_a`$ in Eq. (13) are the hard scattering functions. They have an expansion in powers of $`\alpha _s`$. Keeping the first two terms in this expansion gives

$`F`$ $`=`$ $`\widehat{F}_H^{(0)}f_{H/p}+\widehat{F}_{\overline{H}}^{(0)}f_{\overline{H}/p}`$ (15)
$`+\widehat{F}_g^{(1)}f_{g/p}+\widehat{F}_H^{(1)}f_{H/p}+\widehat{F}_{\overline{H}}^{(1)}f_{\overline{H}/p}`$
$`+𝒪(\alpha _s^2)+𝒪(\mathrm{\Lambda }^2/Q^2),`$
where $`\widehat{F}_a^{(n)}`$ is the order $`\alpha _s^n`$ contribution to $`\widehat{F}_a`$.

In the ACOT prescription, we choose the order zero hard scattering function for a heavy quark to be

$$\widehat{F}_{H,\mathrm{ACOT}}^{(0)}=F_H^{(0)}(M_H),$$
(16)
where $`F_H^{(0)}(M_H)`$ is the calculated structure function for scattering from an initial state heavy quark of mass $`M_H`$ that is on its mass-shell. (To be precise, the heavy quark transverse momentum is taken to be zero, and the structure functions $`F_1,F_2/x,F_3`$ are extracted from the tensor $`W^{\mu \nu }`$ using the usual formula but with the quark momentum $`k^\mu `$ replaced by a light-like vector $`\stackrel{~}{k}^\mu =k^\mu [M_H^2/(2uk)]u^\mu `$, where $`u^\mu `$ is a light-like reference vector in the plane of $`q^\mu `$ and $`k^\mu `$.) The function $`F_H^{(0)}(M_H)`$ is rather complicated (see Ref. ), so we do not reproduce it here. The definition of $`\widehat{F}_{\overline{H},\mathrm{ACOT}}^{(0)}`$ for a heavy anti-quark is analogous.

The order $`\alpha _s`$ gluon hard scattering function in Eq. (15) has three pieces:

$$\widehat{F}_{g,\mathrm{ACOT}}^{(1)}=F_g^{(1)}(M_H)F_H^{(0)}(M_H)f_{H/g}^{(1)}(M_H)F_{\overline{H}}^{(0)}(M_H)f_{\overline{H}/g}^{(1)}(M_H).$$
(17)
Here $`F_g^{(1)}(M_H)`$ is the calculated structure function for scattering from an initial state massless gluon that is on its mass-shell, using graphs with a heavy quark loop and a heavy anti-quark loop. The function $`f_{H/g}^{(1)}(M_H)`$ is the calculated $`\alpha _s^1`$ contribution to the distribution of heavy quarks in an on-mass-shell gluon. Similarly, $`f_{\overline{H}/g}^{(1)}(M_H)`$ is the calculated $`\alpha _s^1`$ contribution to the distribution of heavy anti-quarks in an on-mass-shell gluon. Both functions are given by

$$f_{H/g}^{(1)}(\xi ,Q)=f_{\overline{H}/g}^{(1)}(\xi ,Q)=\frac{\alpha _s(Q)}{2\pi }\mathrm{ln}\frac{Q^2}{M_H^2}P_{qg}^{(1)}(\xi ),$$
(18)
where $`P_{q/g}`$ is the usual gluon $``$ quark splitting function $`P_{qg}(\xi )=T_F(\xi ^2+(1\xi )^2)`$.

The first order hard scattering function for a heavy quark has a structure similar to that of the corresponding function for a gluon,

$$\widehat{F}_{H,\mathrm{ACOT}}^{(1)}=F_H^{(1)}(M_H)F_H^{(0)}(M_H)f_{H/H}^{(1)}(M_H).$$
(19)
Here $`F_H^{(1)}(M_H)`$ is the order $`\alpha _s^1`$ contribution to the structure function for scattering from an initial state heavy quark that is on its mass-shell, as given in . The function $`f_{H/H}^{(1)}(M_H)`$ is the calculated $`\alpha _s^1`$ contribution to the distribution of heavy quarks in an on-mass-shell heavy quark,

$$f_{H/H}^{(1)}(\xi ,Q)=C_F\frac{\alpha _s(Q)}{2\pi }\left[\frac{1+\xi ^2}{1\xi }\left\{\mathrm{ln}\left(\frac{Q^2}{(1\xi )^2M_H^2}\right)1\right\}\right]_+,$$
(20)
where the + subscript denotes the usual prescription,

$$_0^1𝑑\xi f(\xi )\left[F(\xi )\right]_+=_0^1𝑑\xi \{f(\xi )f(1)\}F(\xi ).$$
(21)
This result can be calculated easily from the definition  of   MS  parton distribution functions or it can be extracted from the ACOT subtraction terms in .

This defines the ACOT prescription at order $`\alpha _s^1`$. The prescription has two important properties.

* Property 1. For $`Q/M\mathrm{}`$, the hard scattering functions in Eq. (15) approach the hard scattering functions of the ZM-VFN prescription, in which the heavy quark $`H`$ is taken as a parton with zero mass.

To be specific, we have the relations:

$`\widehat{F}_{H,\text{ZM-VFN}}^{(0)}`$ $`=`$ $`F_H^{(0)}(0),`$
$`\widehat{F}_{g,\text{ZM-VFN}}^{(1)}`$ $`=`$ $`F_g^{(1)}(0)F_H^{(0)}(0)f_{H/g}^{(1)}(0)F_{\overline{H}}^{(0)}(0)f_{\overline{H}/g}^{(1)}(0),`$
$`\widehat{F}_{H,\text{ZM-VFN}}^{(1)}`$ $`=`$ $`F_H^{(1)}(0)F_H^{(0)}(0)f_{H/H}^{(1)}(0).`$ (22)
We should note that in a calculation with $`M_H=0`$, there are infrared divergences. Thus the calculations are performed in $`42ϵ`$ dimensions of space-time and $`1/ϵ`$ poles appear. The pole terms cancel in Eq. (22). Instead of using dimensional regulation and taking $`ϵ0`$, one could use an infrared regulator mass $`m`$ and let $`m0`$.

* Property 2. For $`Q`$ of order $`M_H`$, the structure function is that of the fixed-flavor-number prescription, up to corrections of order $`\alpha _s^2`$.

In the fixed-flavor-number prescription we have

$$F=F_g^{(1)}(M_H)f_{g/p}+𝒪(\alpha _s^2)+𝒪(\mathrm{\Lambda }^2/Q^2).$$
(23)
Property 2 follows because when the factorization scale, $`\mu =Q`$, is of order $`M`$, the heavy quark distribution function is given by the perturbative formula

$$f_{H/p}=f_{H/g}^{(1)}f_{g/p}+𝒪(\alpha _s^2).$$
(24)
First of all, since $`f_{H/p}`$ is of order $`\alpha _s`$, the terms

$$\widehat{F}_H^{(1)}f_{H/p}+\widehat{F}_{\overline{H}}^{(1)}f_{\overline{H}/p}$$
(25)
in Eq. (15) are of order $`\alpha _s^2`$ and can be dropped. This leaves

$`F`$ $`=`$ $`F_H^{(0)}(M_H)f_{H/p}+F_{\overline{H}}^{(0)}(M_H)f_{\overline{H}/p}`$ (26)
$`+`$ $`[F_g^{(1)}(M_H)F_H^{(0)}(M_H)f_{H/g}^{(1)}(M_H)F_{\overline{H}}^{(0)}(M_H)f_{\overline{H}/g}^{(1)}(M_H)]f_{g/p}`$
$`+`$ $`𝒪(\alpha _s^2)+𝒪(\mathrm{\Lambda }^2/Q^2).`$
Inserting Eq. (24) into Eq. (26), we obtain Eq. (23).

We can summarize properties 1 and 2 by saying that the ACOT prescription interpolates between the zero-mass prescription for $`M_HQ`$ and the fixed-flavor-number prescription for $`M_HQ`$.

What about the simplified ACOT prescription? We use the same formulas as for the ACOT prescription, but set $`M_H=0`$ in the hard scattering functions for heavy quarks:

$`\widehat{F}_{H,\text{S-ACOT}}^{(0)}`$ $`=`$ $`F_H^{(0)}(0),`$
$`\widehat{F}_{g,\text{S-ACOT}}^{(1)}`$ $`=`$ $`F_g^{(1)}(M_H)F_H^{(0)}(0)f_{H/g}^{(1)}(M_H)F_{\overline{H}}^{(0)}(0)f_{\overline{H}/g}^{(1)}(M_H),`$
$`\widehat{F}_{H,\text{S-ACOT}}^{(1)}`$ $`=`$ $`F_H^{(1)}(0)F_H^{(0)}(0)f_{H/H}^{(1)}(0).`$ (27)
Repeating the derivation just given, we see that properties 1 and 2 hold for the S-ACOT prescription. That is, the S-ACOT prescription also interpolates between the zero-mass prescription for $`M_HQ`$ and the fixed-flavor-number prescription for $`M_HQ`$. However the S-ACOT prescription has the advantage that the functions $`\widehat{F}_H`$ are easier to calculate.

## 6 Comparison of different prescriptions

In this section, we investigate how well the matching properties among the different prescriptions work. In the plots presented, the masses and couplings have been fixed to be consistent with the values used in the CTEQ4L/M fits, i.e. $`m_c=1.6`$ GeV, $`m_b=5`$ GeV . Below the bottom threshold $`\mathrm{n}_\mathrm{f}=4`$ active flavors are used for $`\alpha _s`$, and the scale has been chosen as $`\mu =Q`$.

### 6.1 Parton distribution matching

In the threshold region, the structure function calculated in the S-ACOT prescription matches that calculated in the FFN prescription, in which the heavy quark does not appear as a parton. This matching, Property 2, results from the fact that the heavy quark distribution function $`f_{H/p}(x,\mu )`$, matches the approximate function

$$\stackrel{~}{f}_{H/p}(\mu )=f_{H/g}^{(1)}(\mu )f_{g/p}(\mu ).$$
(28)
Both functions vanish at $`\mu =M_H`$ and, for $`\mu >M_H`$, the difference between them is of order $`\alpha _s^2`$. In this subsection, we study the quality of this matching.

In Figs. 1–3 we plot $`f_{H/p}(x,\mu )`$ and $`\stackrel{~}{f}_{H/p}(x,\mu )`$ for the case of the charm quark. In Fig. 1, we use CTEQ4L parton distributions, which are based on lowest order evolution. Fig. 1 reveals that the matching works very well when the order $`\alpha _s^1`$ evolution kernel is used for $`f_{H/p}(x,\mu )`$ and the order $`\alpha _s^1`$ perturbative expression (Eq. (28)) is used for $`\stackrel{~}{f}`$. In Fig. 2, we use CTEQ4M parton distributions, which are based on NLO evolution. We observe a mismatch because the full NLO evolution kernel is used for $`f_{H/p}`$ while the order $`\alpha _s^1`$ perturbative expansion is used for $`\stackrel{~}{f}`$. The difference between $`f`$ and $`\stackrel{~}{f}`$ is of order $`\alpha _s^2`$, but this difference is numerically quite large. Finally, we see in Fig. 3 that a close match is restored when the NLO evolution kernel is used for $`f`$ while $`\stackrel{~}{f}`$ is defined using a modified version of Eq. (28) in which we replace the leading order $`gH`$ evolution kernel in Eq. (18) by the next-to-leading order $`gH`$ evolution kernel.222For the NLO perturbative expansion, we have included the $`\alpha _s^2`$ splitting kernels, $`P_{j/i}^{(2)}(\xi )`$. For simplicity, we have ignored iterated terms such as $`P_{k/j}^{(1)}P_{j/i}^{(1)}`$ which contribute as $`\mathrm{ln}^2(\mu /M_H)`$, and hence only play a role away from threshold.

We can draw two conclusions. First, the threshold matching discussed in the previous section will work order by order for the perturbation expansion.333 At order $`\alpha _s^2`$, Property 2 is still preserved even though the matching conditions on the parton distribution functions are modified. The $`\alpha _s^2`$ matching conditions shift the evolved parton distribution functions $`f_{H/p}`$, but they also shift the “perturbative” parton distribution functions $`\stackrel{~}{f}_{H/p}`$ leaving the difference unchanged up to order $`\alpha _s^3`$.  Second, the order $`\alpha _s^2`$ terms in the evolution kernel are quite large, so that the leading order calculations illustrated in this paper may not be sufficient for obtaining accurate predictions.444 One might ask whether retaining the heavy quark mass in the kinematics (in the spirit of the “slow-rescaling” correction) might prove beneficial. The answer is that as long as the order of the parton distribution functions and the hard scattering coefficients are matched, this simply amounts to a shuffling of $`M_H/Q`$ terms between the quark and gluon initiated terms, cf., Eq. (6).

### 6.2 Structure function matching

In this subsection, we examine predictions for $`F_2^c(x,Q)`$, which we define here to be the contribution to $`F_2(x,Q)`$ from graphs in which the current couples to a charm quark. We compare $`F_2^c(x,Q)`$ calculated with the S-ACOT prescription at order $`\alpha _s`$ with that calculated with the original ACOT prescription, the ZM-VFN prescription in which the charm quark can appear as a parton but has zero mass, and the FFN prescription in which the charm quark has its proper mass but does not appear as a parton. For simplicity, we take $`\mu =Q`$. In our calculations, we evaluate the hard scattering coefficients $`\widehat{F}`$ at order $`\alpha _s`$. Thus for the ZM-VFN, ACOT, and S-ACOT prescriptions, the terms displayed in Eq. (15) are included. The functions $`\widehat{F}`$ are given by Eq. (22) for ZM-VFN, by Eqs. (16,17,19) for ACOT, and by Eq. (27) for S-ACOT. For the FFN prescription, there is only one term at order $`\alpha _s`$, as displayed in Eq. (23).

In Fig. 4 we show $`F_2^c(x,Q)`$ as a function of $`Q`$ for $`x=0.1`$. Then in Fig. 5 we show $`F_2^c(x,Q)`$ as a function of $`Q`$ for $`x=0.001`$. In each case we display results using both the CTEQ4L and CTEQ4M parton distributions.

When we use the CTEQ4L parton distributions, we notice that there is a close match between the S-ACOT result and the FFN result near $`Q=m_c`$. Based on the results of the previous subsection, we expect this matching to be degraded when we use CTEQ4M parton distributions because of the important role played by the order $`\alpha _s^2`$ term in the evolution kernel that is not matched in the lowest order calculation of the hard scattering function. This degradation is seen in the figures.

In the asymptotic regime, $`QM`$, we find the S-ACOT result approximates the ZM-VFN result, as expected.

We observe that the ACOT and S-ACOT prescriptions are effectively identical throughout the kinematic range. There is a slight difference in the threshold region, but this is small in comparison to the size of the $`\mu `$-variation (not shown). Hence the difference between the ACOT and S-ACOT results is of no physical consequence. The fact that the ACOT and S-ACOT prescriptions match extremely well throughout the full kinematic range provides explicit numerical verification that the S-ACOT prescription fully contains the physics.

## 7 Conclusions and outlook

We have performed a numerical study of different prescriptions for dealing with the quark mass in heavy quark leptoproduction. We have seen that the simplest prescription , S-ACOT, is numerically equivalent to the earlier ACOT prescription .

The S-ACOT prescription is extensible order by order in $`\alpha _s`$. At $`𝒪(\alpha _s^1)`$ we already find a significant simplification in the S-ACOT prescription as compared with the ACOT prescription. We expect that this simplification becomes even more dramatic at higher orders.

We have not attempted to implement the S-ACOT prescription at order $`\alpha _s^2`$, but we note that the NLO corrections in the the FFN prescription have been calculated , and that the leading (collinear) logarithms of the type $`\alpha _s^i\mathrm{log}^i(Q^2/M_H^2)`$ have been extracted in analytic form . Thus the S-ACOT subtraction term can be constructed as well.

Finally, we emphasize that the choice of a prescription for dealing with quark masses in the hard scattering coefficients for deeply inelastic scattering is a separate issue from the choice of definition of the parton distribution functions. For all of the prescriptions discussed here, one uses the standard   MS  definition of parton distributions.

## Acknowledgments

We thank J.C. Collins, R.J. Scalise, R.S. Thorne, and W.-K. Tung for valuable discussions and we thank S. Kretzer for making his NLO code available to us. We thank the CERN Theory Division, Fermilab Theory Group, and the University of Oregon for their kind hospitality during the period in which part of this research was carried out. This work is supported by the U.S. Department of Energy, the Lightner-Sams Foundation, and in part by the EU TMR contract FMRX-CT98-0194 (DG 12 - MIHT).