format: text

text:
# Calculation of 𝜀'/𝜀

## 1 Introduction

The results obtained in the last two years by the NA48 NA48  and the KTeV KTeV  collaborations have marked a great experimental achievement, establishing some 35 years after the discovery of CP violation in the neutral kaon system Christenson  the existence of a much smaller violation acting directly in the decays:

$$\text{Re}(\epsilon ^{}/\epsilon )=\{\begin{array}{cc}(15.3\pm 2.6)\times 10^4\hfill & (\mathrm{NA48})\hfill \\ (20.7\pm 2.8)\times 10^4\hfill & (\mathrm{KTeV}).\hfill \end{array}$$
(1)
The average of these results with the previous measurements by the NA31 collaboration at CERN and by the E731 experiment at Fermilab gives

$$\text{Re}(\epsilon ^{}/\epsilon )=(17.2\pm 1.8)\times 10^4.$$
(2)
While the Standard Model (SM) of strong and electroweak interactions provides an economical and elegant understanding of indirect ($`\epsilon `$) and direct ($`\epsilon ^{}`$) CP violation in term of a single phase, the detailed calculation of the size of these effects implies mastering strong interactions at a scale where perturbative methods break down. In addition, direct CP violation in $`K\pi \pi `$ decays arises from a detailed balance of two competing sets of contributions, which may hopelessly inflate the uncertainties related to the relevant hadronic matrix elements in the final outcome. All that makes predicting $`\epsilon ^{}/\epsilon `$ a complex and challenging task review .

Just from the onset of the calculation the presence in the definition of $`\epsilon ^{}/\epsilon `$ , written as

$$\frac{\epsilon ^{}}{\epsilon }=\frac{1}{\sqrt{2}}\left\{\frac{(\pi \pi )_{I=2}|_W|K_L}{(\pi \pi )_{I=0}|_W|K_L}\frac{(\pi \pi )_{I=2}|_W|K_S}{(\pi \pi )_{I=0}|_W|K_S}\right\},$$
(3)
of given ratios of isospin amplitudes warns us of a longstanding and still unsolved theoretical “problem”: the explanation of the $`\mathrm{\Delta }I=1/2`$ selection rule.

The $`\mathrm{\Delta }I=1/2`$ selection rule in $`K\pi \pi `$ decays is known since 45 years Pais-Gell-Mann  and it states the experimental evidence that kaons are 400 times more likely to decay in the $`I=0`$ two-pion state than in the $`I=2`$ component ($`\omega A_2/A_01/22`$). This rule is not justified by any symmetry argument and, although it is common understanding that its explanation must be rooted in the dynamics of strong interactions, there is no up to date derivation of this effect from first principle QCD.

Given the possibility that common systematic uncertainties may a-priori affect the calculation of $`\epsilon ^{}/\epsilon `$ and the $`\mathrm{\Delta }I=1/2`$ rule (see for instance the present difficulties in calculating on the lattice the “penguin contractions” for CP violating as well as for CP conserving amplitudes Romapost ) a convincing calculation of $`\epsilon ^{}/\epsilon `$ must involve at the same time a reliable explanation of the $`\mathrm{\Delta }I=1/2`$ selection rule. Both observables indicate the need of large corrections to factorization in the evaluation of the four-quark hadronic transitions. Among these corrections Final State Interactions (FSI) play a substantial role. However, FSI alone are not enough to account for the large ratio of the $`I=0`$ over $`I=2`$ amplitudes. Other sources of large non-factorizable corrections are therefore needed for the CP conserving amplitudes, which might affect the determination of $`\epsilon ^{}/\epsilon `$ as well. As a consequence, a self-contained calculation of $`\epsilon ^{}/\epsilon `$ should also address the determination of the CP conserving $`K\pi \pi `$ amplitudes.

## 2 OPE: an “effective” approach

The Operator Product Expansion (OPE) provides us with a very effective way to address the calculation of hadronic transitions in gauge theories. The integration of the “heavy” gauge and matter fields allows us to write the relevant amplitudes in terms of the hadronic matrix elements of effective quark operators and of the corresponding Wilson coefficients (at a scale $`\mu `$), which encode the information about those dynamical degrees of freedom which are heavier than the chosen renormalization scale. With the SM flavor structure $`\mathrm{\Delta }S=1`$ transitions are effectively described by

$$_{\mathrm{\Delta }S=1}=\frac{G_\mathrm{F}}{\sqrt{2}}V_{ud}V_{us}^{}\underset{i}{}\left[z_i(\mu )+\tau y_i(\mu )\right]Q_i(\mu ).$$
(4)
The entries $`V_{ij}`$ of the $`3\times 3`$ Cabibbo-Kobayashi-Maskawa (CKM) matrix describe the flavour mixing in the SM and $`\tau =V_{td}V_{ts}^{}/V_{ud}V_{us}^{}`$. For $`\mu <m_c`$ ($`q=u,d,s`$), the relevant quark operators are:

$$\begin{array}{ccc}\begin{array}{ccc}Q_1\hfill & =& \left(\overline{s}_\alpha u_\beta \right)_{\mathrm{V}\mathrm{A}}\left(\overline{u}_\beta d_\alpha \right)_{\mathrm{V}\mathrm{A}}\hfill \\ Q_2\hfill & =& \left(\overline{s}u\right)_{\mathrm{V}\mathrm{A}}\left(\overline{u}d\right)_{\mathrm{V}\mathrm{A}}\hfill \end{array}\}\hfill & & \text{Current-Current}\hfill \\ \begin{array}{ccc}Q_{3,5}\hfill & =& \left(\overline{s}d\right)_{\mathrm{V}\mathrm{A}}_q\left(\overline{q}q\right)_{\mathrm{V}\mathrm{A}}\hfill \\ Q_{4,6}\hfill & =& \left(\overline{s}_\alpha d_\beta \right)_{\mathrm{V}\mathrm{A}}_q(\overline{q}_\beta q_\alpha )_{\mathrm{V}\mathrm{A}}\hfill \end{array}\}\hfill & & \text{Gluon “penguins”}\hfill \\ \begin{array}{ccc}Q_{7,9}\hfill & =& \frac{3}{2}\left(\overline{s}d\right)_{\mathrm{V}\mathrm{A}}_q\widehat{e}_q\left(\overline{q}q\right)_{\mathrm{V}\pm \mathrm{A}}\hfill \\ Q_{8,10}\hfill & =& \frac{3}{2}\left(\overline{s}_\alpha d_\beta \right)_{\mathrm{V}\mathrm{A}}_q\widehat{e}_q(\overline{q}_\beta q_\alpha )_{\mathrm{V}\pm \mathrm{A}}\hfill \end{array}\}\hfill & & \text{Electroweak “penguins”}\hfill \end{array}$$
(5)
Current-current operators are induced by tree-level W-exchange whereas the so-called penguin (and “box”) diagrams are generated via an electroweak loop. Only the latter “feel” all three quark families via the virtual quark exchange and are therefore sensitive to the weak CP phase. Current-current operators control instead the CP conserving transitions. This fact suggests already that the connection between $`\epsilon ^{}/\epsilon `$ and the $`\mathrm{\Delta }I=1/2`$ rule is by no means a straightforward one.

Using the effective $`\mathrm{\Delta }S=1`$ quark Hamiltonian we can write $`\epsilon ^{}/\epsilon `$ as

$$\frac{\epsilon ^{}}{\epsilon }=e^{i\varphi }\frac{G_\mathrm{F}\omega }{2|ϵ|\text{Re}A_0}\text{Im}\lambda _t\left[\mathrm{\Pi }_0\frac{1}{\omega }\mathrm{\Pi }_2\right]$$
(6)
where

$$\begin{array}{ccc}\mathrm{\Pi }_0\hfill & =& \frac{1}{\mathrm{cos}\delta _0}_iy_i\text{Re}Q_i_0(1\mathrm{\Omega }_{\mathrm{IB}})\hfill \\ \mathrm{\Pi }_2\hfill & =& \frac{1}{\mathrm{cos}\delta _2}_iy_i\text{Re}Q_i_2,\hfill \end{array}$$
(7)
and $`Q_i\pi \pi |Q_i|K`$. The rescattering phases $`\delta _{0,2}`$ can be extracted from elastic $`\pi \pi `$ scattering data FSIphases  and are such that $`\mathrm{cos}\delta _00.8`$ and $`\mathrm{cos}\delta _21`$. Given that the phase of $`\epsilon `$ ($`\theta _\epsilon `$) is approximately $`\pi /4`$, as well as $`\delta _0\delta _2`$, $`\varphi =\frac{\pi }{2}+\delta _2\delta _0\theta _\epsilon `$ turns out to be consistent with zero. While $`G_\mathrm{F}`$, $`\omega `$, $`|\epsilon |`$ and $`\text{Re}A_0`$ are precisely determined by experimental data, the first source of uncertainty that we encounter in eq. 6 is the value of $`\text{Im}\lambda _t\text{Im}(V_{ts}^{}V_{td})`$, the combination of CKM elements which measures CP violation in $`\mathrm{\Delta }S=1`$ transitions. The determination of $`\text{Im}\lambda _t`$ depends on B-physics constraints and on $`\epsilon `$ Ciuchinietal . In turn, the fit of $`\epsilon `$ depends on the theoretical determination of $`B_K`$, the $`\overline{K}^0K^0`$ hadronic parameter, which should be self-consistently determined within every analysis. The theoretical uncertainty on $`B_K`$ was in the past the main component of the final uncertainty on $`\text{Im}\lambda _t`$. The improved determination of the unitarity triangle coming from B-factories and hadronic colliders Lacker  has weakened and will eventually lift the dependence of $`\text{Im}\lambda _t`$ on $`B_K`$, allowing for an experimental measurement of the latter from $`\epsilon `$ . Within kaon physics, the decay $`K_L\pi ^0\nu \overline{\nu }`$ gives the cleanest “theoretical” determination of $`\text{Im}\lambda _t`$, albeit representing a great experimental challenge. At present, a typical range of values for $`\text{Im}\lambda _t`$ is $`(0.941.60)\times 10^4`$ BurasLP01 .

We come now to the quantities in the square brackets. While the calculation of the Wilson coefficients is well under control, thanks primarily to the work done in the early nineties by the Munich MunichNLO  and Rome RomaNLO  groups, the evaluation of the “long-distance” factors in eq. 7 is the crucial issue for the ongoing calculations. The isospin breaking (IB) parameter $`\mathrm{\Omega }_{\mathrm{IB}}`$, gives at the leading-order (LO) in the chiral expansion a $`positive`$ correction to the $`A_2`$ amplitude (proportional to $`A_0`$ via the $`\pi ^0\eta `$ mixing) of about 0.13 LOpieta . At the next-to-leading order (NLO) the full inclusion of the $`\pi ^0\eta \eta ^{}`$ mixing lift the value of $`\mathrm{\Omega }_{\mathrm{IB}}`$ to $`0.16\pm 0.03`$ NLOpieta . On the other hand, the complete NLO calculation of IB effects beyond the $`\pi ^0\eta \eta ^{}`$ mixing (of strong and electromagnetic origin, among which the presence of $`\mathrm{\Delta }I=5/2`$ transitions) involves a number of unknown NLO chiral couplings and is presently quite uncertain. Dimensional estimates show that IB effects may be large and affect $`\epsilon ^{}/\epsilon `$ sizeably in both directions GardnerValencia . Although a partial cancellation of the indirect ($`\mathrm{\Delta }\omega `$) and direct $`\mathrm{\Delta }\mathrm{\Omega }_{\mathrm{IB}}`$ NLO isospin breaking corrections in eq. 6 may reduce their final numerical impact on $`\epsilon ^{}/\epsilon `$ , we must await for further analyses in order to confidently assess their relevance. At present one may use $`\mathrm{\Omega }_{\mathrm{IB}}=0.10\pm 0.20`$ omegaNLOstr ; omegaNLOmodel ; omegaNLOem  as a conservative estimate of the IB effects.

The final basic ingredient for the calculation of $`\epsilon ^{}/\epsilon `$ is the evaluation of the $`K\pi \pi `$ hadronic matrix elements of the quark operators in eq. 5. A simple albeit naive approach to the problem is the Vacuum Saturation Approximation (VSA), which is based on two drastic assumptions: the factorization of the four quark operators in products of currents and densities and the saturation of the intermediate states by the vacuum state. As an example:

$`\pi ^+\pi ^{}|Q_6|K^0`$ $`=`$ $`2\pi ^{}|\overline{u}\gamma _5d|0\pi ^+|\overline{s}u|K^02\pi ^+\pi ^{}|\overline{d}d|00|\overline{s}\gamma _5d|K^0`$ (8)
$`+2\left[0|\overline{s}s|00|\overline{d}d|0\right]\pi ^+\pi ^{}|\overline{s}\gamma _5d|K^0`$
The VSA does not exhibit a consistent matching of the renormalization scale and scheme dependences of the Wilson coefficients and it carries potentially large systematic uncertainties review . On the other hand it provides useful insights on the main features of the problem.

A pictorial summary of the relative weights of the contributions of the various operators to $`\epsilon ^{}/\epsilon `$ , as obtained in the VSA, is shown in Fig. 1.

As we have already mentioned, CP violation involves loop-induced operators ($`Q_3Q_{10}`$). From Fig. 1 one clearly notices the potentially large cancellation among the strong and electroweak sectors and the leading role played by the gluonic penguin operator $`Q_6`$ and the electroweak operator $`Q_8`$. Tipical range of values for $`\epsilon ^{}/\epsilon `$ , obtained using the VSA, are shown in Fig. 2 together with the three most updated predictions available before 1999 (when the first KTeV and NA48 results became known) Munichpre ; Romapre ; ts98b . The fact that the cancellation among the strong and electroweak sectors turns out to be quite effective (at least in the VSA) warns us about the possibility that the uncertainties in the determination of the relevant hadronic matrix elements may be largely amplified in the calculation of $`\epsilon ^{}/\epsilon `$ . It is therefore important to asses carefully the approximations related to the various parts of the calculations. In particular, the analysis of the problem suggests that factorization may be highly unreliable.

## 3 Beyond Factorization

The dark gray bars in Fig. 2 represent the results of three calculations of $`\epsilon ^{}/\epsilon `$ which are representative of approaches that in principle allow us to go beyond naive factorization. They are based from left to right on the large $`N_c`$ expansion Munichpre ; BBG , on lattice regularization Romapre ; Roma , and on phenomenological modelling of low-energy QCD (the chiral quark model) chiQM ; PichEdR ; ts98a ; ts98b .

The experimental and theoretical scenarios have changed substantially after the first KTeV data and the subsequent NA48 results. Fig. 3 shows the present experimental world average for $`\epsilon ^{}/\epsilon `$ compared with the revised or new theoretical calculations that appeared during the last year. Without entering into the details of the results (for a short summary see radcor00 ) they all represent attempts to incorporate non-perturbative information into the calculation of the hadronic matrix elemnts, whether it is based on the large $`N_c`$ expansion (München Munichpost , Dortmund Dortmund , Beijing Beijing , Taipei Taipei , Valencia Valencia ), phenomenological modelling of low-energy QCD (Dubna Dubna , Trieste ts00a ; ts00b , Lund Lund ), QCD Sum Rules (Montpellier Montpellier ) or, finally, lattice regularization (Roma Romapost ).

Overall most of the theoretical calculations are consistent with a non-vanishing positive effect in the SM (with the exception of some very recent lattice results on which I will comment shortly).

At a closer look however, if we focus our attention on the central values, many of the predictions prefer the $`10^4`$ regime, whereas only a few of them stand above $`10^3`$. Is this just “noise” in the theoretical calculations? Actually, without entering the many details on which the estimates are based, most of the aforementioned difference is explained in terms of a single effect: the different size of the hadronic matrix element of the gluonic penguin $`Q_6`$ obtained in the various approaches. In turn this can be understood in terms of sizeable higher order chiral contributions (NLO in the $`1/N_c`$ expansion) to the $`I=0`$ amplitudes.

This effect was stigmatized before the latest experimental round by the work of the Trieste group ts98a ; ts98b , and appears clearly in the comparison of the leading $`1/N_c`$ and lattice results with the chiral quark model analysis in Fig. 2. The chiral quark model approach, together with the fit of the CP conserving amplitudes to phenomenologically normalize matching and model parameters, allows one to carry a complete calculation of the hadronic matrix elements beyond the leading order in the chiral expansion (including the needed local counterterms). Non-factorizable chiral contributions (not included in leading $`1/N_c`$ or lattice calculations) were shown to produce a substantial enhancement of the $`I=0`$ transitions lifting the expectation of $`\epsilon ^{}/\epsilon `$ at the $`10^3`$ level.

Since then a number of groups have attempted to improve the calculation of $`K\pi \pi `$ matrix elements in a model independent way. The present status of the art is summarized in Table 1. Due to the leading role played by $`Q_6`$ and $`Q_8`$ we may write a simplified version of eq. 6, namely

$$\frac{\epsilon ^{}}{\epsilon }13\left(\frac{\mathrm{\Lambda }_{\overline{MS}}^{(4)}}{340\mathrm{MeV}}\right)\text{Im}\lambda _t\left[\frac{110\mathrm{MeV}}{m_s(2\mathrm{GeV})}\right]^2\left[B_6(1\mathrm{\Omega }_{\mathrm{IB}})0.4B_8^{(2)}\right],$$
(9)
which in spite of the fact that it should “not be used for any serious analysis” Munichpost  gives an effective and practical way to test and compare different calculations.

The B-factors $`B_iQ_i/Q_i_{\mathrm{VSA}}`$ represent a convenient parametrization of the hadronic matrix elements, albeit tricky, in that their values are in general scale and renormalization-scheme dependent, and a spurious dependence on the quark masses is introduced in the result whenever quark densities are involved. The latter is the case for the $`Q_6`$ and $`Q_8`$ penguins. As a consequence the VSA normalization may vary from author to author thus introducing systematic ambiguities. By defining the VSA matrix elements at the scale $`\mu =2`$ GeV we obtain review

$$\begin{array}{ccc}(\pi \pi )_{I=2}|Q_8|K^0_{\mathrm{VSA}}=\sqrt{6}fm_K^4(m_s+m_d)^21.1\mathrm{GeV}^3,\hfill & & \\ \pi \pi |Q_6|K^0_{\mathrm{VSA}}/(\pi \pi )_{I=2}|Q_8|K^0_{\mathrm{VSA}}=2\sqrt{2}(f_Kf_\pi )/f_\pi 0.63,\hfill & & \end{array}$$
(10)
where I have used $`(m_s+m_d)(2\mathrm{GeV})=110\mathrm{MeV}`$ and the chiral value $`f=86`$ MeV for the octet decay constant.

It is known that $`B_6`$ and $`B_8^{(2)}`$ are perturbatively very weakly dependent on the renormalization scale review . Therefore it makes sense to compare the $`B`$’s obtained in different approaches, where the matrix elements $`Q_i`$ are computed at different scales. In Table 1 I collected some recent results for the relevant penguin matrix elements coming from various approaches, having care of normalizing the results in a homogeneous way (as far as detailed information on definitions and renormalization schemes was available).

As a guiding information, taking $`\text{Im}\lambda _t=1.3\times 10^4`$, the present experimental central value of $`\epsilon ^{}/\epsilon `$ is reproduced by $`B_6(1\mathrm{\Omega }_{\mathrm{IB}})0.4B_8^{(2)}1`$.

The most important fact is the first evidence of a signal in lattice calculations of $`\pi \pi |Q_6|K`$, obtained by the CP-PACS CP-PACS  and RBC RBC  collaborations. Both groups use the Domain Wall Fermion approach which allows to control the chiral symmetry on the lattice as a volume effect in a fifth dimension. This approach softens in principle the problem of large power subtractions which affects the lattice extraction of $`I=0`$ amplitudes (penguin contractions). Still the $`\pi |Q_i|K`$ transition is directly computed on the lattice and chiral perturbation is needed to extrapolate to the physical amplitude. The two groups obtain comparable values of the $`Q_6`$ and $`Q_8`$ matrix elements leading both to a negative $`\epsilon ^{}/\epsilon `$ . On the other hand the calculations are at an early stage and do not fully include higher order chiral contributions which may be responsible for the dynamical enhancement of the gluonic penguin, as large $`N_c`$ approaches beyond LO and the Chiral Quark Model suggest.

Important results have been obtained using data on spectral functions in connection with QCD sum rules and dispersive relations in the attempt to obtain model independent information on the relevant matrix elements. These approaches have produced as of today calculations of $`Q_8`$ (in the chiral limit) which are subtantially larger than the factorization (and lattice) results. In addition there is still substantial disagreement among different methods, in particular between the result of ref. PerisEdR , which is based on large $`N_c`$ and Lowest Meson Dominance (LMD), and that of ref. BGP  obtained via dispersive relations and data.

Calculations which sofar have gone beyond LO in the large $`N_c`$ expansion, based on chiral perturbation theory and/or models of low-energy QCD, have shown the crucial role of higher order non-factorizable corrections in the enhancement of the $`I=0`$ matrix elements ts98b ; Dortmund ; Valencia ; Lund . Chiral loop corrections drive the final value of $`\epsilon ^{}/\epsilon `$ in the ballpark of the present data. However the calculation of higher order chiral corrections cannot be fully accomplished in a model independent way due to the many unknown NLO local couplings. In the chiral quark model approach all needed local coupling are computed in terms of quark masses, meson decay constants and a few non-perturbative parameters as quark and gluon condensates. The latter are determined self-consistently in a phenomenological way via the fit of the CP conserving $`K\pi \pi `$ amplitudes, thus encoding the $`\mathrm{\Delta }I=1/2`$ rule in the calculation ts98a ; ts00a . The analysis shows that the role of local counterterms is subleading to the chiral logs when using the Modified Minimal Subtraction as opposed to the commonly used Gasser-Leutwyler prescription.

Among non-factorizable chiral corrections FSI play a leading role. As a matter of fact, one should in general expect an enhancement of $`\epsilon ^{}/\epsilon `$ with respect to the naive VSA due to FSI. As Fermi first argued Fermi , in potential scattering the isospin $`I=0`$ two-body states feel an attractive interaction, of a sign opposite to that of the $`I=2`$ components thus affecting the size of the corresponding amplitudes. This feature is at the root of the enhancement of the $`I=0`$ amplitude over the $`I=2`$ one and of the corresponding enhancement of $`\epsilon ^{}/\epsilon `$ beyond factorization. An attempt to resum these effects in a model independent way has been worked out by the authors of ref. Valencia , using a dispersive approach a la Omnès-Mushkelishvili Omnes . Their analysis shows that resummation does not substantially modify the one-loop perturbative result and, as it appears from Table 1, a 50% enhancement of the gluonic penguin matrix element is found over the factorized result. However, the analysis suffers from a sistematic uncertainty due to the value of the factorized amplitude at the chosen subtraction point which is identified with the large $`N_c`$ result FSIcritics . Even when the authors in the most recent work match the dispersive resummation with the perturbative one-loop calculation, thus including $`1/N_c`$ effects, again a systematic uncertainty remains hidden in the unknown polinomial parts of the local chiral counterterms. In this respect it is important to stress that gauging the scheme dependence of the result by varying the renormalization scale in the chiral logs cannot account for the arbitrariness related to the chiral-loop subtraction (a model-independent complete calculation is missing).

Finally, it has been recently emphasized dim8  that cut-off based approaches should pay attention to higher-dimension operators which become relevant for matching scales below 2 GeV and may represent one of the largest sources of uncertainty in present calculations. The results of refs. CDGM ; BGP  include these effects. The calculations based on dimensional regularization may be safe if phenomenological input is used in order to encode in the relevant hadronic matrix elements the physics at all scales (this is done in the Trieste approach).

In summary, while some model dependent calculations suggest no conflict between theory and experiment for $`\epsilon ^{}/\epsilon `$ , a precise and ”pristine” prediction of the observable is still quite ahead of us.

## 4 Outlook and Conclusions

Higher-order chiral corrections are taking the stage of $`K\pi \pi `$ physics. They are needed in order to asses the size of crucial parameters (as $`\mathrm{\Omega }_{\mathrm{IB}}`$) and the effect of non-factorizable contributions in the penguin matrix elements.

Lattice, as a regularization of QCD, is the first-principle approach to the problem. However, lattice calculations still heavily depend on chiral perturbation theory Golterman . Presently, very promising developments are being undertaken to circumvemt the technical and conceptual shortcomings related to the calculation of weak matrix elements (for a recent survey see ref. Sachrajda ). Among those are the Domain Wall Fermion approach DWF  which allows us to decouple the chiral symmetry from the continuum limit, and the very interesting observation that the Maiani-Testa theorem MaianiTesta  can be overcomed using the fact that lattice calculations are performed in finite volume LellouchLuscher , thus allowing for the direct calculation of the physical $`K\pi \pi `$ amplitude on the lattice. All these developments need a tremendous effort in machine power and in devising faster algorithms. Preliminary results for lattice calculations of both $`\epsilon ^{}/\epsilon `$ and the $`\mathrm{\Delta }I=1/2`$ selection rule are already available and others are currently under way SPQcdR .

In the meantime analytical and semi-phenomenological approaches have been crucially helpful in driving the attention of the community on some systematic short-comings of ”first-principle” calculations. The amount of theoretical work triggered by the NA48 and KTeV data promises rewarding and perhaps exciting results within the next few years.

I wish to thank the organizers of Heavy Flavours 9 for the excellent organization of the meeting at a time when world events played the lead in everyone’s mind.


text:
# Universal Equation for Efimov States

thanks:  Present address: Helmholtz-Institut für Strahlen- und Kernphysik (Abt. Theorie), Universität Bonn, 53115 Bonn, Germany (November 22, 2002)

## Abstract

Efimov states are a sequence of shallow 3-body bound states that arise when the 2-body scattering length is large. Efimov showed that the binding energies of these states can be calculated in terms of the scattering length and a 3-body parameter by solving a transcendental equation involving a universal function of one variable. We calculate this universal function using effective field theory and use it to describe the three-body system of 4He atoms. We also extend Efimov’s theory to include the effects of deep 2-body bound states, which give widths to the Efimov states.

Efimov states, universality, 4He, effective field theory

The interactions of nonrelativistic particles (such as atoms) with short-range interactions at extremely low energies are determined primarily by their $`S`$-wave scattering length $`a`$. If $`|a|`$ is much larger than the characteristic range $`l`$ of their interaction, low-energy atoms exhibit universal properties that are insensitive to the details of the interaction potential. In the 2-body sector, the universal properties are simple and familiar. The differential cross section for two identical bosons with relative wave number $`k1/l`$ and mass $`m`$ is $`d\sigma /d\mathrm{\Omega }=4a^2/(1+k^2a^2)`$. If $`a>0`$, there is also a shallow 2-body bound state (the dimer) with binding energy $`B_2=\mathrm{}^2/ma^2`$. In the 3-body sector, there are also universal properties that were first deduced by Efimov Efi71 . The most remarkable is the existence of a sequence of 3-body bound states with binding energies geometrically spaced in the interval between $`\mathrm{}^2/ma^2`$ and $`\mathrm{}^2/ml^2`$. The number of these “Efimov states” is roughly $`\mathrm{ln}(|a|/l)/\pi `$ if $`|a|`$ is large enough. In the limit $`|a|\mathrm{}`$, there is an accumulation of infinitely many 3-body bound states at threshold (the “Efimov effect”). The knowledge of the Efimov binding energies is essential for understanding the energy-dependence of low-energy 3-body observables. For example, Efimov states can have dramatic effects on atom-dimer scattering if $`a>0`$ Efi71 ; BHK99  and on 3-body recombination if $`a<0`$ EGB99 ; BBH00 .

A large 2-body scattering length can be obtained by fine-tuning a parameter in the interatomic potential to bring a real or virtual 2-body bound state close to the 2-atom threshold. The fine tuning can be provided accidentally by nature. An example is the 4He atom, whose scattering length $`a=104`$ Å Gri00  is much larger than the effective range $`l7`$ Å. Another example is the 2-nucleon system in the $`{}_{}{}^{3}S_{1}^{}`$ channel, for which the deuteron is the shallow bound state. This system provided the original motivation for Efimov’s investigations Efi71 . In the case of atoms, the fine tuning can also be obtained by tuning an external electric field NFJ99  or by tuning an external magnetic field to the neighbourhood of a Feshbach resonance Festh . Such resonances have, e.g., been observed for 23Na and 85Rb atoms Fesex  and are used to tune the interactions in Bose-Einstein condensates. An important difference from He is that the interatomic potentials for Na and Rb support many deep 2-body bound states.

Efimov derived some powerful constraints on low-energy 3-body observables for systems with large scattering length Efi71 . They follow from the approximate scale-invariance at length scales $`R`$ in the region $`lR|a|`$ and the conservation of probability. He introduced polar variables $`H`$ and $`\xi `$ in the plane whose axes are $`1/a`$ and the energy variable $`\mathrm{sgn}(E)|mE|^{1/2}/\mathrm{}`$, and showed that low-energy 3-body observables are determined by a few universal functions of the angle $`\xi `$. In particular, the binding energies of the Efimov states are solutions to a transcendental equation involving a single universal function of $`\xi `$ Efi71 . In this paper, we calculate this universal function for the case of 3 identical bosons and apply Efimov’s equation to the 4He trimer. We also extend Efimov theory to atoms with deep 2-body bound states.

The existence of Efimov states can easily be understood in terms of the equation for the radial wave function $`f(R)`$ in the adiabatic hyperspherical representation of the 3-body problem Jen93 ; Nie01 . The hyperspherical radius for 3 identical atoms with coordinates $`𝐫_1`$, $`𝐫_2`$, and $`𝐫_3`$ is $`R^2=(r_{12}^2+r_{13}^2+r_{23}^2)/3`$, where $`r_{ij}=|𝐫_i𝐫_j|`$. If $`|a|l`$, the radial equation for 3 atoms with total angular momentum zero reduces in the region $`lR|a|`$ to

$$\frac{\mathrm{}^2}{2m}\left[\frac{^2}{R^2}+\frac{s_0^2+1/4}{R^2}\right]f(R)=Ef(R),$$
(1)
where $`s_01.00624`$. This looks like the Schrödinger equation for a particle in a one-dimensional, scale-invariant $`1/R^2`$ potential. If we impose a boundary condition on $`f(R)`$ at short-distances of order $`l`$, the radial equation (1) has solutions at discrete negative values of the eigenvalue $`E=B_3`$, with $`B_3`$ ranging from order $`\mathrm{}^2/ml^2`$ to order $`\mathrm{}^2/ma^2`$. The corresponding eigenstates are called Efimov states. As $`|a|\mathrm{}`$, their spectrum approaches the simple law $`B_3^{(n)}515^n\mathrm{}^2/ma^2`$.

Efimov’s constraints can be derived by constructing a solution to Eq. (1) that is valid in the region $`lR|a|`$. In the case of a bound state with energy $`E=B_3`$, the radial variable is $`H^2=mB_3/\mathrm{}^2+1/a^2`$ and the angular variable $`\xi `$ is

$`\xi `$ $`=`$ $`\mathrm{arctan}(a\sqrt{mB_3}/\mathrm{})\pi \theta (a),`$ (2)
where $`\theta (x)`$ is the unit step function. Since we are interested in low energies $`|E|\mathrm{}^2/ma^2`$, the energy eigenvalue in (1) can be neglected. The most general solution therefore has the form Efi71

$$f(R)=\sqrt{HR}\left[Ae^{is_0\mathrm{ln}(HR)}+Be^{is_0\mathrm{ln}(HR)}\right],$$
(3)
which is the sum of outgoing and incoming hyperspherical waves. The dimensionless coefficients $`A`$ and $`B`$ can depend on $`\xi `$. At shorter distances $`Rl`$ and longer distances $`R|a|`$, the wavefunction becomes very complicated. Fortunately, we can avoid solving for it by using simple considerations based on unitarity Efi71 .

We first consider the short-distance region. Efimov assumed implicitly that there are no deep 2-body bound states with binding energies $`B\mathrm{}^2/ma^2`$. Thus the 2-body potential supports no bound states at all if $`a<0`$ and only the dimer with binding energy $`B_2=\mathrm{}^2/ma^2`$ if $`a>0`$. We will address the complication of deep 2-body bound states later. The probability in the incoming wave must then be totally reflected by the potential at short distances, so we can set $`B=Ae^{i\theta }`$. The phase $`\theta `$ can be specified by giving the logarithmic derivative $`R_0f^{}(R_0)/f(R_0)`$ at any point $`lR_0|a|`$. The resulting expression for $`\theta `$ has a simple dependence on $`H`$:

$`\theta /2`$ $`=`$ $`s_0\mathrm{ln}(H/c\mathrm{\Lambda }_{}).`$ (4)
The denominator $`c\mathrm{\Lambda }_{}`$ is a complicated function of $`R_0`$ and $`R_0f^{}(R_0)/f(R_0)`$. It differs by an unknown constant $`c`$ from the 3-body parameter $`\mathrm{\Lambda }_{}`$ introduced in Ref. BHK99 .

We next consider large distances $`R|a|`$. In general, an outgoing hyperspherical wave incident on the $`R|a|`$ region can either be reflected or else transmitted to $`R\mathrm{}`$ as a 3-atom or atom-dimer scattering state. The reflection and transmission amplitudes are described by a dimensionless unitary $`3\times 3`$ matrix that is a function of $`\xi `$ only. For $`\pi <\xi <\pi /4`$, scattering states with $`R\mathrm{}`$ are kinematically not allowed. The probability is therefore totally reflected, so we must have $`B=Ae^{i\mathrm{\Delta }(\xi )}`$, where the phase $`\mathrm{\Delta }`$ depends on the angle $`\xi `$. Compatibility with the constraint from short distances requires $`\theta =\mathrm{\Delta }(\xi )\mathrm{mod}\mathrm{\hspace{0.33em}2}\pi `$. Using Eq. (4) for $`\theta `$ and inserting the expression for $`H`$, we obtain Efimov’s equation Efi71

$$B_3+\frac{\mathrm{}^2}{ma^2}=\frac{\mathrm{}^2\mathrm{\Lambda }_{}^2}{m}e^{2\pi n/s_0}\mathrm{exp}\left[\mathrm{\Delta }\left(\xi \right)/s_0\right],$$
(5)
where $`\xi `$ is given by (2) and the constant $`c`$ was absorbed into $`\mathrm{\Delta }(\xi )`$. Note that we measure $`B_3`$ from the 3-atom threshold and $`\mathrm{\Lambda }_{}`$ is only defined up to factors of $`\mathrm{exp}[\pi /s_0]`$. Once the universal function $`\mathrm{\Delta }(\xi )`$ is known, the Efimov binding energies $`B_3`$ can be calculated by solving Eq. (5) for different integers $`n`$. This equation has an exact discrete scaling symmetry: if there is an Efimov state with binding energy $`B_3`$ for the parameters $`a`$ and $`\mathrm{\Lambda }_{}`$, then there is also an Efimov state with binding energy $`\lambda ^2B_3`$ for the parameters $`\lambda ^1a`$ and $`\mathrm{\Lambda }_{}`$ if $`\lambda =\mathrm{exp}[n^{}\pi /s_0]`$ with $`n^{}`$ an integer.

The universal function $`\mathrm{\Delta }(\xi )`$ could be determined by solving the 3-body equation for the Efimov binding energies in various potentials whose scattering lengths are so large that effective range corrections are negligible. It can be calculated more easily by using the effective field theory (EFT) of Ref. BHK99  in which the effective range can be set to zero. In Ref. BHK99 , the dependence of the binding energy on $`a`$ and $`\mathrm{\Lambda }_{}`$ was calculated for the shallowest Efimov state and $`a>0`$. In order to extract the universal function $`\mathrm{\Delta }(\xi )`$, we have calculated the binding energies of the three lowest Efimov states for both signs of $`a`$.

In Fig. 1, we plot $`(mB_3/\mathrm{}^2\mathrm{\Lambda }_{}^2)^{1/4}`$ as a function of $`\mathrm{sgn}(a)(\mathrm{\Lambda }_{}|a|)^{1/2}`$ for these three branches of Efimov states. The binding energies for deeper Efimov states and for shallower states near $`(\mathrm{\Lambda }_{}|a|)^{1/2}=0`$ can be obtained from the discrete scaling symmetry. A given 2-body potential is characterized by values of $`a`$ and $`\mathrm{\Lambda }_{}`$ and corresponds to a vertical line in Fig. 1, such as the dashed line shown. The intersections of this line with the binding energy curves correspond to the infinitely many Efimov states. Those states with $`B_3>\mathrm{}^2/ml^2`$ are unphysical. For $`B_3\mathrm{}`$, the angle $`\xi `$ goes to $`\pi /2`$. The ratio of the binding energies of successive Efimov states therefore approaches $`\mathrm{exp}[2\pi /s_0]515`$. However, for the shallowest Efimov states, this ratio exhibits significant deviations from the asymptotic value. If $`a>0`$, there is an Efimov state at the atom-dimer threshold $`B_3=\mathrm{}^2/ma^2`$ when $`s_0\mathrm{ln}(a\mathrm{\Lambda }_{})=1.444\mathrm{mod}\pi `$. The sequence of binding energies $`B_3`$ in units of $`\mathrm{}^2/ma^2`$ is 1, 6.8, $`1.4\times 10^3`$,$`\mathrm{}`$. Consequently, the ratio of $`B_3`$ for the two shallowest Efimov states can range from 6.8 to 210. If $`a<0`$, there is an Efimov state at the 3-atom threshold $`B_3=0`$ when $`s_0\mathrm{ln}(a\mathrm{\Lambda }_{})=1.378\mathrm{mod}\pi `$. The sequence of binding energies $`B_3`$ in units of $`\mathrm{}^2/ma^2`$ is 0, $`1.1\times 10^3`$, $`6.0\times 10^5`$,$`\mathrm{}`$. Thus, the ratio of $`B_3`$ for the two shallowest Efimov states can range from $`\mathrm{}`$ to 550.

Using Eq. (5) with $`n=0`$, we have extracted the universal function $`\mathrm{\Delta }(\xi )`$ from the data for the middle branch in Fig. 1. The extracted values of $`\mathrm{\Delta }(\xi )`$ are given in Table 1.

An analytic expression for $`\mathrm{\Delta }(\xi )`$ is not known, but we can obtain parametrizations in various regions of $`\xi `$ by fitting the data:

$`\xi [\frac{3\pi }{8},\frac{\pi }{4}]:\mathrm{\Delta }=3.10x^29.63x2.18,`$ (6)
$`\xi [\frac{5\pi }{8},\frac{3\pi }{8}]:\mathrm{\Delta }=1.17y^3+1.97y^2+2.12y8.22,`$ (7)
$`\xi [\pi ,\frac{5\pi }{8}]:\mathrm{\Delta }=0.25z^2+0.28z9.11,`$ (8)
where $`x=(\pi /4\xi )^{1/2}`$, $`y=\pi /2+\xi `$, and $`z=(\pi +\xi )^2\mathrm{exp}[1/(\pi +\xi )^2]`$. These parametrizations deviate from the numerical results in Table 1 by less than 0.013. The discontinuity at $`\xi =\frac{3\pi }{8}`$ and $`\xi =\frac{5\pi }{8}`$ is less than 0.016. This accuracy is sufficient for most practical calculations using Eq. (5).

Universality can be exploited to greatly reduce the calculational effort required to predict 3-body observables for atoms with large $`a`$. The observables can be calculated once and for all as functions of $`a`$ and $`\mathrm{\Lambda }_{}`$ either by using EFT or by solving the Schrödinger or Faddeev equations with various methods. The binding energies obtained by solving Efimov’s equation (5) are shown in Fig. 2.

Simple expressions can be given for other observables, such as the $`S`$-wave atom-dimer scattering length:

$$a_{12}=a\left(1.462.15\mathrm{tan}[s_0\mathrm{ln}(a\mathrm{\Lambda }_{})+0.09]\right),$$
(9)
as well as the phase shifts and the rate constant for 3-body recombination at threshold BrH02 . Given $`a`$ and a measured or calculated value of $`B_3`$ for one Efimov state as input, one can read off $`\mathrm{\Lambda }_{}`$ from Fig. 2. Predictions for other 3-body observables, such as the atom-dimer scattering length in (9), are then immediate.

One of the most promising systems for observing Efimov states is 4He atoms. The 4He trimer has been observed STo96 , but no quantitative experimental information about its binding energy is available to date. The binding energy has been calculated accurately for various model potentials. They indicate that there are two trimers, a ground state with binding energy $`B_3^{(0)}`$ and an excited state with binding energy $`B_3^{(1)}`$. The most accurate calculations have been obtained by solving the Faddeev equations in the hyperspherical representation NFJ98 , in configuration space RoY00 , and with hard-core boundary conditions MSSK01 . These methods all give consistent results. The results of Ref. MSSK01  for $`B_2`$, $`B_3^{(0)}`$, $`B_3^{(1)}`$, and the atom-dimer scattering length $`a_{12}`$ for four different model potentials are given in Table 2.

The discrepancies between $`B_2`$ and the large-$`a`$ prediction $`\mathrm{}^2/ma^2`$ are about 6-8% MSSK01 . They can be attributed to effective range corrections and provide estimates of the error associated with the large-$`a`$ approximation.

We proceed to illustrate the power of universality by applying it to the 4He trimer. For the scattering length $`a`$, we take the value $`a_B\mathrm{}/\sqrt{mB_2}`$ obtained from the calculated dimer binding energy. We determine $`\mathrm{\Lambda }_{}`$ by demanding that $`B_3^{(1)}`$ satisfy (5) with $`n=1`$. Solving Eq. (5) with $`n=2`$, we obtain the predictions for $`B_3^{(0)}`$ in the second-to-last column of Table 2. The predictions are only 1-4% higher than the calculated values, which is within the expected error for the large-$`a`$ approximation. This demonstrates that the ground state of the 4He trimer can be described by Efimov’s equation (5). If we use the calculated values of $`a`$ as input instead of $`B_2`$, the predicted values of $`B_3^{(0)}`$ are larger than the calculated values by 11-21%.

We can use the value of $`a_B\mathrm{\Lambda }_{}`$ determined from the excited state of the trimer to predict the atom-dimer scattering length $`a_{12}`$ and compare with the calculated values in the fourth column of Table 2. The predictions for the four model potentials are given in the last column of Table 2. They are smaller than the calculated values by about 13%. If the calculated value of $`a`$ is used as input they are smaller by about 28%. It should be possible to account for these differences quantitatively by taking into account higher order corrections BrH02 .

Efimov implicitly assumed in his derivation of Eq. (5) that there were no deep 2-body bound states Efi71 . If such states are present, the Efimov states become resonances that can decay into a deep 2-body bound state and a recoiling atom. Thus their energies are given by complex numbers $`E=B_3i\mathrm{\Gamma }_3/2`$. If a potential supports many 2-body bound states, the direct calculation of the widths $`\mathrm{\Gamma }_3`$ by solving the Schrödinger equation is very difficult NiE00 . However, one can show that the cumulative effect of all deep 2-body bound states on low-energy 3-body observables can be taken into account by including one additional low-energy parameter $`\eta `$. Three-body recombination into a deep 2-body bound state with binding energy $`B\mathrm{}^2/ml^2`$ can only take place if $`Rl`$. It is obvious that the atoms that form the bound state must approach to within a distance of order $`l`$, since the size of the bound state is of order $`l`$. However, the third atom must also approach the pair to within a distance of order $`l`$, because it must recoil with momentum $`\sqrt{4mB/3}\mathrm{}/l`$, and the necessary momentum kick can be delivered only if $`Rl`$. The atoms can approach such short distances only by following the lowest continuum adiabatic hyperspherical potential, because it is attractive in the region $`lR|a|`$, while all other potentials are repulsive. Thus all pathways to final states including a deep 2-body bound state must flow through this lowest continuum adiabatic potential. Note the lowest continuum state can either be a 3-atom state or an atom-dimer state. See, e.g., Fig. 5 in Ref. Nie01  for an explicit calculation of these potentials. The cumulative effects of 3-body recombination into deep 2-body bound states can therefore be described by the reflection probability $`e^\eta `$ for hyperspherical waves entering the region $`Rl`$ of this potential. Up to corrections suppressed by $`l/|a|`$, the low-energy 3-body observables are all determined by $`a`$, $`\mathrm{\Lambda }_{}`$, and $`\eta `$.

We proceed to generalize Efimov’s equation (5) to the case in which there are deep 2-body bound states. Again, we combine the analytic solution (3) to the radial equation in the region $`lR|a|`$ with simple probability considerations at $`R|a|`$ and $`Rl`$. Since the existence of deep 2-body bound states plays no role in the unitarity constraint at $`R|a|`$, we still have $`B=Ae^{i\mathrm{\Delta }(\xi )}`$ with the same universal function $`\mathrm{\Delta }(\xi )`$ given in (6)-(8). In the unitarity constraint at $`Rl`$, we need to take into account that only a fraction of the probability that flows to short distances is reflected back to long distances through the lowest adiabatic hyperspherical potential associated with the 3-atom or atom-dimer continuum. Denoting the reflection probability by $`e^\eta `$, the coefficient $`B`$ in (3) can be written $`B=Ae^{i\theta +\eta /2}`$. The compatibility condition $`\theta i\eta /2=\mathrm{\Delta }(\xi )\mathrm{mod}\mathrm{\hspace{0.33em}2}\pi `$ then becomes

$$B_3+\frac{i}{2}\mathrm{\Gamma }_3+\frac{\mathrm{}^2}{ma^2}=\frac{\mathrm{}^2\mathrm{\Lambda }_{}^2}{m}e^{2\pi n/s_0}e^{(\mathrm{\Delta }(\xi )+i\eta /2)/s_0},$$
(10)
where $`\xi `$ is defined by (2) with $`B_3B_3+i\mathrm{\Gamma }_3/2`$. To solve this equation, we need the analytic continuation of $`\mathrm{\Delta }(\xi )`$ to complex values of $`\xi `$. The parametrization (6)-(8) for $`\mathrm{\Delta }(\xi )`$ should be accurate for complex $`\xi `$ with sufficiently small imaginary parts, except perhaps near $`\xi =\pi `$ where it has an essential singularity. If $`B_3`$ and $`\mathrm{\Gamma }_3`$ for one Efimov state are known, they can be used to determine $`\mathrm{\Lambda }_{}`$ and $`\eta `$. The remaining spectrum of Efimov states and their widths can then be calculated by solving Eq. (10). For infinitesimal $`\eta `$, the widths approach $`\mathrm{\Gamma }_3(\eta /s_0)(B_3+\mathrm{}^2/ma^2)`$. The widths of the deeper Efimov states increase geometrically just like the binding energies, as has been observed in numerical calculations NiE00 . If $`\eta `$ is so large that $`B_3\mathrm{\Gamma }_3`$, the Efimov states cease to exist in any meaningful sense.

We have calculated the universal function $`\mathrm{\Delta }(\xi )`$ that appears in Efimov’s equation (5) for the binding energies $`B_3`$ of Efimov states. This equation can be used as an operational definition of the 3-body parameter $`\mathrm{\Lambda }_{}`$ introduced in Ref. BHK99 . If $`B_3`$ for one Efimov state is known, it can be used to determine $`\mathrm{\Lambda }_{}`$, and then universality predicts all low-energy 3-body observables as functions of $`a`$ and $`\mathrm{\Lambda }_{}`$. In Eq. (10), we have generalized Efimov’s equation to permit deep 2-body bound states. The generalization involves an additional inelasticity parameter $`\eta `$, but the spectrum is determined by the same universal function $`\mathrm{\Delta }(\xi )`$.

This research was supported by DOE Grant No. DE-FG02-91-ER4069 and NSF Grant No. PHY-0098645.


text:
# Phenomenological study of charm photoproduction at HERA

## 1 Introduction

In the present work we discuss phenomenological applications of the formalism presented in ref. , relevant to the production of heavy quarks in photon-hadron collisions. All technical details and theoretical motivations for the formalism are discussed there. In the formalism of ref.  one obtains a cross section which is accurate to $`𝒪(\alpha _{\mathrm{em}}\alpha _\mathrm{s}^2)`$ (that is, to the next-to-leading order level) and which, in the large $`p_\mathrm{T}`$ region, includes consistently all enhanced logarithmic terms of the form $`\alpha _{\mathrm{em}}\alpha _\mathrm{s}(\alpha _\mathrm{s}\mathrm{log}p_\mathrm{T}/m)^i`$ and $`\alpha _{\mathrm{em}}\alpha _\mathrm{s}^2(\alpha _\mathrm{s}\mathrm{log}p_\mathrm{T}/m)^i`$. Terms of order $`𝒪(\alpha _{\mathrm{em}}\alpha _\mathrm{s}^3)`$ with no logarithmic enhancement and terms of order $`\alpha _{\mathrm{em}}\alpha _\mathrm{s}^3(\alpha _\mathrm{s}\mathrm{log}p_\mathrm{T}/m)^i`$ are not included. Thus this approach improves both the fixed order (improperly called massive) and the resummed (improperly called massless) calculations.

Our starting point is the “matched” formula

$$\text{FONLL}=\text{FO}+\left(\text{RS}\text{FOM0}\right)\times G(m,p_\mathrm{T}),$$
(1)
where FONLL stands for fixed-order plus next-to-leading logarithms, and

* FO is the fixed order, $`𝒪(\alpha _{\mathrm{em}}\alpha _\mathrm{s}^2)`$ result;
* RS is the resummed result, which includes all terms of the form $`\alpha _{\mathrm{em}}\alpha _\mathrm{s}(\alpha _\mathrm{s}\mathrm{log}p_\mathrm{T}/m)^i`$ and $`\alpha _{\mathrm{em}}\alpha _\mathrm{s}^2(\alpha _\mathrm{s}\mathrm{log}p_\mathrm{T}/m)^i`$, and neglects all terms suppressed by powers of the heavy quark mass $`m`$;
* FOM0 is the massless limit of FO, in the sense that all terms suppressed by powers of $`m`$ are dropped, while logarithms of the mass are retained; thus FOM0 is the truncation of RS to order $`\alpha _{\mathrm{em}}\alpha _\mathrm{s}^2`$;
* $`G(m,p_\mathrm{T})`$ is an arbitrary dumping function, that must be regular in $`p_\mathrm{T}`$, and must approach 1 up to terms suppressed by powers of $`m/p_\mathrm{T}`$ at large $`p_\mathrm{T}`$; our standard choice is

$$G(m,p_\mathrm{T})=\frac{p_\mathrm{T}^2}{p_\mathrm{T}^2+c^2m^2},$$
(2)
which for $`c0`$ dumps the $`\text{RS}\text{FOM0}`$ term in the small-$`p_\mathrm{T}`$ region.

We remind the reader that in the description given above the hadronic photon contribution is also included; the counting of the orders in $`\alpha _{\mathrm{em}}`$ and $`\alpha _\mathrm{s}`$ remains the same, provided one counts the parton density functions of the photon as carrying a $`\alpha _{\mathrm{em}}/\alpha _\mathrm{s}`$ factor. A more detailed discussion can be found in ref. .

This paper is organized as follows. In sect. 2 we specify our choices for the input parameters and the parton density functions (PDF’s). Furthermore, the specific cuts that the Zeus and H1 collaborations use in their analyses are reported in this section. All the subsequent results will be based upon these choices.

In order to present meaningful phenomenological results, we need to consider a non-perturbative fragmentation function (NPFF) to account for the hadronization of the charm quark into a charmed hadron. The NPFF is the source of several uncertainties affecting the theoretical predictions. First of all, various forms of NPFF have been proposed in the literature (see ref.  and references therein for an updated review). Furthermore, while the application of NPFF is unambiguous in the large $`p_\mathrm{T}`$ region, at moderate and small $`p_\mathrm{T}`$ there are ambiguities, depending upon the definition of the $`z`$ variable (i.e., whether it refers to energy, momentum or “plus” component), and upon the choice of the frame in which it is defined. The impact of these ambiguities will be assessed in sect. 3.

Phenomenological results have also uncertainties depending upon the accuracy of the perturbative calculations (which we estimate by varying the renormalisation and factorisation scales), the allowed range in the value of the heavy quark mass and $`\mathrm{\Lambda }_{\mathrm{QCD}}`$, and the choice of the PDF’s. These uncertainties will be discussed in sect. 4.

The H1  and the Zeus  collaborations have studied photoproduction of $`D^{}`$ mesons. The two experiments have different acceptances and apply different cuts to their data samples. In sect. 5 we compute the distributions presented in refs. , implementing the same cuts that are used there.

In sect. 6 we give our conclusions.

## 2 Input parameters

All the results presented in this paper are relevant to photoproduction experiments where the photon arises from collinear radiation off an incoming electron. The radiation spectrum depends upon the cuts that are applied to define the photoproduction regime, which are typically $`Q^2`$ cuts. Further kinematic cuts are applied either to the the variable $`y`$ (the ratio of the photon energy over the incoming electron energy) or, which is equivalent, to the variable $`W`$ (the photon-proton center-of-mass energy). The cuts used in refs.  are summarised in table 1.

The H1 collaboration uses two different electron taggers, denoted as ETAG33 and ETAG44, and has thus two samples of photoproduction data. The electron and proton energies are fixed to $`27.5`$ and $`820`$ GeV respectively. The $`y`$ range quoted for Zeus corresponds to $`130<W<280`$ GeV.

We compute the photon spectrum in the Weizsäcker-Williams approximation, in the improved form of ref. , which includes non-logarithmic terms enhanced at small $`y`$.

The central values and ranges of the physical parameters that we adopt in the present paper are summarised in table 2. The range in $`\mathrm{\Lambda }_{\mathrm{QCD}}^{(5)}`$ reported in the table corresponds to the range $`\alpha _\mathrm{s}(M_\mathrm{Z})=0.118\pm 0.002`$. We use the two-loop expression of $`\alpha _\mathrm{s}`$ throughout this work.

The factorization and renormalization scales ($`\mu _\mathrm{F}`$ and $`\mu _\mathrm{R}`$) are set equal to $`\mu _0/2`$, $`\mu _0`$ and $`2\mu _0`$, where $`\mu _0=\sqrt{p_\mathrm{T}^2+m^2}`$ is our default choice.

In table 3 we report the kinematic ranges of the cross section measurements used in the Zeus and H1 analyses. Notice that the $`D^{}`$ kinematics are described by $`p_\mathrm{T},\eta `$ (pseudorapidity) in the Zeus analysis, and by $`p_\mathrm{T},y`$ (rapidity) in the H1 analysis. Unfortunately, the same notation ($`y`$) is traditionally used to indicate the Weizsäcker-Williams variable, and the meson rapidity. However, no confusion should possibly arise. In what follows, $`y`$ will always indicate the meson rapidity.

## 3 Non-perturbative fragmentation function

We begin by discussing the ambiguities related to the implementation of the non-perturbative fragmentation in our calculation. A NPFF has to be convoluted with our full FONLL result. The convolution is unambiguously defined only at large transverse momenta (i.e. $`p_\mathrm{T}m`$). One writes

$$\frac{d^3\sigma _H(k)}{d^3k}=d^3\widehat{k}𝑑zD_{\mathrm{NP}}(z)\frac{d^3\sigma _Q(\widehat{k})}{d^3\widehat{k}}\delta ^3\left(\stackrel{}{k}z\stackrel{}{\widehat{k}}\right),$$
(3)
where $`\sigma _H`$ is the cross section for the production of the heavy-flavoured hadron111In what follows, we assume the hadron $`H`$ to be a $`D^{}`$ meson, since the H1 and Zeus data that we use in this paper are relevant to $`D^{}`$ mesons. $`H`$ with momentum $`k`$, and $`\sigma _Q`$ is the cross section for the production of the heavy quark $`Q`$ with momentum $`\widehat{k}`$. $`D_{\mathrm{NP}}(z)`$ is the NPFF. In eq. (3) we assume that fragmentation scales the 3-momentum of the heavy quark. Different prescriptions are however possible: for example, one could assume that either $`\widehat{k}^0`$ or $`\widehat{k}^0+|\widehat{k}|`$ are scaled instead of $`\stackrel{}{k}`$. In all cases, it is assumed that $`\stackrel{}{k}`$ remains parallel to $`\stackrel{}{\widehat{k}}`$, and the mass-shell condition is imposed to fix the remaining components of the meson momentum. Only in the large-$`p_\mathrm{T}`$ limit all these prescriptions coincide. Also, none of these definitions is boost invariant, since the condition of parallelism is frame dependent. Again, in the large-$`p_\mathrm{T}`$ limit boost invariance is recovered.

For the meson momentum $`k`$ we expect the condition $`k^2=m_H^2`$ to be satisfied. However, for the purpose of computing distributions, we can also use $`k^2=m^2`$, and consider the difference to be just another ambiguity in our result. The reader may find unreasonable to adopt this second choice, since physically the on-shell condition of the meson $`H`$ must apply. However, in the context of a perturbative QCD computation, setting $`k^2=m_H^2`$ may appear unreasonable from other points of view. For example, if we produce a quark with maximal $`p_\mathrm{T}`$, and fragment it with $`z=1`$, setting $`k^2=m_H^2`$ would yield a meson with the same $`p_\mathrm{T}`$ but larger energy, which violates energy conservation. Similar examples can be given at small $`p_\mathrm{T}`$ and large rapidities. If we are far from the phase space boundaries, which is always our case, both choices are sensible, and their difference can be considered as an inherent uncertainty of the procedure.

In the following, we will always use the Peterson fragmentation function

$$D_{\mathrm{NP}}(z)\frac{1}{z}\left(1\frac{1}{z}\frac{ϵ}{1z}\right)^2.$$
(4)
For the sake of the studies performed in the present section we will take it normalized to 1, i.e. we do not include the $`D^{}`$ production fraction. Although other choices of fragmentation function are possible, we only use eq. (4) because it is the most commonly adopted, and because recent phenomenological fits of the $`ϵ`$ parameter, adequate to the present work, do exist . In refs.  the reader will find studies of alternative forms.

Our default choices for the non-perturbative fragmentation are as follows:

* the $`z`$ variable is defined by the equation $`\stackrel{}{k}=z\stackrel{}{\widehat{k}}`$;
* the fragmentation frame (i.e., the frame where the equation $`\stackrel{}{k}=z\stackrel{}{\widehat{k}}`$ holds) is the HERA laboratory frame;
* the mass shell condition is $`k^2=\widehat{k}^2=m^2`$;
* the value of the Peterson parameter is $`ϵ=0.02`$.

We now study the effect of the ambiguities in the definition of the non-perturbative fragmentation on the pseudorapidity and on the transverse momentum distributions for the Zeus kinematic cuts. Our conclusions would not change had we used H1 cuts. We consider first the effect of the change of the fragmentation frame. We call scheme (A) our default choice of frame. We call scheme (B) the alternative choice of a longitudinally-boosted fragmentation frame, boosted with rapidity equal to the quark rapidity. Thus, in this frame the quark rapidity is zero. We notice that in scheme (A) the quark and the meson have the same pseudorapidity, while in scheme (B) they have the same rapidity.

In figure 1 we show the difference in the pseudorapidity and transverse momentum distributions when using schemes (A) and (B). Furthermore, the effect of varying the value of the Peterson parameter $`ϵ`$ (in scheme (A)) is also shown. As we can see, the dependence on the fragmentation scheme is visible only for $`p_\mathrm{T}<3`$GeV. The rapidity distribution is mostly affected because of the relatively low $`p_\mathrm{T}`$ cut. The effect is generally modest, below 5% in the range of interest. When increasing $`ϵ`$, the $`p_\mathrm{T}`$ spectrum becomes softer. The spectrum with $`ϵ=0.036`$ is below the one with $`ϵ=0.02`$ for $`p_\mathrm{T}>1.5`$ GeV; in the tail, the difference is of the order of 30%. As in the case of the fragmentation frame dependence, the behaviour of the $`p_\mathrm{T}`$ spectrum affects the $`\eta `$ distribution through the $`p_\mathrm{T}`$ cut. When changing $`ϵ`$ from 0.02 to 0.036, the $`\eta `$ spectrum is almost unchanged in shape, and the normalisation decreases of about 10%.

In figure 2 we show the difference in the pseudorapidity and transverse momentum distributions when using the mass of the quark (our default choice) for the on-mass-shell condition, as opposed to the mass of the $`D^{}`$ meson ($`m_D^{}=2.010`$GeV). The comparison is performed in scheme (B), since in scheme (A) there is no difference in $`p_\mathrm{T}`$ and $`\eta `$ for the two mass choices. The effect is only visible below $`p_\mathrm{T}=5`$GeV. In the case of HERA experiments, only $`p_\mathrm{T}>2`$GeV is relevant; in this region, the variation is below 10%. Furthermore, larger relative variations are observed only for very small $`p_\mathrm{T}`$, where the differential cross section is actually smaller in absolute value.

## 4 Dependence upon input parameters

We now turn to the study of the sensitivity of our result to the chosen range of parameters. We begin by showing in fig. 3 the sensitivity to the choice of $`\mathrm{\Lambda }_{\mathrm{QCD}}`$.

As we can see the sensitivity is quite small, below 10%. This is a consequence of the fact that the allowed range for $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ has become fairly narrow at present.

We do not include in our study the correlation of $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ with the parton densities, since PDF sets with flavour matching conditions consistent with our scheme that allow for $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ variation are not available at present. For the same reason, we do not include a study of the uncertainty due to the choice of the PDF set. In general, the inclusion of such an uncertainty is a highly non trivial problem, and it is far from having a satisfactory solution at present. However, it must be said that current PDF sets have rather similar gluon densities in the $`x`$ region probed by charm production at HERA; therefore, we do not expect the PDF dependence to be a major source of uncertainty in our case.

The mass sensitivity is shown in fig. 4. The cross section is most sensitive to the mass value for $`p_\mathrm{T}`$ of the order of 1 GeV. In the visible range ($`p_\mathrm{T}>2`$ GeV) the sensitivity is modest, of the order of 10%.

It has to be pointed out that our study of mass sensitivity is somewhat incomplete. In fact, we do not have the possibility of varying the charm mass in the evolution of the parton densities. Studies performed in the hadroproduction case have indicated that these effects are small.

The renormalization and factorization scale dependences are displayed in figures 5 and 6 respectively. The factorization scale choice has little impact on our results, except for small transverse momenta. In particular, when $`\mu _\mathrm{F}=\mu _0/2`$, at small $`p_\mathrm{T}`$ the cross section becomes unphysical, since it depends upon the parton densities at low scales, a region where current parametrizations are unreliable. We simply truncate the theoretical predictions at small transverse momenta in this case. The renormalization scale dependence is instead sizable in the whole production range. Since no data are available for $`p_\mathrm{T}<2`$ GeV, the inspection of figures 3-6 allows us to conclude that the renormalization scale dependence is the largest source of uncertainty when comparing the theoretical results to the data.

## 5 Comparison with data

In this section, we compare the theoretical predictions, obtained with our FONLL formalism, with the data relevant to $`D^{}`$ production, presented by the H1 and Zeus collaborations in refs. . For the sake of this comparison, we now need to multiply the NPFF of eq. (4) by a factor of 23.5 %, to take into account the $`cD^{}`$ fraction . Furthermore, we multiply by a factor of 2 in order to get the $`D_{}^{}{}_{}{}^{+}+D_{}^{}{}_{}{}^{}`$ cross section, which is what the experiments quote in refs. . H1 data are presented in ref.  as $`\gamma p`$ cross sections, obtained by dividing the measured photoproduction $`ep`$ cross sections by the Weizsäcker-Williams flux; here, we quote the data as $`ep`$ cross sections.

In QCD, it is very often the case that the theoretical uncertainties are larger than the experimental errors. Thus, an estimate of these uncertainties is mandatory for a meaningful comparison with the data. Usually, one takes the maximum and minimum values of the theoretical predictions, in the allowed range of the input parameters, as an estimate of the upper and lower theoretical errors. In our case, since the computation is numerically intensive, this procedure is too time consuming. We have thus adopted the following procedure. We construct an error band for our predictions by adding linearly all the variations that we have discussed in the previous section. In other words, for each parameter $`x`$ (mass, $`\mathrm{\Lambda }_{\mathrm{QCD}}`$, $`\mu _\mathrm{R}`$, $`\mu _\mathrm{F}`$) that has central value $`x_0`$ and upper and lower values $`x_h`$ and $`x_l`$ respectively, we define an upper and lower error

$`\delta _x^{\mathrm{up}}`$ $`=`$ $`\mathrm{max}[\sigma (x_0),\sigma (x_h),\sigma (x_l)]\sigma (x_0),`$ (5)
$`\delta _x^{\mathrm{down}}`$ $`=`$ $`\sigma (x_0)\mathrm{min}[\sigma (x_0),\sigma (x_h),\sigma (x_l)],`$ (6)
with all the other parameters in the cross section kept fixed to their central values. We then define our upper (lower) end of the band as the central cross section plus (minus) the sum of all upper (lower) errors. By interpreting eqs. (5) and (6) as true errors affecting the theoretical predictions, we could have also adopted the strategy of adding them quadratically. We refrain from doing this, since it is actually very doubtful that they obey a Gaussian law. However, it must be stressed that, in practice, the quadratic sum differs only by a tiny amount from the linear one, since they are both dominated by the effect due to the renormalization scale, as discussed in the previous section.

The lower error due to the factorization scale dependence is undefined for small momenta. Thus, the lower end of the band is undefined for low momenta if the factorization scale dependence is included in the error. We therefore present one band that includes the factorization scale error, and one band that does not include it. Since the factorization scale dependence is in general rather small, the difference between the two bands is non-negligible only in the small-$`p_\mathrm{T}`$ region, basically outside the region visible to experiments.

We begin with Zeus data from ref. , for which we present a comparison with our FONLL predictions in fig. 7. As often stated in the literature, the data lie above the theoretical predictions obtained with the default choice of parameters. On the other hand, we see that the data are marginally consistent with our upper band. However, the data seem to indicate a $`p_\mathrm{T}`$ spectrum harder than that suggested by QCD. As far as the pseudorapidity distribution is concerned, QCD appears to do a decent job, except in the positive-$`\eta `$ region. The rightmost data points, however, are by far the ones affected by the largest errors.

The data from H1 are compared to our calculations in fig. 8, for the ETAG33 sample, and in fig. 9, for the ETAG44 sample. In the case of the $`y`$ spectra relevant to the cuts $`2.5<p_\mathrm{T}<3.5`$ GeV, $`3.5<p_\mathrm{T}<5`$ GeV, and $`5<p_\mathrm{T}<10.5`$ GeV, we multiply the data by the width of the $`p_\mathrm{T}`$ bin size, since they are originally quoted in ref.  as $`d\sigma /dydp_\mathrm{T}`$ cross sections. The $`p_\mathrm{T}`$ spectra are perfectly described by our FONLL predictions; all the data points are consistent with the default predictions within one standard deviation, with the marginal exception of the largest-$`p_\mathrm{T}`$ point in the ETAG44 sample. Shape-wise, the $`y`$ spectrum of the ETAG33 sample is also statistically compatible with QCD predictions; the same can be said for the ETAG44 sample, although in this case the statistical significance is clearly lower.

In order to assess more quantitatively the size of the agreement (or disagreement) between theory and data, we now present a comparison between the two in a different form. We compute, for each data point, the corresponding FONLL default prediction, and we consider the ratio data/FONLL. The results are presented as open points in fig. 10 for the $`p_\mathrm{T}`$ spectra, and in figs. 11 and 12 for the $`\eta `$ and $`y`$ spectra. The figures also display bands, whose pattern is identical to that of the bands in figs. 79. The current bands are in fact obtained by computing the ratios of FONLL bands over FONLL default predictions.

The lower inset of fig. 10 clearly shows that four (out of six) data points relevant to the $`p_\mathrm{T}`$ spectrum measured by Zeus are more than one standard deviation away from the upper end of the theoretical band. The measured spectrum appears to be harder than predicted; however, the effect is probably statistically not significant in the case in which the lowest-$`p_\mathrm{T}`$ measurement is excluded from the data sample. All but the lowest-$`p_\mathrm{T}`$ point are more than 50% away from the default prediction, but still within a factor of 2. On the other hand, there appears to be a nice agreement between H1 ETAG33 and FONLL predictions; ETAG44 data are also reasonably reproduced, but there the large errors prevent from reaching firm conclusions.

The case of the $`\eta `$ and $`y`$ spectra is dealt with in figs. 11 and 12. As far as the Zeus data are concerned, fig. 11 confirms the trend already visible in fig. 7: for all the $`p_\mathrm{T}`$ cuts, the measured distributions are more enhanced in the positive-$`\eta `$ region than FONLL predicts. This trend is however weaker if the rightmost data points, those affected by the largest errors, are excluded from the data sample. Also, by excluding these points all data are within a factor of 2 from the default theory result. In general, most of the experimental results lie inside the theoretical uncertainty band. As already in the case of the $`p_\mathrm{T}`$ spectrum, H1 data do not suggest any statistically-significant deviation from FONLL predictions.

From this study, we see that Zeus and H1 data compare differently to FONLL predictions, the former displaying some disagreements with theory that are not present in the latter. It is difficult to understand whether the two data sets are compatible within errors, since the two experiments have different visible regions, and use different observables (Zeus use $`\eta `$, H1 use $`y`$).

The error bands presented in this section do not include any of the uncertainties relevant to the non-perturbative fragmentation discussed in sect. 3. This appears to be justified for the choice of the fragmentation frame and of the on-mass-shell condition, given the size of the effects induced. The choice of the $`ϵ`$ parameter is phenomenologically more relevant, since it can give an up to 20% effect on the $`p_\mathrm{T}`$ spectrum. However, it has to be kept in mind that $`ϵ`$ is not a physical quantity; in particular, its value crucially depends upon the scheme in which the hard cross section, which is convoluted with the NPFF, is defined. According to ref. , $`ϵ=0.02`$ is the value to be used in the case of a FONLL computation in order to obtain a meaningful comparison with data. The study of sect. 3 simply suggests that a larger (smaller) $`ϵ`$ value gives a softer (harder) $`p_\mathrm{T}`$ spectrum. Thus, the comparison of FONLL predictions to Zeus data appears to imply the necessity of an even smaller $`ϵ`$ value than the LEP data fit performed in ref.  suggests. On the other hand, H1 data appear to be consistent with $`ϵ=0.02`$.

It is interesting to compare the FONLL and FO calculations. The two predictions are displayed in fig. 13 for the Zeus cuts. FONLL results are obtained with the default parameters. We show two sets of FO curves. One set (dotted curves) is obtained with default parameters except for the value of $`ϵ`$, set equal to 0.036. This is in fact the value that gives the best fit to $`e^+e^{}`$ $`D^{}`$ production data when a FO calculation is used . The other set (dashed curves) has $`ϵ=0.02`$. We include it in order to emphasise, in comparison to the FONLL result, the effect of the large-log resummation, which is only present in the latter. As expected (see sect. 1 and ref. ), FONLL and FO predictions (with the same $`ϵ`$) coincide at small $`p_\mathrm{T}`$’s; the difference between the two becomes visible for $`p_\mathrm{T}>4`$ GeV. Where the two predictions differ, FONLL must be considered superior. Notice that although none of the two predictions agree well with data, the shape is better described by the FONLL calculation. We interpret this result as evidence of the onset of resummation effects in the data. In fact, the FO result requires larger value of $`ϵ`$ in order to compensate for the lack of resummation effects.

We do not present here RS predictions, since we believe that in most of the visible range the transverse momentum is simply too small for them to apply. The RS result is below our FONLL result by 20% to 30% in the region $`4<p_\mathrm{T}<20`$ GeV, the difference becoming smaller as $`p_\mathrm{T}`$ increases. This difference is due to our choice of the $`c`$ parameter that appears in eq. (2), which we fixed to 5 on the basis of the fact that, in the hadronic component of the cross section at the NLO level, the massless limit of the cross section is a good approximation of the massive cross section only for $`p_\mathrm{T}5m`$ . We show our prediction for the extreme choice $`c=0`$ in fig. 14. The $`c=0`$ choice has a softer spectrum, and it undershoots our default choice by at most 25% in the large $`p_\mathrm{T}`$ region.

The definition of a photoproduction event in Zeus and H1 is slightly different, since the allowed photon virtuality is much smaller in H1. Thus, one is tempted to ask whether the differences of H1 and Zeus data compared to QCD calculations are due to the fact that Zeus is not completely in the photoproduction regime. The Weizsäcker-Williams approximation is valid up to non-factorizable terms of order

$$\delta =\frac{Q_{\mathrm{max}}^2}{S_{\mathrm{min}}}=\frac{Q_{\mathrm{max}}^2}{4m^2},$$
(7)
which yields $`\delta 0.1`$ for Zeus, and $`\delta 0.001`$ for H1. This does not mean that one can expect non-factorizable terms of order 10% in the Zeus case, since the leading Weizsäcker-Williams term is logarithmically enhanced by a factor of $`\mathrm{log}4m^2/m_e^217`$. Thus, one expects non-factorizable terms to be of the order of a percent in the Zeus case, and one hundred times smaller in the H1 case. It is thus unlikely that they could be the cause of the differences. A direct comparison of Zeus and H1 data, using the same kinematical cuts, would thus be useful to understand whether the two data sets are consistent among each other.

## 6 Conclusions

In this paper, we have presented comparisons between theoretical FONLL predictions and data relevant to the photoproduction of $`D^{}`$ mesons, as measured at HERA by the Zeus and H1 collaborations. We have studied the uncertainties affecting QCD predictions, and thus estimated the largest possible range for theoretical cross sections. Generally speaking, the agreement between theory and data is acceptable, after a tuning of the input parameters. While H1 data do not display any significant discrepancies with FONLL predictions, Zeus data seem to suggest a harder transverse momentum spectrum, and an excess towards the positive-pseudorapidity region. A detailed comparison between the two sets of data is not possible due to the different observables and visible regions used by the two experiments.

Having established that QCD does a reasonable job in describing charm data, it is interesting to notice that Zeus and H1 collaborations have now presented several results on bottom production , that compare much worse to QCD predictions than analogous charm results. Since QCD expectations for charm agree reasonably with data, one would expect even better agreement for the bottom quark, because of the larger mass. It is thus unlikely that the bottom excess may be attributed to a failure of QCD predictions.

###### Acknowledgments.

We wish to thank the CERN theory division for kind hospitality during the preparation of this work. We also would like to thank M. Cacciari, L. Gladilin and C. Grab for useful discussions.


text:
# A Monitor of the Focusing Strength of Plasma Lenses using MeV Synchrotron Radiation

## 1 Introduction

An experiment has been carried out at SLAC to study the focusing of electron and positron beams by thin targets of plasma . Focusing systems of this type may be used to boost luminosity at future linear colliders by sharpening their final focus. If so, a means of monitoring the operation of the lenses must be developed, at least during single beam tuning. As part of the experimental study, a monitor based on the synchrotron radiation from the plasma lens has been developed. It was used both online, as a tuning tool, and in data analysis to validate the measurements of focusing by wire scanning of the beam spots.

The SLAC Final Focus Test Beam (FFTB) , where the experiment was done, is a model of the final focus of a linear collider and is presently operated up to 28.5 GeV. With these energies, plasma lenses are strong enough that their synchrotron radiation exhibits critical energies, $`\mathrm{E}_\mathrm{c}`$, in the MeV range. The radiation is therefore quite penetrating, and the depth of penetration was used as a measure of $`\mathrm{E}_\mathrm{c}`$.

Above the absorption edges, the gamma ray absorption length increases with increasing energy until the onset of pair production which reverses this trend. Common particle calorimeter materials, like Pb, have high atomic number, Z, and so their pair production cross section rises quickly above threshold, becoming dominant above a few MeV. For these materials, the slope of the graph of penetration depth against energy quickly changes sign and higher energy gamma rays tend to interact at shallower depths.

However, the synchrotron radiation spectrum has a significant component up to $`5\times \mathrm{E}_\mathrm{c}`$. To accommodate $`\mathrm{E}_\mathrm{c}`$ in the MeV range one must use materials of low Z, such as carbon, where the dominance of pair production is postponed until energies of 50 MeV or above are reached . This is illustrated in Fig. 1.

A calibration beam of synchrotron radiation with a critical energy in the range of interest — or even a suitable monochromatic photon beam of tunable energy in this range — is rarely available, and so one must rely on particle transport codes such as EGS4  to model the detector response. At critical energies above $``$50 MeV, as could be the case with a future linear collider, these codes allow the technique to be extended into the range where the shower profile from the synchrotron radiation lengthens logarithmically with increasing energy, and low-Z material would not be advantageous.

The study of plasma lenses at presently available energies is hampered by the physical limits of practical beam-spot size measuring techniques. Wire scanning with carbon fibers has been the tool of choice, but plasma focusing can easily produce spots that instantaneously vaporize the fiber . Monitoring of the focusing strength by the use of synchrotron radiation will permit a substantial extension of the work on plasma lenses.

## 2 Experimental Setup

The experiment tested the focusing of electron and positron beams at 28.5 GeV, produced by the SLAC linac and delivered to the FFTB at $`1.5\times 10^{10}`$ particles per pulse, with pulse duration, typical of linear colliders, of 4 ps FWHM. The beam was intercepted by a 3 mm thick jet of nitrogen whose density, for the data reported here, was set to correspond to 0.25 Atm. The resulting ionization formed a plasma that focused the beam with an effective focal length in the range of millimeters , accompanied by the emission of a burst of synchrotron radiation.

The beam dimensions downstream of the jet were measured in a standard procedure  by a carbon fiber wire scanner. The procedure was to use magnetic dipoles to scan the beam, in micron steps, across a carbon fiber. The resulting Bremsstrahlung continued down the beam line along with the charged beam. Eventually the charged beam was deflected to its dump, and the gammas were allowed to leave the vacuum pipe, 30 meters from the gas jet, through a thin window. There they entered the front of a stack of absorber material in which were embedded a set of detectors. The separation between charged and neutral beams was 26 cm at this point. The final detector, after 4 radiation lengths of material, was an air Cherenkov counter , and its signal, proportional to the overlap of the beam and the fiber, was used as the wire scanner signal.

In addition to the Bremsstrahlung, there was the MeV synchrotron radiation from the plasma lens, and some synchrotron radiation, in the $``$100 keV range, from conventional beam line elements. The material before the air Cherenkov counter absorbed this radiation, and each of the detectors interspersed through the material was sensitive to the surviving flux at its own depth.

As a convenient absorber material we chose polyethylene, n(CH2), although water would have similar absorption properties. It was available in 2.5 cm thick plates. We had these cut into 25 cm squares which we stacked, along the beam direction, in blocks 22.5 cm thick. Fig. 2 illustrates the long profile of the stack. Between each block was a planar ion chamber. In total there were 8 layers, each with one absorber block and its ion chamber. Two additional chambers were inserted to allow the accelerator data collection system to monitor conditions independently. The complete stack corresponded to approximately 4 radiation lengths.

Ion chambers were chosen as the sensitive elements primarily because of the intensity of the signal. The electron density from the converted gamma rays could reach $`2\times 10^7\mathrm{cm}^2`$ per pulse, and radiation doses in the range of tens of Megarads were accumulated during the experiment. Ion chambers held out the promise of rugged performance, while requiring minimal intervention, and have excellent linearity. Preliminary tests showed pulse linearity, with a nitrogen filling, above 2 nanocoulombs $`\mathrm{cm}^2`$, beyond what would be needed for the plasma focusing experiment.

The ion chambers were of a simple double-gap window-frame construction. The bodies were constructed of three acrylic plates 3.2 mm thick, spaced by acrylic ”window-frames” 6.4 mm thick and 19 mm wide, to form two 6.4 mm collecting gaps. On the internal surfaces of the plates was bonded a laminate of 8 $`\mu \mathrm{m}`$ aluminum on 75 $`\mu \mathrm{m}`$ Mylar, the aluminum forming the electrodes of the chambers. Polyimide tape covered the edge of the aluminum, extending 6 mm in from the acrylic frame. This increased the surface path length between opposite electrodes, and so reduced the likelihood of surface discharges across the frame material. The area remaining active in each gap was then 20 cm square. Negative high voltage was applied to the outer plates, through a 1 megohm resistor for security. Signals from collected electrons were picked up from the electrodes on either side of the middle plate and taken to a preamplifier in a shielded box at the side of the chamber. The preamplifier, with a gain of 11, was based on a KH300 chip .

The gas chosen was nitrogen, and this was fed in parallel to the two gaps at a rate of a few cc per minute, the gas lines penetrating a Faraday enclosure around both the chamber and its preamplifier. The preamplifier box was also in the Faraday shield, and its signals emerged by way of a BNC connector to be carried to standard 11-bit CAMAC ADCs over about 75 m of cable, mostly RG 214.

Prototype tests indicated that the chambers could operate at 1000V with a collection efficiency that was independent of the applied high voltage, and so all were operated at the same setting, and minimal monitoring of the H.V. was needed. At this voltage, the signal at the ADC was observed to reach peak amplitude at 50 ns, and to fall within 10% of peak at about 360 ns.

## 3 Operation and Results

The carbon fiber scanning of the beam spot profile was carried out at various positions along the beam line (z coordinate) close to the focal point of the beam. The intention was to map the convergence of the beam to the waist, and its divergence beyond there. The rms spot widths were in the range 1.5 to 15 $`\mu \mathrm{m}`$. At any given z position, the beam profile was scanned in about 50 steps, of which perhaps one third would be in the peak. Although the beam was present at 10 Hz, the vacuum pumping system would permit the jet of gas to be injected only at 2 Hz. Consequently, at each scan step, data were taken for four beam pulses without any gas injection. This provided a baseline measurement of the undisturbed beam. On the fifth pulse, the gas, and its plasma focusing effect, were present, and the records from this pulse could be compared with those of the previous four. Since wire scanning is not a fast procedure, this pattern, repeated for every step of the scan, allowed us to minimize the effects of drifts of the beam parameters.

An example of the profiles obtained from an ion chamber during a scan is shown in Fig. 3. The upper graph is in the absence of nitrogen (four pulses are averaged per step). The peak is, of course, from the beam-fiber overlap, and the lack of smoothness is caused by fluctuations of the beam parameters. Details of the shape are strongly consistent at different depths in the polyethylene. The fitted peak amplitude can be used as a relative measurement of the Bremsstrahlung shower intensity at each chamber. The lower graph is from the same scan, but from pulses with the gas jet firing. The peak again is from the carbon fiber, although with only a single pulse per point, but in this plot the baseline is considerably higher. It has contributions from the beam background (also visible in the baseline of the gas-off curve), from Bremsstrahlung from the nitrogen, and from the plasma-focus synchrotron radiation. The analysis task is to separate these.

The beam background is easily removed by subtracting the baseline measured in the gas-off plot. For the first ion chamber, the “beam background” is dominated by the low energy synchrotron radiation from the conventional beam line dipoles. In fact, its ADC would saturate when strong plasma focusing was added to this. For this reason, its data are discarded from the procedure discussed below.

The Bremsstrahlung depth profile can be obtained from the ion chamber signals corresponding to the carbon fiber peak amplitudes. The shapes of the Bremsstrahlung depth profiles are consistent from run to run, within 0.4% in layer 2, and 0.2% for deeper layers. Beam intensity and spot size variations control the overall amplitude, but do not affect the profile.

To validate the procedure, an EGS4 simulation was made of a Bremsstrahlung spectrum incident on a model of the detector stack. A comparison is shown (Fig. 4) of this simulation against the results from wire scans. The average of three scans is used, and the ordinates are simply scaled to equate the average signal strength in simulated and real chambers — no fit is made. The discrepancies are taken to be estimates of the relative sensitivities of the chambers and electronics channels, with small contributions from nonuniformity in the polyethylene. They are used to obtain correction factors for each layer. An estimate of the uncertainties in the plotted points is made from variations between scans, and among different EGS4 runs.

In order to address the performance of the system for synchrotron radiation, the EGS4 code was run on the same model of the apparatus, but with a set of individual energies from 200 keV to 150 MeV. The output for each energy and layer was incorporated in a code which simulated a dipole synchrotron radiation spectrum. It interpolated the EGS4 results at each spectrum point. Results from critical energies in the range 0.5 to 10 MeV were examined.

In this experiment the final two layers saw very little signal from synchrotron radiation. They were dominated by Bremsstrahlung, and this fact allows for a straightforward subtraction of the Bremsstrahlung contribution to the other chambers. A simulated depth profile of the five relevant layers is shown in Fig. 5. Superimposed is the fit of a simple three parameter exponential absorption function $`y=y_0+a\mathrm{exp}(bx)`$, which is evidently quite adequate. The absorption coefficient, $`b`$, is obviously the inverse of the effective absorption length. The constant term, $`y_0`$, is retained empirically since it improves the fits. However, it is always found to be very small, positive, and just significant statistically. A plot of the values of the absorption coefficient obtained from fits to the simulations is shown in Fig. 6 against the critical energies of the simulated spectra. This is effectively a calibration curve.

The experimental depth profile obtained from the scan baseline (caused by Bremsstrahlung from the gas plus plasma synchrotron radiation), is illustrated in Fig. 7. The gas-off background has already been subtracted and the relative chamber sensitivities corrected for. The lines in the figure are discussed below.

The Bremsstrahlung part is now to be removed. For this purpose we note that, for critical energies in the range 1 to 10 MeV, the ratio of synchrotron radiation in layers 8 and 7 should be in the range 0.66 to 0.72 (from simulation). For Bremsstrahlung the layer 8 / layer 7 ratio should be 1.027. Equating the observed signals in each of these layers to the sum of the synchrotron radiation and Bremsstrahlung contributions, and coupling these equations with the ratios just quoted, one can simply solve for the intensities of the two components. We obtain upper and lower limits to the Bremsstrahlung intensity corresponding to the full 0.66 to 0.72 range of the synchrotron radiation ratio. With the Bremsstrahlung intensity determined for the last two layers, its full depth profile can then be subtracted from the observed signals. The separate contributions from Bremsstrahlung and plasma-lens synchrotron radiation are indicated as lines in Fig. 7.

Plasma synchrotron radiation results from a scan are illustrated in Fig. 8. The exponential absorption fit is indicated. The effect of the relative sensitivity corrections of the layers is illustrated by comparing the residuals from the exponential fits, with and without the corrections, in Fig. 9. It is evident that the corrections, made using the Bremsstrahlung depth profile, substantially improve the fit to the synchrotron radiation profile. The approximation of using a dipole synchrotron radiation spectrum to parametrize the radiation from the lens is also seen to be good.

Evaluation of $`\mathrm{E}_\mathrm{c}`$ for two scans is illustrated graphically in Fig. 10. For one of the scans, the plasma was pre-excited by a laser, and the focusing was somewhat stronger in that case. The curve is a magnified part of the calibration curve. The synchrotron radiation absorption coefficients were determined from the fits to the data as above, and are drawn as horizontal lines intersecting the calibration curve. The uncertainties obtained from the fits and from the upper and lower limits of the Bremsstrahlung subtraction are also indicated. One can read off the equivalent dipole critical energies of the two cases: $`3.86\pm 0.13`$ MeV for the less strongly focused case, and $`4.25\pm 0.20`$ MeV for the pre-excited gas. The respective absorption lengths are 42.1 cm and 43.4 cm polyethylene equivalent. In Ref. , it is concluded that $`\mathrm{E}_\mathrm{c}`$ values in this range are consistent with the focusing strengths, 0.7 T/$`\mu `$m in the x plane and 4 T/$`\mu `$m in y, estimated from the beam waist profile measurement.

## 4 Conclusions

An initial demonstration has been made of a monitoring technique for the very strong focusing of high energy particle beams made possible by plasma lenses. The technique is not difficult or expensive, is robust, and could be modified for a range of applications. It makes use of the synchrotron radiation emitted by the lens, and should work over a range $`1<\mathrm{E}_\mathrm{c}<50`$ MeV as implemented in this example. At higher energies, the logarithmic lengthening of the electromagnetic shower with gamma ray energy would be used to monitor the profile. The principal requirement is simply that the synchrotron radiation flux be separated sufficiently from the charged beam that it becomes accessible.

## 5 Acknowledgments

We thank the SLAC Plasma Lens Collaboration for their efforts in running the experiment. We acknowledge especially the work of R. Kirby, F. King, G. Collet and Y.-Y. Sung on the carbon fiber scanner, W. Craddock on the gas jet, D. Walz on support systems and R. Iverson on the beam. And we particularly thank Pisin Chen who proposed the problem of monitoring the synchrotron radiation.