text:
# Azimuthal Asymmetries in Heavy Quark Leptoproduction as a Test of pQCD

## I Introduction

In the framework of perturbative QCD, the basic spin-averaged characteristics of heavy flavor hadro-, photo- and electroproduction are known exactly up to the next-to-leading order (NLO). During the last ten years, these NLO results have been widely used for a phenomenological description of available data (for a review see ). At the same time, the key question remains open: How to test the applicability of QCD at fixed order to heavy quark production? The problem is twofold. On the one hand, the NLO corrections are large; they increase the leading order (LO) predictions for both charm and bottom production cross sections by approximately a factor of two. For this reason, one could expect that higher-order corrections, as well as nonperturbative contributions, can be essential, especially for the $`c`$-quark case. On the other hand, it is very difficult to compare pQCD predictions for spin-averaged cross sections with experimental data directly, without additional assumptions, because of a high sensitivity of the theoretical calculations to standard uncertainties in the input QCD parameters. The total uncertainties associated with the unknown values of the heavy quark mass, $`m`$, the factorization and renormalization scales, $`\mu _F`$ and $`\mu _R`$, $`\mathrm{\Lambda }_{QCD}`$ and the parton distribution functions are so large that one can only estimate the order of magnitude of the pQCD predictions for total cross sections at fixed target energies .

In recent years, the role of higher-order corrections has been extensively investigated in the framework of the soft gluon resummation formalism. For a review see Ref.. Soft gluon (or threshold) resummation is based on the factorization properties of the cross section near the partonic threshold and makes it possible to resum to all orders in $`\alpha _s`$ the leading (Sudakov double) logarithms (LL) and the next-to-leading ones (NLL) . Formally resummed cross sections are ill-defined due to the Landau pole contribution, and a few prescriptions have been proposed to avoid the renormalon ambiguities . Unfortunately, numerical predictions for the heavy quark production cross sections can depend significantly on the choice of resummation prescription . Another open question, also closely related to convergence of the perturbative series, is the role of subleading logarithms which are not, in principle, under control of the resummation procedure .

For this reason, it is of special interest to study those observables that are well-defined in pQCD. A nontrivial example of such an observable is proposed in , where the single spin asymmetry (SSA) in charm and bottom production by linearly polarized photons, $`\gamma ^{}+NQ+X[\overline{Q}]`$, was calculated.The well-known examples are the shapes of differential cross sections of heavy flavor production which are sufficiently stable under radiative corrections. It was shown that, contrary to the production cross section, the single spin asymmetry in heavy flavor photoproduction is quantitatively well defined in pQCD: it is stable, both parametrically and perturbatively, and insensitive to nonperturbative corrections. Therefore, measurements of this asymmetry would provide an ideal test of pQCD. As was shown in Ref. , the SSA in open charm photoproduction can be measured with an accuracy of about ten percent in the approved E160/E161 experiments at SLAC  using the inclusive spectra of secondary (decay) leptons.

In the present paper we continue the studies of perturbatively stable observables and calculate the radiative and nonperturbative corrections to the azimuthal asymmetry (AA) in heavy quark leptoproduction:

$$l(\mathrm{})+N(p)l(\mathrm{}q)+Q(p_Q)+X[\overline{Q}](p_X).$$
(1)
In the case of unpolarized initial states and neglecting the contribution of $`Z`$boson, the cross section of the reaction (1) may be written as

$`{\displaystyle \frac{\text{d}\sigma _{lN}}{\text{d}x\text{d}Q^2\text{d}\phi }}={\displaystyle \frac{\alpha _{em}}{(2\pi )^2}}{\displaystyle \frac{1}{xQ^2}}`$ $`\{[1+(1y)^2]\sigma _T(x,Q^2)+2(1y)\sigma _L(x,Q^2)`$ (3)
$`+2(1y)\sigma _A(x,Q^2)\mathrm{cos}2\phi +(2y)\sqrt{1y}\sigma _I(x,Q^2)\mathrm{cos}\phi \}.`$
The kinematic variables are defined by

$`\overline{S}=\left(\mathrm{}+p\right)^2,`$ $`Q^2=q^2,`$ $`x={\displaystyle \frac{Q^2}{2pq}},`$ (4)
$`y={\displaystyle \frac{pq}{p\mathrm{}}},`$ $`Q^2=xy\overline{S},`$ $`\rho ={\displaystyle \frac{4m^2}{\overline{S}}}.`$ (5)
In Eq. (3), $`\sigma _T(\sigma _L)`$ is the usual $`\gamma ^{}N`$ cross section describing heavy quark production by a transverse (longitudinal) virtual photon. The third cross section, $`\sigma _A`$, comes about from interference between transverse states and is responsible for the SSA which occurs in real photoproduction using linearly polarized photons . The fourth cross section, $`\sigma _I`$, originates from interference between longitudinal and transverse components . In the nucleon rest frame, the azimuth $`\phi `$ is the angle between the lepton scattering plane and the heavy quark production plane, defined by the exchanged photon and the detected quark $`Q`$ (see Fig. 1). The covariant definition of $`\phi `$ is

$`\mathrm{cos}\phi `$ $`=`$ $`{\displaystyle \frac{rn}{\sqrt{r^2}\sqrt{n^2}}},\mathrm{sin}\phi ={\displaystyle \frac{Q^2\sqrt{1/x^2+4m_N^2/Q^2}}{2\sqrt{r^2}\sqrt{n^2}}}n\mathrm{},`$ (6)
$`r^\mu `$ $`=`$ $`\epsilon ^{\mu \nu \alpha \beta }p_\nu q_\alpha \mathrm{}_\beta ,n^\mu =\epsilon ^{\mu \nu \alpha \beta }p_\nu q_\alpha p_{Q\beta }.`$ (7)
In Eqs. (5) and (6), $`m`$ and $`m_N`$ are the masses of the heavy quark and the target, respectively.

In leading order pQCD, the $`\mathrm{cos}\phi `$ dependence vanishes, $`\sigma _I^{\mathrm{Born}}(x,Q^2)=0`$. For this reason, in this paper we restrict ourselves to the azimuthal asymmetry, $`A(\rho ,x,Q^2)`$, associated with the $`\mathrm{cos}2\phi `$ distribution:

$$A(\rho ,x,Q^2)=\frac{\text{d}^3\sigma _{lN}(\phi =0)+\text{d}^3\sigma _{lN}(\phi =\pi )2\text{d}^3\sigma _{lN}(\phi =\pi /2)}{\text{d}^3\sigma _{lN}(\phi =0)+\text{d}^3\sigma _{lN}(\phi =\pi )+2\text{d}^3\sigma _{lN}(\phi =\pi /2)}=\frac{\epsilon \sigma _A(x,Q^2)}{\sigma _T(x,Q^2)+\epsilon \sigma _L(x,Q^2)},$$
(8)
where $`\epsilon ={\displaystyle \frac{2(1y)}{1+(1y)^2}}`$ and $`\text{d}^3\sigma _{lN}(\phi ){\displaystyle \frac{\text{d}^3\sigma _{lN}}{\text{d}x\text{d}Q^2\text{d}\phi }}(\rho ,x,Q^2,\phi )`$. Note that the asymmetry defined by Eq. (8) is simply related to the mean value of $`\mathrm{cos}2\phi `$:

$$A(\rho ,x,Q^2)=2\mathrm{cos}2\phi (\rho ,x,Q^2),\mathrm{cos}2\phi (\rho ,x,Q^2)=\frac{\underset{0}{\overset{2\pi }{}}\text{d}\phi \mathrm{cos}2\phi {\displaystyle \frac{\text{d}^3\sigma _{lN}}{\text{d}x\text{d}Q^2\text{d}\phi }}(\rho ,x,Q^2,\phi )}{\underset{0}{\overset{2\pi }{}}\text{d}\phi {\displaystyle \frac{\text{d}^3\sigma _{lN}}{\text{d}x\text{d}Q^2\text{d}\phi }}(\rho ,x,Q^2,\phi )}.$$
(9)
In this paper, we calculate the NLO corrections to the $`\mathrm{cos}2\phi `$ asymmetry to the next-to-leading logarithmic accuracy (so-called soft-gluon approximation). Also we analyze the nonperturbative contributions to the AA due to the gluon transverse motion in the target. Our main results can be formulated as follows:

* The azimuthal asymmetry defined by Eq. (8) is of leading twist; at energies sufficiently above the production threshold, it is predicted to be about $`15\%`$ for both charm and bottom quark production.
* The soft-gluon approximation provides a good description of the available exact NLO results on leptoproduction in the region of relatively low virtualities, $`Q^2m^2`$; when $`Q^2m^2`$, the quality of the NLL approximation becomes worse.
* Contrary to the production cross sections, the $`\mathrm{cos}2\phi `$ asymmetry in azimuthal distributions of both charm and bottom quark is practically insensitive to radiative corrections at $`Q^2m^2`$. This implies that large soft-gluon contributions to the $`\phi `$-dependent and $`\phi `$-integrated cross sections cancel each other in Eq. (8) with a good accuracy.
* pQCD predictions for the $`\mathrm{cos}2\phi `$ asymmetry are parametrically stable; to within few percent, they are insensitive to standard uncertainties in the QCD input parameters: $`\mu _R`$, $`\mu _F`$, $`\mathrm{\Lambda }_{QCD}`$ and in the gluon distribution function.
* Nonperturbative corrections to the $`b`$-quark azimuthal asymmetry due to the gluon transverse motion in the target are negligible. Because of the smallness of the $`c`$-quark mass, the analogous corrections to $`A(\rho ,x,Q^2)`$ in the charm case are larger; they are of the order of 20% at $`Q^2m^2`$.

We conclude that, in contrast with the production cross sections, the $`\mathrm{cos}2\phi `$ asymmetry in heavy quark leptoproduction, $`A(\rho ,x,Q^2)`$, is an observable quantitatively well defined in pQCD: it is stable, both parametrically and perturbatively, and insensitive (in the case of bottom production) to nonperturbative corrections. Measurements of the AA in bottom production would provide an ideal test of the conventional parton model based on pQCD.

Concerning the experimental aspects, azimuthal asymmetries in bottom production can, in principle, be measured at HERA using the angular distributions of secondary (decay) leptons . AAs in charm leptoproduction can also be measured in the COMPASS and HERMES experiments. Due to the relatively low $`c`$-quark mass, data on the $`D`$-meson azimuthal distributions would make it possible to clarify the role of subleading twist contributions.

The paper is organized as follows. In Section II we analyze the LO and NLO parton level predictions for $`\phi `$-dependent leptoproduction of heavy flavor in the single-particle inclusive kinematics. We check the quality of the soft-gluon approach against available exact results and discuss the region of applicability of the NLL approximation. Hadron level predictions for $`A(\rho ,x,Q^2)`$ are given in Section III. We consider in detail the pQCD contributions and nonperturbative corrections to the $`\mathrm{cos}2\phi `$ asymmetry at the HERMES, SLAC, COMPASS and HERA energies.

## II  Partonic Cross Sections

### A Born level predictions

At leading order, $`𝒪(\alpha _{em}\alpha _s)`$, the only partonic subprocess which is responsible for heavy quark leptoproduction is the two-body photon-gluon fusion:

$$\gamma ^{}(q)+g(k_g)Q(p_Q)+\overline{Q}(p_{\stackrel{\mathrm{\_}}{Q}}).$$
(10)
The $`\gamma ^{}g`$ cross sections, $`\widehat{\sigma }_k^{\mathrm{Born}}`$ ($`k=T,L,A,I`$), corresponding to the Born diagrams are :

$`\widehat{\sigma }_T^{\mathrm{Born}}(\widehat{\rho },\widehat{x})`$ $`=`$ $`{\displaystyle \frac{\pi e_Q^2\alpha _{em}\alpha _s}{2m^2}}\widehat{\rho }\{[(1\widehat{x})^2+\widehat{x}^2+\widehat{\rho }(1\widehat{x}\widehat{\rho }/2)]\mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }}`$ (12)
$`[(1\widehat{x})^2\widehat{x}(23\widehat{x})+(1\widehat{x})\widehat{\rho }]\beta \},`$
$`\widehat{\sigma }_L^{\mathrm{Born}}(\widehat{\rho },\widehat{x})`$ $`=`$ $`{\displaystyle \frac{\pi e_Q^2\alpha _{em}\alpha _s}{m^2}}\widehat{\rho }\widehat{x}\left\{\widehat{\rho }\mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }}+2(1\widehat{x})\beta \right\},`$ (13)
$`\widehat{\sigma }_A^{\mathrm{Born}}(\widehat{\rho },\widehat{x})`$ $`=`$ $`{\displaystyle \frac{\pi e_Q^2\alpha _{em}\alpha _s}{2m^2}}\widehat{\rho }^2\left\{(12\widehat{x}\widehat{\rho }/2)\mathrm{ln}{\displaystyle \frac{1+\beta }{1\beta }}(1\widehat{x})(12\widehat{x}/\widehat{\rho })\beta \right\},`$ (14)
$`\widehat{\sigma }_I^{\mathrm{Born}}(\widehat{\rho },\widehat{x})`$ $`=`$ $`0.`$ (15)
In Eqs. (13), $`e_Q`$ is the quark charge in units of electromagnetic coupling constant and we use the following definition of partonic kinematic variables:

$`s=\left(q+k_g\right)^2,`$ $`\widehat{x}={\displaystyle \frac{Q^2}{2qk_g}},`$ (16)
$`\beta =\sqrt{1{\displaystyle \frac{4m^2}{s}}},`$ $`\widehat{\rho }={\displaystyle \frac{4m^2}{s+Q^2}}.`$ (17)
Note that the $`\mathrm{cos}\phi `$ dependence vanishes due to the $`Q\overline{Q}`$ symmetry which, at leading order, requires invariance under $`\phi \phi +\pi `$ .

The hadron level cross sections, $`\sigma _k(x,Q^2)`$ ($`k=T,L,A,I`$), have the form

$$\sigma _k(x,\xi )=\underset{x+4x/\xi }{\overset{1}{}}\text{d}zg(z,\mu _F)\widehat{\sigma }_k(\frac{4x}{z\xi },\frac{x}{z}),k_g=zp,\xi =\frac{Q^2}{m^2},$$
(18)
where $`g(z,\mu _F)`$ describes gluon density in a nucleon $`N`$ evaluated at a factorization scale $`\mu _F`$. The partonic cross sections, $`\widehat{\sigma }_k`$, are functions of $`\widehat{\rho }`$ and $`\widehat{x}`$ defined by Eq. (17).

### B Soft-gluon corrections at NLO

To take into account the NLO contributions, one needs to calculate the virtual $`𝒪(\alpha _{em}\alpha _s^2)`$ corrections to the Born process (10) and the real gluon emission:

$$\gamma ^{}(q)+g(k_g)Q(p_Q)+\overline{Q}(p_{\stackrel{\mathrm{\_}}{Q}})+g(p_g).$$
(19)
The partonic invariants describing the single-particle inclusive (1PI) kinematics are

$`s^{}=2qk_g=s+Q^2=zS^{},`$ $`t_1=\left(k_gp_Q\right)^2m^2=zT_1,`$ (20)
$`s_4=s^{}+t_1+u_1,`$ $`u_1=\left(qp_Q\right)^2m^2=U_1,`$ (21)
where $`s_4`$ measures the inelasticity of the reaction (19). The corresponding 1PI hadron level variables describing the reaction (1) are

$`S^{}=2qp=S+Q^2,`$ $`T_1=\left(pp_Q\right)^2m^2,`$ (22)
$`S_4=S^{}+T_1+U_1,`$ $`U_1=\left(qp_Q\right)^2m^2.`$ (23)
We neglect the photon-(anti)quark fusion subprocesses. This is justified as their contributions vanish at LO and are small at NLO .

The exact NLO calculations of the unpolarized heavy quark production in $`\gamma g`$ , $`\gamma ^{}g`$ , and $`gg`$  collisions show that, near the partonic threshold, a strong logarithmic enhancement of the cross sections takes place in the collinear, $`\stackrel{}{p}_{g,T}`$ $`0`$, and soft, $`\stackrel{}{p}_g0`$, limits. This threshold (or soft-gluon) enhancement has universal nature in the perturbation theory and originates from incomplete cancellation of the soft and collinear singularities between the loop and the bremsstrahlung contributions. Large leading and next-to-leading threshold logarithms can be resummed to all orders of perturbative expansion using the appropriate evolution equations . The analytic results for the resummed cross sections are ill-defined due to the Landau pole in the coupling strength $`\alpha _s`$. However, if one considers the obtained expressions as generating functionals of the perturbative theory and re-expands them at fixed order in $`\alpha _s`$, no divergences associated with the Landau pole are encountered.

Soft-gluon resummation for the photon-gluon fusion has been performed in Ref. and checked in Refs.. To NLL accuracy, the perturbative expansion for the partonic cross sections, d$`{}_{}{}^{2}\widehat{\sigma }_{k}^{}/`$d$`t_1`$d$`u_1`$ ($`k=T,L,A,I`$), can be written in a factorized form as

$$s^2\frac{\text{d}^2\widehat{\sigma }_k}{\text{d}t_1\text{d}u_1}(s^{},t_1,u_1)=B_k^{\text{Born}}(s^{},t_1,u_1)\left\{\delta (s^{}+t_1+u_1)+\underset{n=1}{\overset{\mathrm{}}{}}\left(\frac{\alpha _sC_A}{\pi }\right)^nK^{(n)}(s^{},t_1,u_1)\right\},$$
(24)
with the Born level distributions $`B_k^{\text{Born}}`$ given by

$`B_T^{\text{Born}}(s^{},t_1,u_1)`$ $`=`$ $`\pi e_Q^2\alpha _{em}\alpha _s\left[{\displaystyle \frac{t_1}{u_1}}+{\displaystyle \frac{u_1}{t_1}}+4\left({\displaystyle \frac{s}{s^{}}}{\displaystyle \frac{m^2s^{}}{t_1u_1}}\right)\left({\displaystyle \frac{s^{}(m^2Q^2/2)}{t_1u_1}}+{\displaystyle \frac{Q^2}{s^{}}}\right)\right],`$ (25)
$`B_L^{\text{Born}}(s^{},t_1,u_1)`$ $`=`$ $`\pi e_Q^2\alpha _{em}\alpha _s\left[{\displaystyle \frac{8Q^2}{s^{}}}\left({\displaystyle \frac{s}{s^{}}}{\displaystyle \frac{m^2s^{}}{t_1u_1}}\right)\right],`$ (26)
$`B_A^{\text{Born}}(s^{},t_1,u_1)`$ $`=`$ $`\pi e_Q^2\alpha _{em}\alpha _s\left[4\left({\displaystyle \frac{s}{s^{}}}{\displaystyle \frac{m^2s^{}}{t_1u_1}}\right)\left({\displaystyle \frac{m^2s^{}}{t_1u_1}}+{\displaystyle \frac{Q^2}{s^{}}}\right)\right],`$ (27)
$`B_I^{\text{Born}}(s^{},t_1,u_1)`$ $`=`$ $`\pi e_Q^2\alpha _{em}\alpha _s\left[4\sqrt{Q^2}\left({\displaystyle \frac{t_1u_1s}{s^2}}m^2\right)^{1/2}{\displaystyle \frac{t_1u_1}{t_1u_1}}\left(1{\displaystyle \frac{2Q^2}{s^{}}}{\displaystyle \frac{2m^2s^{}}{t_1u_1}}\right)\right].`$ (28)
Note that the functions $`K^{(n)}(s^{},t_1,u_1)`$ in Eq. (24) originate from the collinear and soft limits. Since the azimuthal angle $`\phi `$ is the same for both $`\gamma ^{}g`$ and $`Q\overline{Q}`$ center-of-mass systems in these limits, the functions $`K^{(n)}(s^{},t_1,u_1)`$ are also the same for all $`k=T,L,A,I`$. At NLO, the soft-gluon corrections to NLL accuracy in the $`\overline{\text{MS}}`$ scheme are

$`K^{(1)}(s^{},t_1,u_1)`$ $`=`$ $`2\left[{\displaystyle \frac{\mathrm{ln}\left(s_4/m^2\right)}{s_4}}\right]_+\left[{\displaystyle \frac{1}{s_4}}\right]_+\left\{1+\mathrm{ln}\left({\displaystyle \frac{u_1}{t_1}}\right)\left(1{\displaystyle \frac{2C_F}{C_A}}\right)\left(1+\text{Re}L_\beta \right)+\mathrm{ln}\left({\displaystyle \frac{\mu ^2}{m^2}}\right)\right\}`$ (30)
$`+\delta (s_4)\mathrm{ln}\left({\displaystyle \frac{u_1}{m^2}}\right)\mathrm{ln}\left({\displaystyle \frac{\mu ^2}{m^2}}\right),`$
where we use $`\mu =\mu _F=\mu _R`$. In Eq. (30), $`C_A=N_c`$ and $`C_F=(N_c^21)/(2N_c)`$, where $`N_c`$ is number of colors, while $`L_\beta =(12m^2/s)\{\mathrm{ln}[(1\beta )/(1+\beta )]+`$i$`\pi \}`$. The single-particle inclusive ”plus“ distributions are defined by

$$\left[\frac{\mathrm{ln}^l\left(s_4/m^2\right)}{s_4}\right]_+=\underset{ϵ0}{lim}\left\{\frac{\mathrm{ln}^l\left(s_4/m^2\right)}{s_4}\theta (s_4ϵ)+\frac{1}{l+1}\mathrm{ln}^{l+1}\left(\frac{ϵ}{m^2}\right)\delta (s_4)\right\}.$$
(31)
For any sufficiently regular test function $`h(s_4)`$, Eq. (31) gives

$$\underset{0}{\overset{s_4^{\mathrm{max}}}{}}\text{d}s_4h(s_4)\left[\frac{\mathrm{ln}^l\left(s_4/m^2\right)}{s_4}\right]_+=\underset{0}{\overset{s_4^{\mathrm{max}}}{}}\text{d}s_4\left[h(s_4)h(0)\right]\frac{\mathrm{ln}^l\left(s_4/m^2\right)}{s_4}+\frac{1}{l+1}h(0)\mathrm{ln}^{l+1}\left(s_4^{\mathrm{max}}/m^2\right).$$
(32)
In Eq. (30) , we have preserved the NLL terms for the scale-dependent logarithms too. We have checked that the results (27) and (30) agree to NLL accuracy with the exact $`𝒪(\alpha _{em}\alpha _s^2)`$ calculations of the photon-gluon cross sections $`\widehat{\sigma }_T`$ and $`\widehat{\sigma }_L`$ given in Ref..

To perform a numerical investigation of the results (27) and (30), it is convenient to introduce for the fully inclusive (integrated over $`t_1`$ and $`u_1`$) cross sections, $`\widehat{\sigma }_k`$ ($`k=T,L,A,I`$), the dimensionless coefficient functions $`c_k^{(n,l)}`$,

$$\widehat{\sigma }_k(\eta ,\xi ,\mu ^2)=\frac{e_Q^2\alpha _{em}\alpha _s(\mu ^2)}{m^2}\underset{n=0}{\overset{\mathrm{}}{}}\left(4\pi \alpha _s(\mu ^2)\right)^n\underset{l=0}{\overset{n}{}}c_k^{(n,l)}(\eta ,\xi )\mathrm{ln}^l\left(\frac{\mu ^2}{m^2}\right),$$
(33)
where the variable $`\eta `$ measures the distance to the partonic threshold:

$$\eta =\frac{s}{4m^2}1,\xi =\frac{Q^2}{m^2}.$$
(34)
Concerning the NLO scale-independent coefficient functions, only $`c_T^{(1,0)}`$ and $`c_L^{(1,0)}`$ are known exactly . As to the $`\mu `$-dependent coefficients, they can by calculated explicitly using the renormalization group equation:

$$\frac{\text{d}\widehat{\sigma }_k(s^{},Q^2,\mu ^2)}{\text{d}\mathrm{ln}\mu ^2}=\underset{z_{\mathrm{min}}}{\overset{1}{}}\text{d}z\widehat{\sigma }_k(zs^{},Q^2,\mu ^2)P_{gg}(z),$$
(35)
where $`z_{\mathrm{min}}=(4m^2+Q^2)/s^{}`$, $`\widehat{\sigma }_k(s^{},Q^2,\mu )`$ are the cross sections resummed to all orders in $`\alpha _s`$ and $`P_{gg}(z)`$ is the corresponding (resummed) Altarelli-Parisi gluon-gluon splitting function. Expanding Eq. (35) in $`\alpha _s`$, one can find

$$c_k^{(1,1)}(s^{},\xi )=\frac{1}{4\pi ^2}\underset{z_{\mathrm{min}}}{\overset{1}{}}\text{d}z\left[b_2\delta (1z)P_{gg}^{(0)}(z)\right]c_k^{(0,0)}(zs^{},\xi ),$$
(36)
where $`b_2=(11C_A2n_f)/12`$ is the first coefficient of the $`\beta (\alpha _s)`$-function expansion and $`n_f`$ is the number of active quark flavors. The one-loop gluon splitting function is :

$$P_{gg}^{(0)}(z)=\underset{ϵ0}{lim}\left\{\left(\frac{z}{1z}+\frac{1z}{z}+z(1z)\right)\theta (1zϵ)+\delta (1z)\mathrm{ln}ϵ\right\}C_A+b_2\delta (1z).$$
(37)
With Eq. (36) in hand, we are able to check the quality of the NLL approximation against exact answers. In Figs. 2 and 3 we plot the functions $`c_T^{(n,l)}(\eta ,\xi )`$ and $`c_A^{(n,l)}(\eta ,\xi )`$ $`(n,l=0,1)`$ at $`\xi =10^2`$ and $`\xi =3.16`$, respectively. Predictions of the NLL approximation (30) are given by dotted curves. The available exact results are given by solid lines. One can see a reasonable agreement up to energies $`\eta 2`$. As to the $`Q^2`$-dependence, we have found that the soft-gluon approach reproduces satisfactorily the exact results at $`\xi 1.`$ At high values of $`\xi `$, $`Q^2m^2`$, the quality of the NLL approximation becomes worse.Our analysis shows that the same situation takes also place for the energy and $`Q^2`$ behavior of the functions $`c_L^{(1,l)}(\eta ,\xi )`$, $`l=0,1`$. We do not give corresponding plots for $`c_L^{(1,l)}(\eta ,\xi )`$ because the contribution of the longitudinal cross section to the $`\mathrm{cos}2\phi `$ asymmetry is small numerically.

## III Hadron Level Results

### A pQCD predictions

Let us now analyze the impact of the approximate NLO perturbative corrections on the AA at hadron level. We will consider the parameters $`A(Q^2)`$, $`A(x)`$ and $`A(y)`$,

$`A(Q^2)={\displaystyle \frac{2\underset{0}{\overset{2\pi }{}}\text{d}\phi \mathrm{cos}2\phi {\displaystyle \frac{\text{d}^2\sigma _{lN}}{\text{d}Q^2\text{d}\phi }}(\rho ,Q^2,\phi )}{\underset{0}{\overset{2\pi }{}}\text{d}\phi {\displaystyle \frac{\text{d}^2\sigma _{lN}}{\text{d}Q^2\text{d}\phi }}(\rho ,Q^2,\phi )}},A(x)={\displaystyle \frac{2\underset{0}{\overset{2\pi }{}}\text{d}\phi \mathrm{cos}2\phi {\displaystyle \frac{\text{d}^2\sigma _{lN}}{\text{d}x\text{d}\phi }}(\rho ,x,\phi )}{\underset{0}{\overset{2\pi }{}}\text{d}\phi {\displaystyle \frac{\text{d}^2\sigma _{lN}}{\text{d}x\text{d}\phi }}(\rho ,x,\phi )}},`$

$$A(y)=\frac{2\underset{0}{\overset{2\pi }{}}\text{d}\phi \mathrm{cos}2\phi {\displaystyle \frac{\text{d}^2\sigma _{lN}}{\text{d}y\text{d}\phi }}(\rho ,y,\phi )}{\underset{0}{\overset{2\pi }{}}\text{d}\phi {\displaystyle \frac{\text{d}^2\sigma _{lN}}{\text{d}y\text{d}\phi }}(\rho ,y,\phi )},$$
(38)
which describe the dependence of the $`\mathrm{cos}2\phi `$ asymmetry on $`Q^2`$, Bjorken $`x`$ and $`y`$, respectively. Unless otherwise stated, the CTEQ5M  parametrization of the gluon distribution function is used. The default values of the charm and bottom mass are $`m_c=`$ 1.5 GeV and $`m_b=`$ 4.75 GeV, $`\mathrm{\Lambda }_3=`$ 260 MeV and $`\mathrm{\Lambda }_4=`$ 200 MeV. For the factorization scale we use $`\mu _F=\sqrt{m_b^2+Q^2/4}`$ in the case of bottom production and $`\mu _F=2\sqrt{m_c^2+Q^2/4}`$ in the charm case .

Our results for the $`Q^2`$\- and $`x`$-distributions of the AA in charm leptoproduction at several values of initial energy are presented in Fig. 4 and Fig. 5, respectively. The LO and NLO predictions are given by solid and dotted lines, correspondingly. The lines with label $`\mathrm{`}\mathrm{`}1\mathrm{"}`$ correspond to $`\rho _1=0.2`$, $`\mathrm{`}\mathrm{`}2\mathrm{"}\rho _2=0.1`$, $`\mathrm{`}\mathrm{`}3\mathrm{"}\rho _3=0.05`$ and $`\mathrm{`}\mathrm{`}4\mathrm{"}\rho _4=0.025`$, where $`\rho =4m^2/\overline{S}`$. So, in the charm case, we have: $`E_1=24`$ GeV, $`E_2=47`$ GeV, $`E_3=95`$ GeV and $`E_4=190`$ GeV, where $`E`$ is the lepton energy in the lab (nucleon rest) frame: $`E=(\overline{S}m_N^2)/(2m_N)`$. The $`y`$-distribution of the asymmetry at the COMPASS energies is given in Fig. 6. In Fig. 7 we plot the $`A(Q^2)`$, $`A(x)`$ and $`A(y)`$ distributions of the asymmetry in bottom production at the same values of $`\rho _i`$ $`=\{0.2,0.1,0.05,0.025\}`$ which correspond to the following set of initial energies: $`E_i=\{240,480,960,1920\}`$ in units of GeV.

Our calculations given in Figs. 4$``$7 represent the central result of this paper. One can see from Figs. 4 and 7 that soft-gluon corrections to $`A(Q^2)`$ are about few percent at not large $`Q^2m^2`$. At fixed values of $`x`$, the kinematical restriction $`Q^2m^2`$ leads to $`xm^2/\overline{S}`$. For this reason, radiative correction to $`A(x)`$ are small in the region of $`x\rho `$ (see Figs. 5 and 7). For comparison, we plot in Fig. 8 the so-called $`K`$-factors for $`\phi `$-integrated cross sections: $`K(Q^2)=\frac{\text{d}\sigma _{lN}^{\mathrm{NLO}}}{\text{d}Q^2}/\frac{\text{d}\sigma _{lN}^{\mathrm{LO}}}{\text{d}Q^2}`$ and $`K(x)=\frac{\text{d}\sigma _{lN}^{\mathrm{NLO}}}{\text{d}x}/\frac{\text{d}\sigma _{lN}^{\mathrm{LO}}}{\text{d}x}`$. One can see that large soft-gluon corrections to the production cross sections practically (to within few percent) do not affect the Born predictions for the $`\mathrm{cos}2\phi `$ asymmetry at $`Q^2m^2`$ and $`x\rho `$.

At fixed values of $`y`$, the allowed region of $`Q^2`$ is $`m_l^2y^2/(1y)Q^2y\overline{S}4m^2`$, where $`m_l`$ is the initial lepton mass. Since production cross sections rapidly vanish with growth of $`Q^2`$, practically whole contribution to $`A(y)`$ originates from the low-$`Q^2`$ region. For this reason, radiative corrections to $`A(y)`$ are negligible practically in the whole region of $`y`$ (see Figs. 6 and 7).

Let us now discuss the region of applicability of the adopted soft-gluon approximation. As noted in previous Section, soft radiation reproduces satisfactorily the existing exact NLO results when $`\xi `$ is not large, $`Q^2m^2`$. At large $`Q^2m^2`$, the hard $`(\stackrel{}{p}_{g,T}0)`$ contributions becomes sizeable and the quality of the NLL approximation becomes worse for both $`\phi `$-dependent and $`\phi `$-independent cross sections. Moreover, we have observed that soft-gluon approximation overestimates the exact results for $`c_T^{(1,0)}(\eta ,\xi )`$ and $`c_T^{(1,1)}(\eta ,\xi )`$ at $`\xi 1`$ and, simultaneously, underestimates the corresponding ones for $`c_A^{(1,1)}(\eta ,\xi )`$. For this reason, in the high-$`Q^2`$ region, the full NLO corrections to the $`\mathrm{cos}2\phi `$ asymmetry may be smaller than the soft-gluon ones.

As to the energy dependence, one can see from Figs. 2 and 3 that soft radiation describes very well the exact NLO results on $`\phi `$-independent photon-gluon fusion at partonic energies up to $`\eta 2`$. Since the gluon distribution function supports just the threshold region, the soft-gluon contribution dominates the photon-hadron cross sections approximately up to $`S/4m^210`$ (and, correspondingly, up to $`\rho 0.1`$ for $`\sigma _{lN}`$). Using the exact results for the $`\gamma ^{}g`$ cross sections , we have verified that the contribution originating from the region $`\eta >2`$ makes only few percent from the NLO hadron-level predictions for $`\sigma _T(S,Q^2)`$ and $`\sigma _L(S,Q^2)`$ at $`S/4m^210`$ (and $`Q^2m^2`$). Results of Ref. on soft-gluon corrections to $`F_2^{\mathrm{charm}}(x,\xi )`$ confirm our conclusion.

Presently, exact NLO calculations of the $`\phi `$-dependent cross section of heavy flavor production are not completed. However we can be sure that, at energies not so far from the production threshold, the soft radiation is the dominant perturbative mechanism in the case of $`\sigma _A(S,Q^2)`$ too. First, LO predictions for the $`\mathrm{cos}2\phi `$-dependent cross section are large and the Sudakov logarithms have universal, $`\phi `$-independent structure. For this reason, $`\sigma _A(S,Q^2)`$ has also a strong threshold enhancement. Second, our analysis of the exact scale-dependent cross section $`c_A^{(1,1)}(\eta ,\xi )`$ given in Figs. 2 and 3 confirms with a good accuracy the dominance of the soft-gluon contribution. These facts argue that hard and virtual corrections to the $`\mathrm{cos}2\phi `$-dependent cross section cannot affect significantly the soft-gluon predictions for the azimuthal asymmetry at low $`Q^2`$ in the energy region up to $`S/4m^210`$.

Let us briefly discuss the origin of perturbative stability of the $`\mathrm{cos}2\phi `$ asymmetry. Note that the mere $`\phi `$-independent structure of the Sudakov logarithms cannot explain our results since perturbative stability does not take place at the parton level. In fact, the ratios $`\frac{c_A^{(1,0)}}{c_T^{(1,0)}}(\eta ,\xi )`$ and $`\frac{c_A^{(0,0)}}{c_T^{(0,0)}}(\eta ,\xi )`$ differ essentially from each other even at $`\eta ,\xi 1`$. This is due to the fact that the physical soft-gluon corrections (24) are determined by a convolution of the Born cross sections with the Sudakov logarithms which, apart from factorized $`\delta (s_4)`$-terms, contain also non-factorizable ones (see Eq. (32)). Kinematically, sizeable values of $`\eta 1`$ allow $`s_4/m^21`$ that leads to significant non-factorizable corrections. In other words, collinear bremsstrahlung carries away a large part of initial energy. Since the $`\phi `$-dependent and $`\phi `$-independent Born level partonic cross sections have different energy behavior, the soft radiation has different impact on these quantities.

Our analysis shows that two more factors are responsible for perturbative stability of the hadron level asymmetry. First, one can see from Figs. 2 and 3 that both $`\phi `$-dependent and $`\phi `$-independent Born level cross sections take their maximum values practically at the same values of $`\eta `$. Second, at fixed target energies, the gluon distribution function supports the contribution of the threshold region. In other words, sufficiently soft gluon distribution function makes the collinear gluon radiation effectively soft at the hadron level. In detail, the role of the gluon distribution function in perturbative stability of the azimuthal asymmetry in heavy quark photoproduction is discussed in Ref..

Another remarkable property of the azimuthal asymmetry closely related to fast perturbative convergence is its parametric stability.<sup>§</sup><sup>§</sup>§Of course, parametric stability of the fixed order results does not imply a fast convergence of the corresponding series. However, a fast convergent series must be parametrically stable. In particular, it must be $`\mu _R`$\- and $`\mu _F`$-independent. Our analysis shows that the pQCD predictions for the $`\mathrm{cos}2\phi `$ asymmetry are less sensitive to standard uncertainties in the QCD input parameters than the corresponding ones for the production cross sections. For instance, changes of $`\mu _F`$ in the range $`m_c<\mu _F<2\sqrt{m_c^2+Q^2/4}`$ affect the quantity $`A(Q^2)`$ in charm production by less than $`7\%`$ at $`\rho =0.025`$ and $`\xi 4`$. For the $`\phi `$-integrated cross section, such changes of $`\mu _F`$ lead to $`30\%`$ variations in the same kinematics. We have also verified that all the NLO CTEQ5 versions of gluon density as well as the LO parametrization  lead to asymmetry predictions which coincide with each other with an accuracy of better than $`1.5\%`$.

Parametric stability of the azimuthal asymmetry leads to the scaling: with a good accuracy the quantity $`A(\rho ,x,\xi )`$ in Eq. (8) is a function of three variables, so that

$$A^{\mathrm{Charm}}(\rho ,x,\xi )A^{\mathrm{Bottom}}(\rho ,x,\xi )$$
(39)
at the same values of $`\rho `$, $`x`$ and $`\xi `$. To illustrate this property, in Fig. 9 we plot the asymmetry parameter $`A(Q^2)`$ defined by Eqs. (38) as a function of $`\xi `$ at $`\rho =10^3`$ and $`10^4`$. Since the soft-gluon approximation is inapplicable to heavy flavor production at high energies, we give only the LO predictions for $`A(\xi )`$. The LO results for the charm and bottom cases are plotted by solid and dash-dotted lines, respectively. One can see that both curves coincide with each other with an accuracy of better than $`1\%`$.

It is also seen from Fig. 9 that, at the HERA energy ($`\rho 10^3`$ and $`10^4`$ for bottom and charm quark, respectively), the pQCD predictions for the AA in heavy quark leptoproduction are large and can be tested experimentally.

### B Nonperturbative corrections

Let us discuss how the pQCD predictions for azimuthal asymmetry are affected by nonperturbative contributions due to the intrinsic transverse motion of the gluon in the target. Because of the relatively low $`c`$-quark mass, these contributions are especially important in the description of the cross sections for charmed particle production .

To introduce $`k_T`$ degrees of freedom, $`\stackrel{}{k}_gz\stackrel{}{p}+\stackrel{}{k}_T`$, one extends the integral over the parton distribution function in Eq. (18) to $`k_T`$-space,

$$\text{d}zg(z,\mu _F)\text{d}z\text{d}^2k_Tf\left(\stackrel{}{k}_T\right)g(z,\mu _F).$$
(40)
The transverse momentum distribution, $`f\left(\stackrel{}{k}_T\right)`$, is usually taken to be a Gaussian:

$$f\left(\stackrel{}{k}_T\right)=\frac{\mathrm{e}^{k_T^2/k_T^2}}{\pi k_T^2}.$$
(41)
In practice, an analytic treatment of $`k_T`$ effects is usually used. According to , the $`k_T`$-smeared differential cross section of the process (1) is a 2-dimensional convolution:

$$\frac{\text{d}^4\sigma _{lN}^{\mathrm{kick}}}{\text{d}x\text{d}Q^2\text{d}p_{QT}\text{d}\phi }\left(\stackrel{}{p}_{QT}\right)=\text{d}^2k_T\frac{\mathrm{e}^{k_T^2/k_T^2}}{\pi k_T^2}\frac{\text{d}^4\sigma _{lN}}{\text{d}x\text{d}Q^2\text{d}p_{QT}\text{d}\phi }\left(\stackrel{}{p}_{QT}\frac{1}{2}\stackrel{}{k}_T\right).$$
(42)
The factor $`\frac{1}{2}`$ in front of $`\stackrel{}{k}_T`$ in the r.h.s. of Eq. (42) reflects the fact that the heavy quark carries away about one half of the initial energy in the reaction (1).

Values of the $`k_T`$-kick corrections to LO predictions for the $`\mathrm{cos}2\phi `$ asymmetry in the charm production are shown in Figs. 4$``$6 and 9 by dashed curves. Calculating the $`k_T`$-kick effects we use $`k_T^2=0.5`$ GeV$`^\text{2}`$. At fixed target energies, $`k_T`$-smearing for $`A(Q^2)`$ and $`A(x)`$ is about $`20`$-$`25\%`$ in the region of low $`Q^2`$ and $`x`$, respectively, and decreases at large $`Q^2`$ and $`x`$. In the HERA range, expected values of the $`k_T`$-corrections are systematically smaller (see Fig. 9).

Analogous calculations for the case of bottom production are presented in Figs. 7 and 9. It is seen that $`k_T`$-kick corrections to the $`b`$-quark AA are practically negligible in the whole region of $`Q^2`$ and $`x`$.

## IV Conclusion

In this paper we have investigated the impact of soft-gluon radiation on the $`\mathrm{cos}2\phi `$ asymmetry in heavy flavor leptoproduction. The NLL approximation provides a good description of the available exact NLO results for $`Q^2m^2`$ at energies of the fixed target experiments. Our calculations show that the azimuthal asymmetry is practically insensitive to soft-gluon corrections in this kinematics. We conclude that, unlike the $`\phi `$-integrated cross sections, the $`\mathrm{cos}2\phi `$ asymmetry in heavy quark leptoproduction is an observable quantitatively well defined in pQCD: it is stable both parametrically and perturbatively, and insensitive (in the case of bottom production) to nonperturbative contributions. This asymmetry is of leading twist and predicted to be about $`15\%`$ at energies sufficiently above the production threshold for both charm and bottom quark. Measurements of the $`\mathrm{cos}2\phi `$ asymmetry in bottom production at HERA would provide an ideal test of pQCD.

Data on the charm azimuthal distributions from the COMPASS and HERMES experiments would make it possible to clarify the role of subleading twist contributions. Our analysis shows that, in the low-$`x`$ region, the AA is sensitive to the gluon transverse motion in the target. At high $`x`$, the intrinsic charm contribution  to the structure functions may be significant . In detail, the possibility of measuring the intrinsic charm content of the proton using the $`\mathrm{cos}2\phi `$ asymmetry will be considered in a forthcoming publication.

###### Acknowledgements.

The author would like to thank S.J. Brodsky, A.B. Kaidalov, A.Kotzinian and A.G. Oganesian for useful discussions. I am grateful to High Energy Section of ICTP for hospitality while this work has been completed.


format: text

text:
# Approximation Thresholds for Combinatorial Optimization Problems

## 1. Introduction

Given an instance $`I`$ of a combinatorial optimization problem $`\mathrm{\Pi }`$, one wishes to find a feasible solution of optimal value to $`I`$. For example, in the maximum clique problem, the input instance is a graph $`G(V,E)`$, a feasible solution is a clique (a set of vertices $`SV`$ such that for all vertices $`u,vS`$, edge $`(u,v)E`$), the value of the clique is its size, and the objective is to find a feasible solution of maximum value. By introducing an extra parameter $`k`$ to a combinatorial maximization problem, one can formulate a decision problem of the form “does $`I`$ have a feasible solution of value at least $`k`$?” (or “at most $`k`$”, for a minimization problem). In this work we consider only combinatorial optimization problems whose decision version is NP-complete. Hence solving these problems is NP-hard, informally meaning that there is no polynomial time algorithm that is guaranteed to output the optimal solution for every input instance (unless P=NP). (For information on the theory of NP-completeness, see , or essentially any book on computational complexity.)

One way of efficiently coping with NP-hard combinatorial optimization problems is by using polynomial time approximation algorithms. These algorithms produce a feasible solution that is not necessarily optimal. The quality of the approximation algorithm is measured by its approximation ratio: the worst case ratio between the value of the solution found by the algorithm and the value of the optimal solution. Numerous approximation algorithms have been designed for various combinatorial optimization problems, applying a diverse set of algorithmic techniques (such as greedy algorithms, dynamic programming, linear programming relaxations, applications of the probabilistic method, and more). See for example .

When one attempts to design an approximation algorithm for a specific problem, it is natural to ask whether there are limits to the best approximation ratio achievable. The theory of NP-completeness is useful in this context, allowing one to prove NP-hardness of approximation results. Namely, for some value of $`\rho `$, achieving an approximation ratio better than $`\rho `$ is NP-hard. Of particular interest are threshold results. A combinatorial optimization problem $`\mathrm{\Pi }`$ is said to have an approximation threshold at $`t`$ if there is a polynomial time algorithm that approximates the optimal value within a ratio of $`t`$, and a hardness of approximation result that shows that achieving approximation ratios better than $`t`$ is NP-hard.

The notion of an approximation threshold is often not fully appreciated, so let us discuss it in more detail. First, let us make the point that thresholds of approximation (in the sense above) need not exist at all. One may well imagine that for certain problems there is a gradual change in the complexity of achieving various approximation ratios: solving the problem exactly is NP-hard, approximating it within a ratio of $`\rho `$ can be done in polynomial time, and achieving approximation ratios between $`\rho `$ and 1 is neither in P nor NP-hard, but rather of some intermediate complexity. The existence of an approximation threshold says that nearly all approximation ratios (except for ratios that differ from the threshold only in low order terms) are either NP-hard to achieve, or in P. This is analogous to the well known empirical observation that “most” combinatorial optimization problems that people study turn out to be either in P or NP-hard, with very few exceptions. But note that in the context of approximation ratios, the notion of “most” is well defined. The second point that we wish to make is that an NP-hardness of approximation result is really a polynomial time algorithm. This algorithm reduces instances of 3SAT (or of some other NP-complete language) to instances of $`\mathrm{\Pi }`$, and an approximation within a ratio better than $`\rho `$ for $`\mathrm{\Pi }`$ can be used in order to solve 3SAT. Hence the threshold is a meeting point between two polynomial time algorithms: the reduction and the approximation algorithm. Here is another way of looking at it. For problem $`\pi `$, we can say that approximation ratio $`\rho _1`$ reduces to $`\rho _2`$ if there is a polynomial time algorithm that can approximate instances of $`\pi `$ within ratio $`\rho _1`$ by invoking a subroutine that approximates instances of $`\pi `$ within ratio $`\rho _2`$. (Each call to the subroutine counts as one time unit.) Two approximation ratios are equivalent if they are mutually reducible to each other. The existence of a threshold of approximation for problem $`\mathrm{\Pi }`$ says that essentially all approximation ratios for it fall into two equivalent classes (those above the threshold and those below the threshold). A-priori, it is not clear why the number of equivalent classes should be two.

Hardness of approximation results are often (though not always) proved using “PCP techniques”. For more details on these techniques, see for example the survey of Arora . We remark that in all cases the hardness results apply also to the problem of estimating the optimal value for the respective problem. An estimation algorithm is a polynomial time algorithm that outputs an upper bound and a lower bound on the value of the optimal solution (without necessarily producing a solution whose value falls within this range). The estimation ratio of the algorithm is the worst case ratio between the upper bound and the lower bound. An approximation algorithm is stronger than an estimation algorithm in a sense that it supports its estimate by exhibiting a feasible solution of the same value. Potentially, designing estimation algorithms is easier than designing approximation algorithms. See for example the remark in Section 2.3.

For many combinatorial optimization problems, approximation thresholds are known. In Section 2. we survey some of these problems. For each such problem, we sketch an efficient algorithm for estimating the optimal value of the solution. The estimation ratios of these algorithms (the analog of approximation ratios) meet the thresholds of approximation for the respective problems (up to low order additive terms), and hence are best possible (unless P=NP). One would expect these estimation algorithms to be the “state of the art” in algorithmic design. However, as we shall see, for every problem above there is a core version for which the best possible estimation algorithm is elementary: it bases its estimate only on easily computable properties of the input instance, such as the number of clauses in a formula, without searching for a solution. The core versions that we present often have the property that Hastad  characterizes as “non-approximable beyond the random assignment threshold”. In these cases (e.g., max 3SAT), the estimates on the value of the optimal solution are derived from analysing the expected value of a random solution.

In contrast, for problems that are known to have a polynomial time approximation scheme (namely, that can be approximated within ratios arbitrarily close to 1), their approximation algorithms do perform an extensive search for a good solution, often using dynamic programming.

The empirical findings are discussed in Section 3.

## 2. A survey of some threshold results

In this section we survey some of the known threshold of approximation results. For each problem we show simple lower bounds and upper bounds on the value of the optimal solution. Obtaining any better bounds is NP-hard. (There are certain exceptions to this. For set cover and domatic number the matching hardness result are under the assumption that NP does not have slightly super-polynomial time algorithms. The hardness results for clique and chromatic number assume that NP does not have randomized algorithms that run in expected polynomial time.).

Conventions. For a graph, $`n`$ denotes the number of vertices, $`m`$ denotes the number of edges, $`\delta `$ denotes the minimum degree, $`\mathrm{\Delta }`$ denotes the maximum degree. For a formula, $`n`$ denotes the number of variables, $`m`$ denotes the number of clauses. For each problem we define its inputs, feasible solutions, value of solutions, and objective. We present the known approximation ratios and hardness results, which are given up to low order additive terms. We then present a core version of the problem. For the core version, we present an upper bound and a lower bound on the value of the optimal solution. Improving over these bounds (in the worst case) is NP-hard. In most cases, we give hints to the proofs of our claims. We also cite references were full proofs can be found.

### 2.1. Max coverage

Input. A collection $`S_1,\mathrm{},S_m`$ of sets with $`_{i=1}^mS_i=\{1,\mathrm{},n\}`$. A parameter $`k`$.

Feasible solution. A collection $`I`$ of $`k`$ indices.

Value. Number of covered items, namely $`|_{iI}S_i|`$.

Objective. Maximize.

Algorithm. Greedy. Iteratively add to $`I`$ the set containing the maximum number of yet uncovered items, breaking ties arbitrarily.

Approximation ratio. $`11/e`$.

Hardness. $`11/e`$.

Core. $`d`$-regular, $`r`$-uniform. Every set is of cardinality $`d`$. Every item is in $`r`$ sets. $`k=n/d`$.

Upper bound. $`n`$.

Lower bound. $`(11/e)n`$. (Pick $`k=n/d=m/r`$ sets at random.)

Hardness of core. $`(11/e+ϵ)`$ for every $`ϵ>0`$, when $`d,r`$ are large enough.

Remarks. For the special case in which each item belongs to exactly 4 sets, and $`k=m/2`$, there always is a choice of sets covering at least $`15n/16`$ items. For very $`ϵ>0`$, a $`(1+ϵ)15/16`$ approximation ratio is NP-hard .

### 2.2. Min set cover

Input. A collection $`S_1,\mathrm{},S_m`$ of sets with $`_{i=1}^mS_i=\{1,\mathrm{},n\}`$.

Feasible solution. A set of indices $`I`$ such that $`_{iI}S_i=\{1,\mathrm{},n\}`$.

Value. Number of sets used in the cover, namely, $`|I|`$.

Objective. Minimize.

Algorithm. Greedy. Iteratively add to $`I`$ the set containing the maximum number of yet uncovered items, breaking ties arbitrarily.

Approximation ratio. $`\mathrm{ln}n`$.

Hardness. $`\mathrm{ln}n`$.

Core. $`d`$-regular, $`r`$-uniform. Every set is of cardinality $`d`$. Every item is in $`r`$ sets.

Lower bound. $`n/d`$.

Upper bound. $`(n/d)\mathrm{ln}n`$, up to low order terms. (Include every set in $`I`$ independently with probability $`\frac{\mathrm{ln}n}{r}`$.)

Hardness of core. $`\mathrm{ln}n`$.

Remarks. The hardness results for set cover in  assume that NP does not have deterministic algorithms that run in slightly super-polynomial time (namely, time $`n^{O(\mathrm{log}\mathrm{log}n)}`$).

### 2.3. Domatic number

Input. A graph.

Feasible solution. A domatic partition of the graph. That is, a partition of the vertices of the graph into disjoint sets, where each set is a dominating set. (A dominating set is a set $`S`$ of vertices that is adjacent to every vertex not in $`S`$.)

Value. Number of dominating sets in the partition.

Objective. Maximize.

Algorithm. Let $`\delta `$ be the minimum degree in the graph. Partition the vertices into $`(1ϵ)(\delta +1)/\mathrm{ln}n`$ sets at random, where $`ϵ`$ is arbitrarily small when $`\delta /\mathrm{ln}n`$ is large enough. Almost all these sets will be dominating. The sets that are not dominating can be unified with the first of the dominating sets to give a domatic partition. The algorithm can be derandomized to give (a somewhat unnatural) greedy algorithm.

Approximation ratio. $`1/\mathrm{ln}n`$.

Hardness. $`1/\mathrm{ln}n`$.

Core. $`\delta /\mathrm{ln}n`$ is large enough.

Upper bound. $`\delta +1`$.

Lower bound. $`(1ϵ)(\delta +1)/\mathrm{ln}n`$.

Hardness of core. $`(1+ϵ)/\mathrm{ln}n`$, for every $`ϵ>0`$.

Remarks. The hardness results assume that NP does not have deterministic algorithms that run in time $`n^{O(\mathrm{log}\mathrm{log}n)}`$. The lower bound can be refined to $`(1ϵ)(\delta +1)/\mathrm{ln}\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ is the maximum degree in the graph. This is shown using a nonconstructive argument (based on a two phase application of the local lemma of Lovasz). It is not known how to find such a domatic partition in polynomial time.

### 2.4. $`k`$-center

Input. A metric on a set of $`n`$ points and a parameter $`k<n`$.

Feasible solution. A set $`S`$ of $`k`$ of the points.

Value. Distance between $`S`$ and point furthest away from $`S`$.

Objective. Minimize.

Algorithm. Greedy. Starting with the empty set, iteratively add into $`S`$ the point furthest away from $`S`$, resolving ties arbitrarily.

Approximation ratio. 2.

Hardness. 2.

Core. The $`n`$ points are vertices of a graph with minimum degree $`\delta `$, equipped with a shortest distance metric. $`kn/(\delta +1)`$.

Lower bound. 1. (Because $`k<n`$.)

Upper bound. 2. (As long as there is a vertex of distance more than 2 from the current set $`S`$, every step of the greedy algorithm covers at least $`\delta +1`$ vertices.)

Hardness of core. 2. (Because it is NP-hard to tell if the underlying graph has a dominating set of size $`k`$.)

### 2.5. Min sum set cover

Input. A collection $`S_1,\mathrm{},S_m`$ of sets with $`_{i=1}^mS_i=\{1,\mathrm{},n\}`$.

Feasible solution. A linear ordering. That is, a one to one mapping $`f`$ from the collection of sets to $`\{1,\mathrm{},m\}`$.

Value. The average time by which an item is covered. Namely, $`\frac{1}{n}_i\mathrm{min}_{\{j|iS_j\}}f(j)`$.

Objective. Minimize.

Algorithm. Greedy. Iteratively pick the set containing the maximum number of yet uncovered items, breaking ties arbitrarily.

Approximation ratio. 4.

Hardness. 4.

Core. $`d`$-regular, $`r`$-uniform. Every set is of cardinality $`d`$. Every item is in $`r`$ sets.

Lower bound. $`n/2d=m/2r`$. (Because every set covers $`d`$ items.)

Upper bound. $`m/(r+1)`$. (By ordering the sets at random.)

Hardness of core. $`2ϵ`$, for every $`ϵ>0`$, when $`r`$ is large enough.

Remarks. Note that the core has a lower threshold of approximation than the general case.

### 2.6. Min bandwidth

Input. A graph.

Feasible solution. A linear arrangement. That is, a one to one mapping $`f`$ from the set of vertices to $`\{1,\mathrm{},n\}`$.

Value. Longest stretch of an edge. Namely, the maximum over all edges $`(i,j)`$ of $`|f(i)f(j)|`$.

Objective. Minimize.

Core. The input graph $`G`$ is a unit length circular arc graph. Namely, the vertices represent arcs of equal length on a circle. Two vertices are connected by an edge if their respective arcs overlap. Let $`\omega (G)`$ denote the size of the maximum clique in $`G`$. (In circular arc graphs, a clique corresponds to a set of mutually intersecting arcs, and $`\omega (G)`$ can be computed easily.)

Lower bound. $`\omega (G)1`$.

Upper bound. $`2\omega (G)2`$.

Hardness of core. 2.

Remarks. No threshold of approximation is known for the problem on general graphs.

### 2.7. Max 3XOR

Input. A logical formula with $`n`$ variables and $`m`$ clauses. Every clause is the exclusive or of three distinct literals.

Feasible solution. An assignment to the $`n`$ variables.

Value. Number of clauses that are satisfied.

Objective. Maximize.

Algorithm. Gaussian elimination can be used in order to test if the formula is satisfiable. If the formula is not satisfiable, pick a random assignment to the variables.

Approximation ratio. $`1/2`$. (In expectation a random assignments satisfies $`m/2`$ clauses, whereas no assignment satisfies more than $`m`$ clauses.)

Hardness. $`1/2`$.

Core. Same as above.

Upper bound. $`m`$.

Lower bound. $`m/2`$.

Hardness of core. $`1/2`$.

### 2.8. Max 3AND

Input. A logical formula with $`n`$ variables and $`m`$ clauses. Every clause is the and of three literals.

Feasible solution. An assignment to the $`n`$ variables.

Value. Number of clauses that are satisfied.

Objective. Maximize.

Algorithm. Based on semidefinite programming.

Approximation ratio. $`1/2`$.

Hardness. $`1/2`$.

Core. Pairwise independent version. Every variable appears at most once in each clause. For a pair of variables $`x_i`$ and $`x_j`$, let $`n_{ij}(0,0)`$ ($`n_{ij}(0,1)`$, $`n_{ij}(1,0)`$, $`n_{ij}(1,1)`$, respectively) denote the number of clauses in which they both appear negated ($`i`$ negated and $`j`$ positive, $`i`$ positive and $`j`$ negated, both positive, respectively). Then for every $`i`$ and $`j`$, $`n_{ij}(0,0)=n_{ij}(0,1)=n_{ij}(1,0)=n_{ij}(1,1)`$.

Upper bound. $`m/4`$. (Consider the pairs of literals in the satisfied clauses. There must be at least three times as many pairs in unsatisfied clauses.)

Lower bound. $`m/8`$. (The expected number of clauses satisfied by a random assignment.)

Hardness of core. $`1/2`$.

Remarks. This somewhat complicated core comes up as an afterthought, by analysing the structure of instances that result from the reduction from max 3XOR.

### 2.9. Max 3SAT

Input. A logical formula with $`n`$ variables and $`m`$ clauses. Every clause is the or of three literals.

Feasible solution. An assignment to the $`n`$ variables.

Value. Number of clauses that are satisfied.

Objective. Maximize.

Algorithm. Based on semidefinite programming.

Approximation ratio. $`7/8`$.

Hardness. $`7/8`$.

Core. In every clause, the three literals are distinct.

Upper bound. $`m`$.

Lower bound. $`7m/8`$. (The expectation for a random assignment.)

Hardness of core. $`7/8`$.

Remarks. The approximation ratio of the algorithm of  is determined using computer assisted analysis.

### 2.10. Inapproximable problems

For some problems the best approximation ratios known are of the form $`n^{1ϵ}`$ for every $`ϵ>0`$ (where the range of the objective function is between 1 and $`n`$). This is often interpreted as saying that approximation algorithms are almost helpless with respect to these problems. Among these problems we mention finding the smallest maximal independent set , max clique , and chromatic number . The hardness results in  are proved under the assumption that NP does not have expected polynomial time randomized algorithms.

### 2.11. Thresholds within multiplicative factors

For some problems, the best approximation ratio is of the form $`O(n^c)`$, for some $`0<c<1`$, and there is a hardness of approximation result of the form $`\mathrm{\Omega }(n^c^{})`$, where $`c^{}`$ can be chosen arbitrarily close to $`c`$. We may view these as thresholds of approximation up to a low order multiplicative factor. An interesting example of this sort is the max disjoint paths problem .

Input. A directed graph and a set $`S`$ of pairs of terminals $`\{(s_1,t_1),\mathrm{}(s_k,t_k)\}`$.

Feasible solution. A collection of edge disjoint paths, each connecting $`s_i`$ to $`t_i`$ for some $`i`$.

Value. Number of pairs of terminals from $`S`$ that are connected by some path in the solution.

Objective. Maximize.

Algorithm. Greedy. Iteratively add the shortest path that connects some yet unconnected pair from $`S`$.

Approximation ratio. $`m^{1/2}`$, where $`m`$ is the number of arcs in the graph.

Hardness. $`m^{1/2+ϵ}`$, for every $`ϵ>0`$.

Core. There is a path from $`s_1`$ to $`t_1`$.

Upper bound. $`k`$. At most all pairs can be connected simultaneously by disjoint paths.

Lower bound. 1. There is a path from $`s_1`$ to $`t_1`$.

Hardness of core. $`1/k`$. Here $`k`$ can be chosen as $`m^{1/2ϵ}`$ for arbitrarily small $`ϵ`$.

## 3. Discussion

We summarized the pattern presented in Section 2. Given an NP-hard combinatorial optimization problem that has an approximation threshold, we identify a core version of the problem. For the core version we identify certain key parameters. Thereafter, an upper bound and a lower bound on the value of the optimal solution is expressed by a formula involving only these parameters. Even an algorithm that examines the input for polynomial time cannot output better lower bounds or upper bounds on the value of the optimal solution (in the worst case, and up to low order terms), unless P=NP.

Note that if we do not restrict what qualifies as being a key parameter, then the pattern above can be enforced on essentially any problem with an approximation threshold. We can simply take the key parameter to be the output of an approximation algorithm for the same problem. Likewise, if we do not restrict what qualifies as a core version, we can simply take the core version to be all those input instances on which a certain approximation algorithm outputs a certain value.

Hence we would like some restrictions on what may qualify as a key parameter, or as a core version. One option is to enforce a computational complexity restriction. Namely, the core should be such that deciding whether an input instance belongs to the core is computationally easy, for example, computable in logarithmic space. (Note that in this respect, the core that we defined for the max disjoint paths problem may be problematic, because testing whether there is a directed path from $`s_1`$ to $`t_1`$ is complete for nondeterministic logarithmic space.) Likewise, computing the key parameters should be easy. Another option is to enforce structural restrictions. For example, membership in the core should be invariant over renaming of variables. Ideally, the notions of “core” and “key parameter” should be defined well enough so that we should be able to say that certain classes of inputs do not qualify as being a core (e.g., because the class is not closed under certain operations), and that certain parameters are not legitimate parameters to be used by an estimation algorithm. Also, the definitions should allow in principle the possibility that certain problems have approximation thresholds without having a core.

Hand in hand with suggesting formal definitions to the notion of a core, it would be useful to collect more data. Namely, to find more approximation threshold results, and to identify plausible core versions to these problems. As the case of min bandwidth shows, one may find core versions even for problems for which an approximation threshold is not known.

In this manuscript we mainly addressed the issue of collecting data regarding approximation thresholds, and describing this data using the notion of a core. The issue of uncovering principles that explain why the patterns discussed in this manuscript emerge is left to the reader.


text:
# XMM observation of the dynamically young galaxy cluster CL 0939+4713

## 1 Introduction

The galaxy cluster CL 0939+4713 (Abell 851) at a redshift $`z=0.41`$ is an exceptionally rich cluster. It was extensively observed in morphological and spectroscopic studies of the galaxy population. These studies are based on ground-based optical and near-infrared (Dressler & Gunn 1992; Belloni et al. 1995; Fukugita et al. 1995; Stanford et al. 1995; Belloni & Röser 1996; Dressler et al. 1999; Poggianti et al. 1999; Iye et al. 2000; Martin et al. 2000) and HST observations (Dressler et al. 1993, 1994; Smail et al. 1999; Ziegler et al. 1999; Ferreras & Silk 2000) as well as observations in UV (Buson et al. 2000) and sub-millimeter (Cowie et al. 2002) wavelengths. A high number of star-burst and post-starburst galaxies was found in this distant cluster. To study the evolution of galaxies this number was compared with the corresponding number in nearby clusters. In the surroundings of the clusters a network of galaxy filaments and subclumps was found (Kodama et al. 2001).
The cluster was also the target of several gravitational lensing analyses (Seitz et al. 1996; Geiger & Schneider 1999; Iye et al. 2000) aiming at the determination of the cluster mass distribution. The comparison between the masses determined by the weak lensing analysis by Seitz et al. (1996) and the X-ray analysis by Schindler et al. (1998) showed a discrepancy of a factor of 2-3.
Previous X-ray observations by ROSAT and ASCA showed that the cluster has a very irregular structure indicating that the cluster is not relaxed but dynamically young (Schindler & Wambsganss 1996; Schindler et al. 1998): this result is in good agreement with the conclusion from larger scales by Kodama et al. (2001). The formation and evolution of clusters depends sensitively on cosmological parameters like the mean matter density in the universe $`\mathrm{\Omega }_m`$ (Thomas et al. 1998; Jenkins et al. 1998; Beisbart et al. 2001). Therefore it is important to determine the dynamical state of clusters at different redshifts, i.e. at different evolutionary states. The X-ray morphology alone is not the best indicator of the dynamical state, but it should be complemented with all other information available, e.g. the temperature map or the galaxy distribution. The combination of all findings gives a detailed picture of the state of a cluster. We do this analysis here for the moderately distant cluster CL 0939+4713, which can be used as a basis for comparisons at lower and higher redshifts.
The capability of XMM to perform spatially resolved spectroscopy can be used also to determine the distribution of the metal abundances. Not only the overall value of metallicity but also its spatial distribution gives important indications on the metal enrichment processes and the cluster formation process. Early enrichment would point preferentially towards galactic winds driven by supernovae (De Young 1978). For this process one would expect a priori no metallicity variations, because there was enough time to distribute the metals uniformly throughout the cluster. However, numerical simulations of this mechanism gave very discordant results. Aguirre et al. (2001) and Metzler & Evrard (1994, 1997) concluded that winds are very efficient in enriching the ICM, but the latter authors found a very steep radial metallicity gradient, which is not observed in any cluster. Murakami & Babul (1999) and Gnedin (1998) on the other hand found that galactic winds play only a minor role.
At a later stage ram pressure stripping (Gunn & Gott 1972) is expected to become increasingly important. If considerable amounts of metals are stripped recently a spatial variation – not necessarily radial gradients – should be observable. To test these effects it is particularly important to analyse clusters at different redshifts to see how the efficiency of the different enrichment mechanisms changes with time.

Throughout this paper we use $`H_0=50\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and $`q_0=0.5`$; all errors are $`90\%`$ confidence levels.

## 2 Observation and data analysis

### 2.1 Observation and data processing

The XMM-Newton observation of the galaxy cluster CL 0939+4713 (ID 0106460101) was performed with the European Photon Imaging Camera (EPIC) on November 6th, 2000 during the satellite’s $`167^{\mathrm{th}}`$ revolution. It consists of one pointing obtained operating the EPIC cameras in the standard Full Frame Mode. For all the three EPIC cameras (MOS1, MOS2, and pn) the thin filter was used for a total exposure time of $`51\mathrm{ks}`$.
We used the SASv5.0.1 processing tasks emchain and epchain to generate calibrated event files from raw data.
The XMM-Newton observations are strongly affected by high internal and external background levels. Efficient ways of subtracting these background components were extensively discussed by Pratt et al. (2002) and Lumb et al. (2002).

### 2.2 Flare rejection

Low energy ($`<1\mathrm{MeV}`$) protons from solar flares are stopped in the first microns of silicon of the EPIC CCDs. They give a signal identical to X-rays, which makes it impossible to recognize them via a pattern selection. We obtained a cleaner set of data by visually inspecting the light curves and applying an intensity filter. In the energy range $`1012\mathrm{keV}`$ for MOS1 & MOS2 cameras ($`1214\mathrm{keV}`$ for pn) we excluded all the intervals of exposure time having a count rate higher than a certain threshold value.

The three light curves, one for each instrument, were analyzed separately in order to find the best threshold values (see Fig. 1). Different threshold values lead to different cuts in the exposure time (see Fig. 2); the selected values were therefore a compromise between maximizing the exposure time and minimizing the presence of events caused by flares.
The chosen threshold values are 15 counts $`100\mathrm{s}^1`$ for MOS1 and MOS2 and 25 counts $`100\mathrm{s}^1`$ for pn; the threshold for the pn camera is higher because of the larger effective area of the instrument compared to MOS. The resulting effective exposure times for the three cameras are

MOS1: $`45800\mathrm{s}`$

MOS2: $`46300\mathrm{s}`$

pn: $`41800\mathrm{s}`$

### 2.3 Vignetting correction

To correct for vignetting the method described in detail by Arnaud et al. (2001a) has been used. A weighting, given by the ratio of the effective area at the centre of the detector -$`A_{0,0}(E_i)`$\- to the effective area at the event position -$`A_{x_j,y_j}(E_i)`$-, is assigned to each photon extracted in each spectrum. The vignetting corrected photon counts -$`C(E_i)`$\- in each energy channel $`E_j`$ ($`E_i\mathrm{\Delta }E_i/2<E_j<E_i+\mathrm{\Delta }E_i/2`$) are

$$C(E_i)=\underset{j}{}w_j(x_j,y_j)$$

with the weights $`w_j`$ being, as described above:

$$w_j=\frac{A_{0,0}(E_i)}{A_{x_j,y_j}(E_i)}$$

which therefore all have values greater than unity.
The errors in the spectra are recomputed taking into account the vignetting correction; the variances in each energy channel are therefore given by

$$\sigma ^2(C(E_i))=\underset{j}{}w_j^2.$$

The above method is easily generalized for imaging, as described by Majerowicz et al. (2002).

### 2.4 Background subtraction

In order to subtract both the cosmic and the non cosmic X-ray background (Arnaud et al. 2002), the blank field observations, collected by D. Lumb were used. 111Event files can be downloaded via anonymous ftp at xmm.vilspa.esa.es in the directory /pub/ccf/constituents/extras/background. These are homogeneous collections of high Galactic latitude observations where no bright source is present, for a total exposure time of about $`400\mathrm{ksec}`$. These observations were all taken in a period of very low flare occurrence and almost exclusively with the use of a thin filter; bad pixels and sources were also removed.
We applied to the background event files the same intensity filter applied to our observations; this yields final exposure times for the blank fields of MOS1: $`373425\mathrm{s}`$, MOS2: $`358904\mathrm{s}`$, and pn: $`332646\mathrm{s}`$. Since the cosmic ray induced background might slightly change with time, the count rate in the $`10.012.0\mathrm{keV}`$ ($`12.014.0\mathrm{keV}`$) energy band for MOS1 & MOS2 (pn) cameras in the blank field was compared with the count rate in our observation (after flare rejection). The obtained normalization factors ($`0.995`$, $`0.992`$ and $`1.038`$ for MOS1, MOS2 and pn, respectively) are then used to renormalize the blank field data.
We then looked for a possible residual soft background due to the soft component ($`\mathrm{E}<1.5\mathrm{keV}`$) of the cosmic X-ray diffuse emission, which varies with the position in the sky. At the position in the sky of CL 0939+4713, no excess of soft cosmic background was observed in previous observations (Snowden et al. 1997).
Nevertheless, we performed a further check. We chose a region in our pointing where no emission from the cluster, or from other point sources, was observed. For all the three detectors, spectra were extracted from this region, together with background spectra from the blank fields by D. Lumb. No evident excess of flux, associated with the soft X-ray background, is observed in any of the three detectors. We then subtracted the blank field spectra from the ones extracted in our observation; the resulting spectra yield no significant signal and are dominated by noise only.
No additional correction for excess of soft background is therefore needed for CL 0939+4713.

## 3 Morphological analysis

Figure 4 shows an image of CL 0939+4713 in the $`0.32.0\mathrm{keV}`$ band. The image is an adaptively smoothed and exposure corrected mosaic of photon events from the three detectors (MOS1, MOS2 and pn).
The cluster emission has an elliptical shape, centred at RA$`=9^\mathrm{h}43^\mathrm{m}01.00^\mathrm{s}`$, DEC$`=+46\text{}59{}_{}{}^{}37.50^{\prime \prime }`$ (J2000). The axis ratio of major to minor axis is $`1.3`$. The position angle of the ellipse is P.A.$`63\text{}`$ (N over E).

### 3.1 Substructure

Figure 4 clearly indicates that the core of CL 0939+4713 does not consist of a single spherically symmetric structure, but is composed of two subclusters: M1 and M2.
In Fig. 3 (left panel) a surface brightness projection of the cluster is shown, which is extracted from the region shown in the right panel. The two extended substructures forming the core are of similar size and of similar brightness, the subcluster in the East –M2– being slightly brighter than the subcluster in the West –M1–). They both have an irregular shape, showing signs of further internal substructure. Their shapes are roughly elliptical with minor axis of $`120\mathrm{kpc}`$ ($`19^{\prime \prime }`$) and major axis of $`180\mathrm{kpc}`$ ($`28^{\prime \prime }`$).222$`1{}_{}{}^{}=389\mathrm{kpc}`$ at the cluster distance of $`\mathrm{z}=0.406`$. Their peaks are $`280\mathrm{kpc}`$ apart and their coordinates (J2000) are:

$`\mathrm{P}_{\mathrm{M1}}`$: RA$`=9^\mathrm{h}42^\mathrm{m}58.9^\mathrm{s}`$, DEC$`=+46\text{}59{}_{}{}^{}36.0^{\prime \prime }`$

$`\mathrm{P}_{\mathrm{M2}}`$: RA$`=9^\mathrm{h}43^\mathrm{m}03.3^\mathrm{s}`$, DEC$`=+46\text{}59{}_{}{}^{}26.5^{\prime \prime }`$

They are slightly offsetted ($`\mathrm{P}_{\mathrm{M1}}`$ $`44^{\prime \prime }`$ north and $`\mathrm{P}_{\mathrm{M2}}`$ $`21^{\prime \prime }`$ south) compared to the position observed by Schindler et al. (1998), and reach a maximum surface brightness of $`5.0\times 10^2\mathrm{counts}\mathrm{s}^1\mathrm{arcmin}^2`$ and $`5.8\times 10^2\mathrm{counts}\mathrm{s}^1\mathrm{arcmin}^2`$, respectively.
The smaller peak in Fig. 3 (left panel), located on the left of the projection at a distance of about $`90^{\prime \prime }`$ from the left border of the region, is caused by the presence of a QSO (Dressler et al. 1993) (see § 3.4).

### 3.2 Surface brightness profile

Although the cluster morphology is not spherically symmetric we perform a standard surface brightness analysis to obtain approximate values, which give a quantitative indication of the size of the cluster and also for comparison with other clusters. We generate an azimuthally averaged surface brightness profile for the cluster in the $`0.32.0\mathrm{keV}`$ energy band for each camera. This energy band was selected to optimize the signal$``$to$``$noise (S/N) ratio. Point sources were removed manually and the photons were binned into concentric annuli with a width of 3.3′′ (3 pixels of the MOS camera). The vignetting correction and the background subtraction was done as described in § 2.3 and § 2.4. The three profiles were then added.
As the cluster is clearly not spherically symmetric the profiles were centred on the positions of the two peaks ($`\mathrm{P}_{\mathrm{M1}}`$ and $`\mathrm{P}_{\mathrm{M2}}`$) of both subclusters. The profile centred on $`\mathrm{P}_{\mathrm{M2}}`$ is obtained excluding a sector of 120 in the direction of the other subcluster. In the case of the profile centred on $`\mathrm{P}_{\mathrm{M1}}`$ a wider sector was excluded (240) in order to avoid the emission coming from the region located at the southern-west part of M1. The resulting surface brightness profiles $`S(r)`$ are shown in Fig. 5 where the overall XMM point spread function (PSF) binned as the observed profiles is plotted as a dashed line. The overall PSF is obtained by adding the PSF of each camera, estimated at an energy of 1 keV and weighted by the respective cluster count rate in the $`0.32.0\mathrm{keV}`$ energy band. The comparison of the surface brightness profile data with the PSF shows that the effect of the PSF on the profile is negligible.

After the subtraction of the corresponding blank field profile we find no residual CXB background left (as already mentioned in Sec. 2.4). We therefore fit the surface brightness with a $`\beta `$-model without the need of an additional background component:

$`S(r)=S_0\left[1+\left(\frac{r}{r_c}\right)^2\right]^{3\beta +1/2}`$.

where $`S_0`$ is the central surface brightness, $`r_c`$ is the core radius, $`\beta `$ is the slope parameter. The profiles for the two subclusters M1 and M2 are fit out to radii of 3.5 and 3.2 respectively, at which the $`S/N>3\sigma `$. Results of the different fits are given in Table 1. We obtain similar fit parameters for the two subclusters; however the central region of the profile centred on $`\mathrm{P}_{\mathrm{M1}}`$ is not well fitted probably because of internal structure of the subcluster that flattens the central part. The derived beta and core radii parameters for both subclusters are larger compared with the ones obtained by Schindler et al. (1998). These discrepancies could be due to the low S/N ratio of the HRI observation that provides artificially low values for the beta-parameters. The HRI profiles decrease rapidly to background level (see Fig. 2 in Schindler et al. (1998)), therefore a beta model fit is actually probing only the central part of the profile. In order to test this possibility decreasing cut off radii were applied to the XMM profiles. Smaller values of the core radius and beta were indeed derived: for a cut off radius of 2’ $`r_c=0.5^{}`$ and $`\beta =0.41`$ were obtained for the profile centred on $`\mathrm{P}_{\mathrm{M1}}`$ and $`r_c=0.8^{}`$ and $`\beta =0.57`$ for the profile centred on $`\mathrm{P}_{\mathrm{M2}}`$. For a more detailed analysis of the substructure, we subtract the emission coming from the subcluster M2 from the original image by using the $`\beta `$-model with the fit parameters in Table 1. The residual image is shown in Fig. 6. In this residual image the emission from the quasar and from the point sources that appear in Fig. 4 is still visible. There is clearly extended emission left from the subcluster M1, but its peak has shifted to the south-west of about 25 arcsec. North of M1 there is also residual emission surrounding some point sources, which is clearly extended and structured. South-west of M1, separated by an artifact minimum caused by the gap in the detector chips is a last piece of residual extended emission.

### 3.3 Count rate

The X-ray emission of the cluster can be traced out to $`3.5`$ arcminutes ($`1.4\mathrm{Mpc}`$). The MOS and pn count rates within the $`R3.5^{}`$ region and in the considered energy band (0.3-2.0 keV) are $`0.125\pm 0.013`$ counts s-1 and $`0.31\pm 0.03`$ counts s-1 respectively.

### 3.4 Other sources in the pointing

In the pointing many other sources are found around the cluster. $`19`$ point sources, located in the area covered by MOS CCD$`\mathrm{\#}1`$, were detected by both the pipeline detection and a peak density finding algorithm. Their coordinates, extension and count rates (after background subtraction) are listed in Table 2. Most of these sources are shown in Fig. 3 (right panel) together with their identification number.

For the sources corresponding to member galaxies, or located within $`4.3\mathrm{arcmin}`$ from the cluster centre, a more detailed analysis has been performed.

For this source no optical counterpart was found. The low statistics of the data do not allow us to constrain the photon index, even using a simple power law model. Assuming a power law spectrum with photon index $`\mathrm{\Gamma }=1.7`$ we obtain an unabsorbed flux of $`f_\mathrm{X}(0.310.0)=5.7\times 10^{15}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$.

No optical counterpart. Assuming a power law spectrum with photon index $`\mathrm{\Gamma }=1.7`$: $`f_\mathrm{X}(0.310.0)=3.5\times 10^{15}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$.

No optical counterpart. Assuming a power law spectrum with photon index $`\mathrm{\Gamma }=1.7`$: $`f_\mathrm{X}(0.310.0)=2.2\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$.

Probable member galaxy Abell 0851:\[DG92\] 124. Assuming a power law spectrum with photon index $`\mathrm{\Gamma }=1.7`$: $`f_\mathrm{X}(0.310.0)=1.0\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, and a bolometric luminosity of $`L_\mathrm{X}(\mathrm{bol})=1.6\times 10^{43}\mathrm{erg}\mathrm{s}^1`$.

Member galaxy Abell 0851:\[DG92\] 311 at a redshift of $`z=0.4007`$. It is also a radio source known as \[CDS90\] R082. Its spectrum was fitted with a power law giving a photon index $`\mathrm{\Gamma }=1.6\pm 0.2`$; the unabsorbed flux is $`f_\mathrm{X}(0.310.0)=4.7\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, and the bolometric luminosity $`L_\mathrm{X}(\mathrm{bol})=7.5\times 10^{43}\mathrm{erg}\mathrm{s}^1`$.

Probable member galaxy Abell 0851:\[SED95\] 074. The spectral fit yields $`\mathrm{\Gamma }=1.3\pm 0.1`$; $`f_\mathrm{X}(0.310.0)=6.5\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$; $`L_\mathrm{X}(\mathrm{bol})=1.1\times 10^{44}\mathrm{erg}\mathrm{s}^1`$.

This source was already detected during a ROSAT PSPC observation as RXJ0943.0+4701 by Schindler & Wambsganss (1996) and had a measured flux in the ROSAT $`\left(0.12.4\mathrm{keV}\right)`$ band of $`f_{\mathrm{ROSAT},\mathrm{PSPC}}=4\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$. During a second ROSAT observation of the cluster done by Schindler et al. (1998) using the HRI, the source was hardly visible with a resulting flux in the $`\left(0.12.4\mathrm{keV}\right)`$ band of $`f_{\mathrm{ROSAT},\mathrm{HRI}}=5\times 10^{15}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, indicating that this source has a very strong variability. The source coordinates observed by PSPC are slightly offset ($`38^{\prime \prime }`$ south) with respect to the coordinates observed with XMM.
We find the galaxy Abell 0851:\[SED95\] 066 as a possible optical counterpart, which is a probable member galaxy. The spectral fit yields $`\mathrm{\Gamma }=1.8\pm 0.1`$; $`f_\mathrm{X}(0.310.0)=6.8\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$; $`L_\mathrm{X}(\mathrm{bol})=1.2\times 10^{44}\mathrm{erg}\mathrm{s}^1`$. We compute a flux in the $`\left(0.12.4\mathrm{keV}\right)`$ band of $`f_{\mathrm{XMM}}=4.7\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, which is only slightly higher than the values computed by Schindler & Wambsganss (1996) in their first X-ray observation of the cluster. Considering the changing values of the luminosity and the results of the spectral analysis, the conclusion drawn by Schindler et al. (1998), of source $`\mathrm{\#}11`$ being a strong X-ray variable AGN, is, up to now, the most reliable one.

No optical counterpart was found for this source. The spectral fit yields $`\mathrm{\Gamma }=1.95\pm 0.2`$; $`f_\mathrm{X}(0.310.0)=2.0\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$.

No optical counterpart. The spectral fit yields $`\mathrm{\Gamma }=1.7_{0.3}^{+0.2}`$; $`f_\mathrm{X}(0.310.0)=2.5\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$.

Probable member galaxy Abell 0851:\[DG92\] 292. The spectral fit yields $`\mathrm{\Gamma }=1.6\pm 0.1`$; $`f_\mathrm{X}(0.310.0)=3.3\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$; $`L_\mathrm{X}(\mathrm{bol})=5.4\times 10^{43}\mathrm{erg}\mathrm{s}^1`$.

Probable member galaxy Abell 0851:\[DG92\] 064. Assuming a power law spectrum with photon index $`\mathrm{\Gamma }=1.7`$ we obtain an unabsorbed flux of $`f_\mathrm{X}(0.310.0)=1.8\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, and a bolometric luminosity of $`L_\mathrm{X}(\mathrm{bol})=2.9\times 10^{43}\mathrm{erg}\mathrm{s}^1`$.

No optical counterpart was found. Assuming a power law spectrum with photon index $`\mathrm{\Gamma }=1.7`$: $`f_\mathrm{X}(0.310.0)=3.7\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$.

This source is a QSO (known as Abell 0851:\[DG92\] 440) lying in the field of view of the cluster and at a redshift of $`z=2.055`$ (Dressler et al. 1993). The spectral fit yields $`\mathrm{\Gamma }=1.72\pm 0.05`$; $`f_\mathrm{X}(0.310.0)=3.1\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$; $`L_\mathrm{X}(\mathrm{bol})=2.3\times 10^{45}\mathrm{erg}\mathrm{s}^1`$. In the ROSAT energy band the QSO has a luminosity of $`L_{\mathrm{XMM}}\left(0.12.4\right)=8.6\times 10^{44}\mathrm{erg}\mathrm{s}^1`$, which is slightly lower than was measured by Schindler et al. (1998) ($`L_{\mathrm{ROSAT}}\left(0.12.4\right)=1.4\times 10^{45}\mathrm{erg}\mathrm{s}^1`$). Due to the low number of counts Schindler et al. (1998) could not fit any spectral model to the data and, to compute the luminosity, they fixed the photon index at $`\mathrm{\Gamma }=2.3`$. Using the same value for the photon index the two luminosities, computed by ROSAT and XMM, are perfectly consistent.

Probable member galaxy Abell 0851:\[DG92\] 105. Assuming a power law spectrum with photon index $`\mathrm{\Gamma }=1.7`$ we obtain an unabsorbed flux of $`f_\mathrm{X}(0.310.0)=1.9\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, and a bolometric luminosity of $`L_\mathrm{X}(\mathrm{bol})=3.1\times 10^{43}\mathrm{erg}\mathrm{s}^1`$.

J0943.2+4700; has no optical counterpart. The spectral fit yields $`\mathrm{\Gamma }=1.8\pm 0.1`$; $`f_\mathrm{X}(0.310.0)=7.8\times 10^{14}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$.

### 3.5 Comparison with the optical waveband

In Fig. 7 an HST image of the central area of the cluster is shown. On it, superposed, are the X-ray contours.
An extensive redshift survey of the galaxies in the field of view of CL 0939+4713 was performed by Dressler et al. (1999) and at least $`60`$ galaxies were spectroscopically proven to belong to the cluster. Among these galaxies, those falling within the field of view of HST, can be seen in Fig. 7 surrounded by rectangles.
Some of these galaxies (e.g. $`\mathrm{\#}9`$) are very luminous also in X-rays but, nevertheless, the X-ray luminosity of the central region of the cluster is strongly dominated by extended hot emission from the intra-cluster gas, peaked in the two subclusters.
Both the galaxy number density distribution and the X-ray emission show two main peaks. The peaks of these two distributions are in similar regions. It is also interesting that the mass distribution (as computed by Seitz et al. (1996) using a weak lensing analysis) has a similar bimodal structure.
The subcluster M1 coincides with the strongest maxima in both the galaxy number density and the mass distribution. The peaks of the galaxy density and of the mass distributions are $`26\mathrm{arcsec}`$ West and $`15\mathrm{arcsec}`$ South of the X-ray peak M1.
The other subcluster, M2, is located near the second maxima in the galaxy number density and in the mass distribution (their peaks are, respectively, $`16\mathrm{arcsec}`$ west and $`20\mathrm{arcsec}`$ south-east of the X-ray one).

## 4 Spectral analysis

For the whole spectral analysis, the ready-made on-axis response files (m1\_r5\_all\_15.rsp for MOS1, m2\_r5\_all\_15.rsp for MOS2 333Both MOS files are for patterns from $`0`$ to $`12`$. and epn\_ff20\_sY9.rmf for pn 444pn response file is for single-pixel events only.) were used. 555All response files can be downloaded via anonymous ftp at xmm.vilspa.esa.es in the directory /pub/ccf/constituents/extras/responses.
Source and background spectra were accumulated separately for each detector. Spectra were re-binned in order to obtain $`S/N>5\sigma `$ in each bin after background subtraction.
As the spectra are already corrected for vignetting, spectra of the same physical region, observed with the three different detectors, can simply be fitted simultaneously in order to maximize the signal to noise ratio.
All the spectral analysis was performed in the $`0.310.0\mathrm{keV}`$ band. Data below $`0.3\mathrm{keV}`$ have been excluded to avoid residual calibration problems in the MOS & pn response matrices below these energies.

### 4.1 Results

The spectra for the cluster analysis were extracted from a circular region centered at $`\mathrm{RA}=09^\mathrm{h}43^\mathrm{m}01.0^\mathrm{s}`$, $`\mathrm{DEC}=+46\text{}59{}_{}{}^{}37.5^{\prime \prime }`$ (J2000) and of radius r=$`3.0^{}`$. The emission coming from all the point sources inside the region was excluded.
X-ray temperatures, metal abundances, fluxes and luminosities, were computed with XSPEC fitting the data to an absorbed isothermal plasma emission code of Kaastra & Mewe (1993) including the FeL calculations of Liedahl et al. (1995) 666MeKaL code from now on., folded through the appropriate response matrices. In Table 3 is a detailed description of the results obtained fitting the theoretical model to data from the three detectors at the same time, and for each detector independently.
Temperatures obtained using MOS and pn detectors independently are largely consistent. The fit results for the redshift are consistent with its spectroscopic value computed by Dressler et al. (1999) ($`z=0.406`$); good agreement is also observed between the fit value of the hydrogen column density and its Galactic value of $`n_\mathrm{H}=1.24\times 10^{20}\mathrm{cm}^2`$ (Dickey & Lockman 1990). When fixing these two parameters the best fit temperature and metal abundance of the cluster based on a simultaneous fit of all (MOS+pn) the spectra (using independent normalizations for the MOS and the pn cameras) are $`kT=4.5\pm 0.2\mathrm{keV}`$ and $`Z=0.20\pm 0.06Z_{}`$, respectively (the corresponding MOS1, MOS2 & pn spectra with the fitted MeKaL model superimposed are plotted in Fig. 8).
Using the best fit model parameters we derived, out to a radius of $`210^{\prime \prime }`$, an unabsorbed flux of $`f_\mathrm{X}(0.310.0)=1.3\pm 0.1\times 10^{12}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, corresponding to a luminosity, in the cluster rest frame, of $`L_\mathrm{X}(0.310.0)=1.1\pm 0.1\times 10^{45}\mathrm{erg}\mathrm{s}^1`$ and a bolometric luminosity of $`L_\mathrm{X}(\mathrm{bol})=1.3\pm 0.1\times 10^{45}\mathrm{erg}\mathrm{s}^1`$.

### 4.2 Temperature map

Since the cluster has been observed to have an irregular morphology, a deeper analysis of its internal structure and dynamical state is required in order to gain a better understanding of it.
The cluster was subdivided into $`18`$ regions; their sizes have been chosen such that each region (after background subtraction) has similar number of photon counts. We fixed the HI column density and the redshift to the average values for the whole cluster and then fitted the spectrum for each region to obtain the values of the temperature.

The resulting temperature map is shown in Fig. 9, where different shades of gray reflect the different temperatures of the various regions. Superposed are linearly spaced X-ray contours.
A clear trend of higher temperatures is observed in the region situated in the centre between the two clumps forming the core. The region of higher temperature also extends in the NW and SE directions.
Average temperatures were also computed from spectra extracted from the central hotter region and from the two regions surrounding it. All three regions have elliptical shapes and the coordinates of their centres, and the major and minor axes are: M2. RA$`=9^\mathrm{h}43^\mathrm{m}08.11^\mathrm{s}`$, DEC$`=+47\text{}00{}_{}{}^{}03.48^{\prime \prime }`$ (J2000), $`a=1.5^{}`$, $`b=1.2^{}`$, P.A.$`=140\text{}`$; Central. RA$`=9^\mathrm{h}43^\mathrm{m}01.58^\mathrm{s}`$, DEC$`=+46\text{}59{}_{}{}^{}18.02^{\prime \prime }`$ (J2000), $`a=2.1^{}`$, $`b=0.4^{}`$, P.A.$`=159\text{}`$; M1. RA$`=9^\mathrm{h}42^\mathrm{m}55.50^\mathrm{s}`$, DEC$`=+46\text{}58{}_{}{}^{}35.50^{\prime \prime }`$ (J2000), $`a=1.4^{}`$, $`b=1.0^{}`$, P.A.$`=179\text{}`$. A MeKaL model was used; HI column density and the redshift are fixed to the average values for the whole cluster. The resulting temperatures, which clearly confirm the trend observed in the temperature map, are $`kT=4.8_{0.4}^{+0.5}\mathrm{keV}`$, $`kT=6.0_{0.6}^{+0.7}\mathrm{keV}`$ and $`kT=4.7_{0.3}^{+0.4}\mathrm{keV}`$ (for M2, central and M1 regions, respectively). Detailed results of the fits are presented in Table 4.

The temperature map in Fig. 9 is similar to the predictions of simulations of cluster merging, for an almost head-on encounter of two subclusters of similar mass, shortly before the merging process (Ricker & Sarazin 2001; Ritchie & Thomas 2002; Takizawa 1999). Similar temperature structure has also been observed in clusters such as A399 and A401 (Markevitch et al. 1998) and A3528, A1750, and A3395 (Donnelly et al. 2001).
The central hotter region therefore delineates an area where the two subclusters are compressing and heating intra-cluster gas, while moving one towards the other in the beginning of a merging process. Simulations by Schindler & Müller (1993) show that this hotter region between the two merging structures begins to appear $`0.51.0\mathrm{Gyr}`$ before the merging. Comparing this result with simulations performed by Takizawa (1999, 2000) allows a determination of the dynamical state: the two subclusters in CL 0939+4713 will collide in a few hundreds of $`\mathrm{Myr}`$.

### 4.3 Metallicity variation

A higher metal abundance is observed in region M1 compared to the abundance in M2 and to the abundance in the central region (see Table 4). While in M1 a metallicity of $`0.33_{0.12}^{+0.14}`$ is found, the metallicity in M2 amounts only to $`0.15\pm 0.11`$ in solar units. The errors are relatively large because only photons from small areas were used. The 90% errors overlap slightly. But there is a clear trend that the metallicity is not constant. The high-metallicity region M1 coincides with a strong peak in the galaxy number density distribution. It seems that the galaxies present in this region have enriched the intra-cluster gas with heavy elements. It is interesting that the high metal content is still localised within this region despite the active dynamical state of the cluster. In order to maintain such a high metallicity during the merger process, the subcluster M1, on its way towards the subcluster M2, did obviously, not accumulate a lot of gas from the front and did not leave large quantities of its gas behind, always retaining the same gas. An additional explanation could be that the gas in this region has been metal enriched only recently, so that the motion of the gas has not had time yet to distribute the metals.

## 5 Mass analysis

The major difficulty to estimate the total mass in this cluster is its active dynamical estate. CL 0939+4713 is undergoing a major merger and the hypothesis of hydrostatic equilibrium and spherical symmetry are in principle not valid. We therefore instead assume a bimodal model where each subcluster follows an isothermal beta-model. To estimate the gas and total mass of the cluster we analyse each subcluster separately assuming for each of them hydrostatic equilibrium and spherical symmetry.
The integrated mass contained within a radius $`r`$ for each subcluster can be calculated from the equation

$$M(<r)=\frac{kr}{\mu m_\mathrm{p}G}T_{\mathrm{gas}}\left(\frac{d\mathrm{ln}\rho _{\mathrm{gas}}(r)}{d\mathrm{ln}r}+\frac{d\mathrm{ln}T_{\mathrm{gas}}(r)}{d\mathrm{ln}r}\right)$$
(1)
where $`\rho _{\mathrm{gas}}`$ and $`T_{\mathrm{gas}}`$ are the density and the temperature of the intracluster gas, and $`r`$,$`k`$,$`\mu `$,$`m_\mathrm{p}`$, and $`G`$ are the radius, the Boltzmann constant, the molecular weight, the proton mass, and the gravitational constant, respectively.
To estimate the gas and total mass for each subcluster, the $`\beta `$ parameters obtained in the fits of the profiles centred on the peaks $`\mathrm{P}_{\mathrm{M1}}`$ and $`\mathrm{P}_{\mathrm{M2}}`$ are used (see Table 1). We assume a constant temperature of 4.7 and 4.8 keV for the subclusters M1 and M2, respectively. The two subclusters are very close each other: the outer radius of $`22^{\prime \prime }`$ (140 kpc) corresponds to half of the distance between the subclusters. Therefore we can compute the gas mass and total cluster mass out to this radius adding the masses corresponding to each subcluster. Fig. 10 shows the integrated gas and total masses out to 140 kpc for the two subclusters independently and the numerical values are given in Table 5.

For an estimation of the total cluster mass out to a larger radius we use the fit centred on the subcluster M2, where the X-ray emission is higher. In Table 6 the gas and total masses found at two different radii are listed. At radius $`r_{500}=0.85\mathrm{Mpc}`$ the total cluster mass and gas mass are: $`M_{\mathrm{tot}}\left(r_{500}\right)=2.3\times 10^{14}M_{}`$ and $`M_{\mathrm{gas}}\left(r_{500}\right)=0.49\times 10^{14}M_{}`$. The gas mass fraction tends to increase slightly outwards ($``$17% at 0.5 Mpc and $``$22% at 1 Mpc).
The tabulated errors on the total mass come from the uncertainty in the temperature and the errors in the $`\beta `$ parameters. The emissivity of the gas in the XMM energy band $`0.32.0\mathrm{keV}`$ is almost constant inside the uncertainties for the estimate temperature. Therefore we can derive the gas density distribution without the uncertainty of the temperature estimate.

The mass distribution of part of CL 0939+4713 has been reconstructed with a weak lensing analysis by Seitz et al. (1996) and a comparison with our cluster mass estimation would be interesting. Unfortunately the restricted region within which the total mass has been calculated with the gravitational lensing analyses makes a proper comparison impossible.

## 6 Discussion of the merger hypothesis

Several observed characteristics of CL 0939+4713, either in the X-rays and in the optical, point to a scenario in which the two subclusters M1 and M2, forming the cluster core, with masses of 4.2$`\times 10^{12}M_{}`$ and 4.8$`\times 10^{12}M_{}`$ respectively, are in the course of a merging process, moving one towards the other along a direction coincident with the cluster major axis. The cores of the subclusters will pass through each other in a few Myr.

The spectral analysis of CL 0939+4713 has revealed an elongated region of heated gas located between the two X-ray peaks, and almost perpendicular to the cluster major axis. Hydrodynamic simulations (Schindler & Müller 1993; Ricker & Sarazin 2001; Ritchie & Thomas 2002; Takizawa 1999) show that the gas located between two subclusters is compressed and heated during the approach of the subclusters. This hotter gas tends to expand in the direction perpendicular to the collision axis, therefore it shows an elongated shape perpendicular to the collision axis. The gas temperature and the elongation of the central region increase as the two subclusters get closer, while the bulk of the X-ray emission still comes from the cores of the two subclusters until the two cores collide.

CL 0939+4713 shows an elongated X-ray morphology. Also this feature is consistent with simulations of merging clusters (Schindler 2002), which show that before the merger the intracluster gas is elongated along the collision axis.

The X-ray maxima of the two subclusters are displaced with respect to the peaks in the galaxy density distribution. In M1 the X-ray peak seems to be ahead of the galaxy concentration. Obviously, the galaxies do not trace the dark matter distribution well in this region. This is not surprising as during a merger the distributions of the different components can be quite different.

CL 0939+4713 has a high velocity dispersion, $`\sigma 1300\mathrm{km}\mathrm{s}^1`$ (Dressler & Gunn 1992; Dressler et al. 1999), indicating a relative motion of two or more components along the line of sight. The velocity dispersions of the two subclusters, computed independently, are instead slightly lower: $`\sigma 1000\mathrm{km}\mathrm{s}^1`$ and $`\sigma 750\mathrm{km}\mathrm{s}^1`$ for regions M2 and M1, respectively.

Dressler et al. (1999) observed a very high fraction of post-starburst galaxies, in CL 0939+4713, compared to other rich clusters and to field galaxies. Fujita et al. (1999) have investigated the effect of ram pressure in merging clusters. They find (in contrast to what was predicted previously by Evrard (1991)) that the increase in the star-formation rate caused by compression is less important than the ram-pressure stripping, which dramatically increases before a merging process. The stripping of the interstellar medium causes a decrease of the star formation rate, with a subsequent increase of the post starburst galaxy fraction.

Andreon et al. (1997) found that early-type galaxies in CL 0939+4713 do not show a radial symmetry, but are distributed along a privileged direction. This direction roughly coincides with the X-ray major axis and therefore with the possible collision axis.

## 7 Summary and Conclusions

We confirm earlier results that CL 0939+4713 is a dynamically young system (Schindler & Wambsganss 1996; Schindler et al. 1998). The X-ray morphology has pronounced substructure. New findings of a hot region between the two main clumps are particularly indicative of a merger process, in which the two major subclusters will collide in a few hundreds of $`\mathrm{Myr}`$.
Particularly interesting are the metallicity variations found in this cluster. Metallicity variations were found before as radial gradients in clusters with cD galaxies (De Grandi & Molendi 2001; Molendi & Gastaldello 2001; Kaastra et al. 2001; Gastaldello & Molendi 2002; Lewis et al. 2002). Clusters without cD galaxies were found to show no radial gradients (De Grandi & Molendi 2001; Arnaud et al. 2001b). The problem when binning the cluster into radial bins is, that one averages over metallicity values coming from very different regions of the clusters. Therefore it is not surprising that the profiles are usually flat for non-cD clusters. We took a different approach here and found a hint that the metal abundances vary over the cluster. Of course it is quite difficult to do this in a relatively distant cluster. We therefore suggest to use this method in nearby clusters, instead of determining radial metallicity profiles.

###### Acknowledgements.

We are very grateful to Bernd Aschenbach for making the XMM data available to us. They are part of the TS/MPE guaranteed time. We thank L.M. Buson and M. Cappellari for kindly providing the HST image (grazie mille!). A. Castillo-Morales acknowledges support by the Marie Curie Training Site grant HPMT-CT-2000-00136 by the European Commission. We also thank Dr. Phil A. James for carefully reading this paper. We also thank the referee Monique Arnaud: her suggestions and comments have greatly improved the reliability of the results obtained.