# 1 Introduction

## 1 Introduction

Charm quarks produced in deep-inelastic scattering have been identified in sizable numbers by the H1 and ZEUS collaborations at HERA , and considerably more charm (and bottom) data are anticipated. At the theoretical level the reaction has already been studied extensively. In the framework where the heavy quark is not treated as a parton, leading order (LO)  and next-to-leading order (NLO)  calculations of the inclusive structure functions exist. Moreover, LO (AROMA, RAPGAP)  and NLO (HVQDIS) Monte-Carlo programs, allowing a much larger class of observables to be compared with data, have been constructed in recent years. The NLO QCD description agrees quite well with the HERA data. Most of the recent theoretical attention for this reaction has been in the context of variable flavor number schemes , which we shall not address here. We shall review here two issues regarding heavy quark production in deep-inelastic scattering. In the next section we discuss two new features of HVQDIS, relevant to recent and future analyses; first, the inclusion of charmed-meson semi-leptonic decays, and second, a switch to describe (LO) $`D`$$`\overline{D}`$$`jet`$ final states. In the third section we review the methods and some key results of the next-to-leading logarithmic Sudakov resummation for DIS production of heavy quarks, and NNLO estimates derived from this resummation.

## 2 The NLO Monte-Carlo HVQDIS

The program HVQDIS  is based on the next-to-leading order fully differential heavy quark contributions to the proton structure functions . The basic components (in terms of virtual-photon-proton scattering) are the 2–to–2 body squared matrix elements through one-loop order and tree level 2–to–3 body squared matrix elements, for both photon-gluon and photon-light-quark initiated subprocesses. It is therefore possible to study single- and semi-inclusive production at next-to-leading order, and three body final states at leading order. The goal of the next-to-leading order calculation is to organize the soft and collinear singularity cancellations without loss of information in terms of observables that can be predicted.

The subtraction method provides a mechanism for this cancellation. It allows one to isolate the soft and collinear singularities of the 2–to–3 body processes within the framework of dimensional regularization without calculating all the phase space integrals in a space-time dimension $`n4`$. Expressions for the three-body squared matrix elements in the limit where an emitted gluon is soft appear in a factorized form where poles $`ϵ=2n/2`$ multiply leading order squared matrix elements. These soft singularities cancel upon addition of the interference of the leading order diagrams with the renormalized one-loop virtual diagrams. The remaining singularities are initial state collinear in origin wherein the three-body squared matrix elements appear in a factorized form, with poles in $`ϵ`$ multiply splitting functions convolved with leading order squared matrix elements. These collinear singularities are removed through mass factorization.

The result of this calculation is an expression that is finite in four-dimensional space time. One can compute all phase space integrations using standard Monte-Carlo integration techniques. The final result is a program which returns parton kinematic configurations and their corresponding weights, accurate to $`𝒪(\alpha \alpha _s^2)`$. The user is free to histogram any set of infrared-safe observables and apply cuts, all in a single histogramming subroutine. Additionally, one may study heavy hadrons using the Peterson et al. model. Detailed physics results from this program are given in .

Below we discuss and give examples of two new options of HVQDIS version 1.3<sup>1</sup><sup>1</sup>1available from harris@hep.anl.gov: electrons from semileptonic decays of the charmed hadron, and a switch for retaining only three body final states.

### 2.1 Semileptonic decays

HVQDIS has been extended to include the electron from the semileptonic decay of the charmed hadron. In the rest frame of the decaying hadron, the electrons are distributed isotropically. Their momentum distribution is the weighted average of multiple decay modes of many different charmed hadrons, and has been deduced  from RAPGAP. The parameterization is shown in fig. 1. The implementation used in HVQDIS1.3 (shown as points) was derived from the fit (dashed line) to the RAPGAP output (histogram). The overall normalization of the cross section is then fixed by the semileptonic branching ratio Br$`(ce)`$ which we take to be $`9.5\%`$.

The inclusive transverse momentum and pseudo-rapidity distributions of the semileptonic decay electrons in the lab frame in deep inelastic scattering of 820 GeV protons with 27.5 GeV electrons in the kinematic range $`0<y<0.7`$ and $`2<Q^2<100`$ GeV<sup>2</sup> are shown in fig. 2. We also show the corresponding distributions for the parent parton ($`c`$-quark) and hadron (charmed-meson). The curves are produced using the next-to-leading order Glück-Reya-Vogt 1994 (GRV94)  parton distribution functions, a two-loop $`\alpha _s`$ with $`n_f=3`$ and $`\mathrm{\Lambda }_{\mathrm{QCD}}^{(n_f=3)}=248`$ MeV, and $`m_c=1.35`$ GeV. The distributions of the charmed partons and hadrons are highly correlated because of the simple Peterson et al. fragmentation model. The semileptonic decay electrons are very soft, taking only a small portion of the hadron $`p_t`$, and more central due to the isotropic nature of the decay.

### 2.2 Three body final states

For speed considerations it pays to add a switch to turn off all two body contributions (primarily the very slow virtual routines) when one is interested only in a manifestly three body observable. Such a switch has been added to HVQDIS1.3. We give here a sample of three body observables.

The final state of interest is $`D`$$`\overline{D}`$jet corresponding to the partonic states $`c`$$`\overline{c}`$$`g`$ and $`c`$$`\overline{c}`$$`q`$. We begin by requiring the $`D`$, $`\overline{D}`$ and jet to be above some minimum transverse momentum ($`P_t^D>1.2`$ GeV, $`P_t^{\overline{D}}>1.2`$ GeV, $`P_t^{jet}>6`$ GeV), to be central ($`|\eta _D|<1.5`$, $`|\eta _{\overline{D}}|<1.5`$, $`|\eta _{jet}|<2.4`$), and to be isolated ($`R_{D\overline{D}}>0.7`$, $`R_{Djet}>0.7`$, $`R_{\overline{D}jet}>0.7`$) in the lab frame. The cone size $`R_{ij}=\sqrt{(\eta _i\eta _j)^2+(\varphi _i\varphi _j)^2}`$.

The total cross section for the deep inelastic production of $`D`$$`\overline{D}`$$`jet`$ as a function of their invariant mass $`M_{D\overline{D}j}`$ is shown in fig. 3 for $`0<y<0.7`$ and $`2<Q^2<100\mathrm{GeV}^2`$. A one loop $`\alpha _s`$ with $`n_f=3`$ and $`\mathrm{\Lambda }_{QCD}^{(n_f=3)}=232`$ MeV, leading order GRV94 parton distributions, and $`m_c=1.35`$ GeV where used. The normalization of this LO curve has a large uncertainty as demonstrated by the various scale choices $`\mu =\{\mu _0/2,\mu _0,2\mu _0\}`$, with $`\mu _0=\sqrt{Q^2+m_c^2+M_{D\overline{D}j}^2}`$. Also shown in the figure is a decomposition into the gluon and light-quark initiated subprocesses. The gluon initiated subprocess dominates due to the relatively large size of the gluon parton distribution function at small $`x`$. As another example, in the $`D\overline{D}j`$ center-of-mass frame we construct the Dalitz energy fractions $`x_i=2E_i/M_{D\overline{D}j}`$, ($`i=D,\overline{D}`$, or $`j`$) that specify how much available energy is shared between the $`D,\overline{D}`$, and jet. They satisfy $`x_D+x_{\overline{D}}+x_j=2`$. The normalized cross section differential in $`x_D`$ and $`x_j`$ is shown in fig. 4.

## 3 Soft-gluon resummation

As already remarked, existing NLO calculations for heavy quark electroproduction provide a solid theoretical perturbative QCD description  for this reaction, and agree well with present data . At moderate $`Q^2`$ and $`x`$ values larger than $`0.01`$, the charm structure function $`F_2^{\mathrm{charm}}`$ is increasingly dominated by partonic processes near the charm quark pair production threshold. The large size of the gluon density $`f_g(x,\mu )`$ for small momentum fractions $`x`$ gives relatively large weight to such processes . Although the QCD corrections at presently accessible $`x`$ values are moderate (about 30-40%), with an increasing number of data to be gathered at higher $`x`$, it is worthwhile to have a closer look at such near-threshold subprocesses. In this kinematic region, the QCD corrections are dominated by large Sudakov double logarithms. Recently , these Sudakov logarithms have been resummed to all orders of perturbation theory, to next-to-leading logarithmic (NLL) accuracy, and, moreover, in single-particle inclusive (1PI) kinematics <sup>2</sup><sup>2</sup>2Analytical results for pair-inclusive kinematics are also given in . Let us recall at this point the main results. First, the quality of the approximation for the next-to-leading logarithmic threshold resummation was found to be clearly superior to leading logarithmic one. Furthermore, the resummation provided NNLO estimates ,which were found to be sizable for $`x0.05`$.

Below we give a synopsis of the analysis presented in . We study electron proton scattering with the exchange of a single virtual photon, $`Q^2=q^2`$, and a detected heavy quark (we concentrate on the charm quark case here) in the final state, i.e. the subprocess

$`\gamma (q)+P(p)`$ $``$ $`\mathrm{Q}(p_1)+X[\overline{\mathrm{Q}}],`$ (1)   
where $`X`$ denotes any additional hadrons, including the heavy anti-quark, in the final state and $`p_1^2=m^2`$. The Mandelstam invariants, $`S^{}=(p+q)^2+Q^2,T_1=(pp_1)^2m^2`$ and $`U_1=(qp_1)^2m^2`$ can be used to define $`S_4=S^{}+T_1+U_1`$, which vanishes at the hadronic threshold. The double differential heavy quark structure function $`dF_2`$ associated to (1) may be written as

$$\frac{d^2F_{2,P}(x,S_4,T_1,U_1,Q^2,m^2)}{dT_1dU_1}=\frac{1}{S^2}\underset{i=q,\overline{q},g}{}\underset{z^{}}{\overset{1}{}}\frac{dz}{z}f_{i/P}(z,\mu ^2)\omega _{2,i}(\frac{x}{z},\frac{s_4}{\mu ^2},\frac{t_1}{\mu ^2},\frac{u_1}{\mu ^2},\frac{Q^2}{\mu ^2},\frac{m^2}{\mu ^2},\alpha _s(\mu )),$$
(2)   
where $`z^{}=U_1/(S^{}+T_1)`$. The $`f_{i/P}`$ denote parton distributions in the proton at momentum fraction $`z`$ and $`\overline{\mathrm{MS}}`$-mass factorization scale $`\mu `$. The functions $`\omega _{2,i}`$ describe the underlying hard parton scattering processes and depend on the partonic Mandelstam variables $`s^{},t_1,u_1`$ and $`s_4`$, which are derived from (1) by replacing the proton $`P`$ by a parton of momentum $`k=zp`$. At $`n`$-th order in perturbation theory, the gluonic hard part $`\omega _{2,g}`$ in eq.(2) typically depends on singular distributions $`\alpha _s^n[\mathrm{ln}^{2n1k}(s_4/m^2)/s_4]_+,k=0,1,\mathrm{}`$, that must be resummed. Contributions from light initial-quark states are neglected, as they are about 5% at NLO.

The resummation of threshold logarithms rests upon the factorization of the kinematics and dynamics of the relevant degrees of freedom near threshold . The dynamical factors involved can be each be assigned a kinematic weight $`w_i`$ that is defined to vanish at threshold. For $`dF_{2,P}`$ in eq.(2), the factorization of the kinematics implies that these weights sum to the overall inelasticity near threshold: $`S_4/m^2w_1(u_1)/m^2+w_s`$, with $`w_1=1z`$ and $`w_s=s_4/m^2`$. Correspondingly, the infrared regulated partonic structure function $`dF_{2,g}`$ factorizes into functions that individually organize contributions from these near-threshold degrees of freedom. Thus, there is here a function $`\psi _{g/g}`$ that sums the singular distributions from incoming collinear gluons, and a soft function $`S`$ that organizes those due to soft gluons not collinear to the incoming parton. Finally, there is a hard function $`H_{2,g}`$ incorporating only regular short-distances corrections. Replacing the proton in eq.(2) by a gluon, and passing to Laplace-moment space, $`\stackrel{~}{f}(N)=_0^{\mathrm{}}𝑑w\mathrm{exp}[Nw]f(w)`$, this gives 

$`\stackrel{~}{\omega }_{2,g}(N,{\displaystyle \frac{t_1}{\mu ^2}},{\displaystyle \frac{u_1}{\mu ^2}},{\displaystyle \frac{Q^2}{\mu ^2}},{\displaystyle \frac{m^2}{\mu ^2}})`$ $`=`$ $`H_{2,g}({\displaystyle \frac{t_1}{\mu ^2}},{\displaystyle \frac{u_1}{\mu ^2}},{\displaystyle \frac{Q^2}{\mu ^2}},{\displaystyle \frac{m^2}{\mu ^2}})\left[{\displaystyle \frac{\stackrel{~}{\psi }_{g/g}(N_u,p\zeta /\mu )}{\stackrel{~}{\varphi }_{g/g}(N_u,\mu )}}\right]\stackrel{~}{S}({\displaystyle \frac{m}{N\mu }},\zeta ),`$ (3)   
where $`\varphi _{g/g}`$ is the usual $`\overline{\mathrm{MS}}`$-distribution from mass factorization and $`N_u=N(u_1)/m^2`$. In moment space, the Sudakov logarithms appear as factors $`\alpha _s^n\mathrm{ln}^{2ni}N`$, with $`i=0,1`$ for NLL accuracy. The $`N`$-dependence in eq.(3) exponentiates for each function individually. All leading logarithms (LL) are exclusively contained in $`\stackrel{~}{\psi }_{g/g}`$, which is a gluon distribution at fixed fraction of $`p\zeta `$ and can be defined as an operator matrix element. It depends on a time-like vector $`\zeta =p_2/m`$ ($`p_2`$ is the heavy antiquark momentum). Its collinear poles are canceled by $`\varphi _{g/g}`$. The threshold logarithms in $`\stackrel{~}{\psi }_{g/g}`$ are resummed in analogy to the Drell-Yan process , while all scale dependence of $`\stackrel{~}{\psi }_{g/g}`$ and $`\stackrel{~}{\varphi }_{g/g}`$ is governed by renormalization group equations (RGE) with anomalous dimensions $`\gamma _\psi =2\gamma _g`$ and $`\gamma _{g/g}`$ .

The soft function $`S`$ requires renormalization, since it is defined as a composite operator, that connects Wilson lines in the direction of the scattering partons . Its RGE, $`\mu (d/d\mu )\mathrm{ln}\stackrel{~}{S}(N)=2\mathrm{Re}\mathrm{\Gamma }_S`$, resums all threshold logarithms in $`\stackrel{~}{S}`$. Its gauge dependence cancels precisely that of $`\psi _{g/g}`$. The soft anomalous dimension $`\mathrm{\Gamma }_S`$ is to order $`\alpha _s`$

$`\mathrm{\Gamma }_S(\alpha _s)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{\pi }}\{({\displaystyle \frac{C_A}{2}}C_F)(L_\beta +1){\displaystyle \frac{C_A}{2}}(\mathrm{ln}\left({\displaystyle \frac{(p\zeta )^2}{m^2}}\right)+\mathrm{ln}{\displaystyle \frac{4m^4}{t_1u_1}})\},`$ (4)   
with $`\beta =\sqrt{14m^2/s}`$ and $`L_\beta =(12m^2/s)\{\mathrm{ln}(1\beta )/(1+\beta )+\mathrm{i}\pi \}/\beta `$.

The final result for the hard scattering function $`\stackrel{~}{\omega }_{2,g}`$ in moment space resums all large logarithms in single-particle inclusive kinematics up to NLL accuracy. Combining the resummed $`\stackrel{~}{\psi }_{g/g}`$ with the integrated RGE for $`\stackrel{~}{S}`$, we obtain for $`\stackrel{~}{\omega }_{2,g}`$ 

$`\stackrel{~}{\omega }_{2,g}(N,{\displaystyle \frac{t_1}{\mu ^2}},{\displaystyle \frac{u_1}{\mu ^2}},{\displaystyle \frac{Q^2}{\mu ^2}},{\displaystyle \frac{m^2}{\mu ^2}})`$ $`=`$   
$`H_{2,g}({\displaystyle \frac{t_1}{m^2}},{\displaystyle \frac{u_1}{m^2}},{\displaystyle \frac{Q^2}{m^2}},1)\stackrel{~}{S}(1,\alpha _s({\displaystyle \frac{m}{N}}))\mathrm{exp}\left\{2{\displaystyle \underset{\mu }{\overset{m}{}}}{\displaystyle \frac{d\mu ^{}}{\mu ^{}}}\gamma _g\left(\alpha _s(\mu ^{})\right)\right\}`$   
$`\times \mathrm{exp}\left\{{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{dw}{w}}(1\mathrm{e}^{N_uw})\left[{\displaystyle \underset{w^2}{\overset{1}{}}}{\displaystyle \frac{d\lambda }{\lambda }}A_{(g)}(\alpha _s(\sqrt{\lambda }m))+{\displaystyle \frac{1}{2}}\nu _{(g)}(\alpha _s(wm))\right]\right\}`$   
$`\times \mathrm{exp}\left\{{\displaystyle \underset{m}{\overset{m/N}{}}}{\displaystyle \frac{d\mu ^{}}{\mu ^{}}}2\mathrm{Re}\mathrm{\Gamma }_S\left(\alpha _s(\mu ^{})\right)\right\}\mathrm{exp}\left\{2{\displaystyle \underset{\mu }{\overset{2p\zeta }{}}}{\displaystyle \frac{d\mu ^{}}{\mu ^{}}}\left(\gamma _g\left(\alpha _s(\mu ^{})\right)\gamma _{g/g}(N_u,\alpha _s(\mu ^{}))\right)\right\}.`$   
The second exponent gives the leading $`N`$-dependence of the ratio $`\stackrel{~}{\psi }_{g/g}/\stackrel{~}{\varphi }_{g/g}`$ with $`\nu _{(g)}(\alpha _s)=2C_A\alpha _s/\pi `$, $`A_{(g)}(\alpha _s)=C_A(\alpha _s/\pi )+(C_AK/2)(\alpha _s/\pi )^2`$ and $`K=C_A(67/18\pi ^2/6)5/9n_f`$ . For NLL Sudakov resummation, the product $`H_{2,g}S`$ on the first line of eq. (3) is determined from matching to the Born cross section at the scale $`\mu =m`$.

The resummed result for $`\stackrel{~}{\omega }_{2,g}`$ in eq. (3) may be used as a generating functional for fixed order approximate perturbation theory by re-expanding $`\stackrel{~}{\omega }_{2,g}`$ to NLO and NNLO and inverting the Laplace transform. After insertion of eq. (3) into eq. (2) and integration over $`T_1,U_1`$, we may expand the structure function as

$`F_2^{\mathrm{charm}}(x,Q^2,m^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _s(\mu )e_\mathrm{c}^2Q^2}{4\pi ^2m^2}}{\displaystyle \underset{ax}{\overset{1}{}}}𝑑zf_{g/P}(z,\mu ^2){\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}(4\pi \alpha _s(\mu ))^k{\displaystyle \underset{l=0}{\overset{k}{}}}c_{2,g}^{(k,l)}(\eta ,\xi )\mathrm{ln}^l{\displaystyle \frac{\mu ^2}{m^2}},`$ (6)   
where $`a=(Q^2+4m^2)/Q^2`$ and $`e_c=2/3`$.

The quality of the NLL approximation eq.(3) can then be investigated by comparing exact, LL and NLL approximation to the gluon coefficient functions $`c_{2,g}^{(k,l)}`$, which is done in fig. 5(a). The functions $`c_{2,g}^{(k,l)}`$ depend on the scaling variables

$`\eta `$ $`=`$ $`{\displaystyle \frac{s}{4m^2}}1,\xi ={\displaystyle \frac{Q^2}{m^2}},`$ (7)   
where $`\eta `$ is a direct measure of the distance to the partonic threshold.

Fig. 5(a) reveals that, although at one loop the LL accuracy provides a good approximation for very small $`\eta `$, the NLL approximation is excellent over a much wider range in $`\eta `$, up to values of about 10 (the same actually holds true for the $`c_{2,g}^{(k,l)},l>0,k2`$). We also show $`c_{2,g}^{(2,0)}`$, which has more structure than in the $`c_{2,g}^{(1,0)}`$ curves.

To address the threshold sensitivity of the integrated charm structure function to threshold processes, we perform the integral over $`z`$ in eq. (2) only up to a value $`z_{\mathrm{max}}`$, and plot the integral then as a function of $`z_{\mathrm{max}}`$. In this way we can see where the integral eq. (2) acquires most of its value. The physical structure function corresponds to $`z_{\mathrm{max}}=1`$. Fig. 5(b) demonstrates that for $`x=0.01`$ $`F_2^{\mathrm{charm}}`$ is mostly determined by partonic processes close to threshold.

In fig. 6a we display at a fixed value of the factorization scale $`\mu =m`$ over a range of $`x`$, $`0.003x0.3`$ the effect of the NNLO corrections. We plot the K-factors $`F_{2(NNLO)}^{\mathrm{charm}}/F_{2(NLO)}^{\mathrm{charm}}`$ and, for comparison, also $`F_{2(NLO)}^{\mathrm{charm}}/F_{2(LO)}^{\mathrm{charm}}`$ <sup>3</sup><sup>3</sup>3For $`F_{2(LO)}^{\mathrm{charm}}`$ we used a two-loop $`\alpha _s`$ and NLO gluon density.. At NNLO we have taken the improved NLL approximation to $`F_2^{\mathrm{charm}}`$ (the exact NLO result plus the NLL approximate NNLO result). We see that particularly for smaller $`x`$, the size of the NNLO corrections is negligible, the K-factor being close to one, whereas for larger $`x`$, their overall size is still quite big, almost a factor of 2 at $`x=0.1`$.

Finally, in fig. 6b we show the NLO results as a function of $`p_T`$ for $`x=0.01`$, $`m=1.6\mathrm{GeV}`$, $`Q^2=10\mathrm{GeV}`$ and $`\mu =\sqrt{Q^2+4(m^2+p_T^2)}`$. At NLO, we compare our LL and NLL approximate results with the exact results of the second Ref. in . We see that the exact curves are reproduced well both in shape and magnitude by our NLL approximations, whereas the curves for LL accuracy systematically underestimate the true result. We also display the improved NLL approximation to $`dF_2^{\mathrm{charm}}/dp_T`$ at NNLO, which contains sizeable contributions to the value of the maximum increases by $`40\%`$ \- $`50\%`$.

## 4 Conclusions

Driven by the ever increasing variety and quantity of deep-inelastic charm production data from the H1 and ZEUS experiments at HERA, we have updated and reviewed two important tools: soft gluon threshold resummation and the next-to-leading order Monte-Carlo HVQDIS. The addition of semileptonic decay electron information and an option for only three body final states to the Monte-Carlo will enhance future physics analysis options. Soft gluon threshold resummation, on the other hand, teaches us about the size of the terms neglected in fixed order calculations.

### Acknowledgments

The work of B.W.H. was supported by the U. S. Department of Energy, High Energy Physics Division, Contract No. W-31-109-Eng-38. The work of J.S. was supported in part by the U. S. National Science Foundation grant PHY-9722101. The work of S.M. and E.L. is part of the research program of the Foundation for Fundamental Research of Matter (FOM) and the National Organization for Scientific Research (NWO).