Heavy Flavour Production at Tevatron and Parton Shower Effects
DESY 10134
Dec 2010
version 2(8)
DESY, Hamburg, Germany
[3mm] University of Antwerp, Antwerp, Belgium
[3mm] D.V. Skobeltsyn Institute of Nuclear Physics,
M.V. Lomonosov Moscow State University, Russia
[3mm]
Abstract
We present hadronlevel predictions from the Monte Carlo generator Cascade and numerical calculations of charm and beauty production at the Fermilab Tevatron within the framework of the factorization QCD approach. Our consideration is based on the CCFMevolved unintegrated gluon densities in a proton. The performed analysis covers the total and differential cross sections of open charm and beauty quarks, and mesons (or rather muons from their semileptonic decays) and the total and differential cross sections of dijet hadroproduction. We study the theoretical uncertainties of our calculations and investigate the effects coming from parton showers in initial and final states. Our predictions are compared with the recent experimental data taken by the D0 and CDF collaborations. Special attention is put on the specific angular correlations between the finalstate particles. We demonstrate that the final state parton shower plays a crucial role in the description of such observables. The decorrelated part of angular separations can be fully described, if the process is included.
PACS number(s): 12.38.t, 13.85.t
1 Introduction
Charm and beauty production at high energies is subject of intense studies from both theoretical and experimental points of view [1, 2, 3, 4, 5, 6]. From the theoretical point, these processes provide an opportunity to test the different predictions based on Quantum Chromodynamics (QCD) since the dominant production mechanism at high energies (i.e. small ) is believed to be quark pair production through the gluongluon fusion subprocess. At present, the problem of description of the charm and beauty production at the Tevatron within QCD is not fully solved. So, the difference between the D0 and CDF measurements [1, 2, 3, 4] of the quark and meson production cross sections and the calculations [7] performed in the framework of fixedorder nexttoleading logarithm scheme (FONLL) is about of a factor of 1.7. The central FONLL predictions [8] also lie below the data on the , , and cross sections which has been measured [5] by the CDF collaboration as a functions of their transverse momenta. Recently the calculations in the generalmass variableflavournumber scheme (GMVFNS) were performed [9, 10] for the transverse momentum distribution of and mesons and have been found to be consistent with the data [4, 5]. However, there is still no GMVFNS predictions for other measured quantities (like distributions in rapidity and azimuthal angle difference between the momenta of final mesons).
Heavy flavour production has been considered [11, 12, 13, 14, 15, 16] also in the framework of the factorization QCD approach [17]. This approach is based on the BalitskyFadinKuraevLipatov (BFKL) [18] or CiafaloniCataniFioraniMarchesini (CCFM) [19] equations for the noncollinear gluon evolution in a proton and gives the possibility to take into account large logarithmic terms (proportional to and, in the case of CCFM, also to ). A detailed description and discussion of the factorization approach can be found in [20]. A reasonable agreement between the factorization predictions and the Tevatron data on the heavy flavour production has been found in [12, 13, 14, 15, 16]. It was demonstrated [11, 15, 16] that studying the specific angular correlations between the transverse momenta of produced quarks can give an unique information about the noncollinear gluon evolution in a proton since taking into account the nonvanishing initial gluon transverse momentum in the factorization approach leads to the violation of backtoback kinematics event at leading order.
Measurements of dijet total and differential cross sections (as a functions of the leading jet transverse energy and the dijet invariant mass) and azimuthal angle correlations between two jets have been performed by the CDF collaboration [21]. These measurements significantly extend the energy range investigated by previous analyses [1, 2, 3, 4, 5]. Studying the jet production is specially interesting since there are no additional assumptions on the fragmentation of the beauty quark into the meson. The CDF collaboration has reported the preliminary data [6] on the charm pair production, where the , pair cross section and the , pair cross section as a function of the azimuthal angle between two charmed mesons have been meausred. The angular correlations in meson production have been also measured [1, 2, 3] by the D0 and CDF collaborations.
The main goal of present paper is to give a systematic analysis of all available experimental data [1, 2, 3, 4, 5, 6, 21] on the heavy flavour production at the Tevatron in the framework of the factorization formalism. We produce the relevant numerical calculations in two ways. First, we will perform analytical partonlevel calculations (labeled as LZ in the following) similar to that done in [15, 16]. In these calculations we will use the CCFMevolved gluon densities [23] as default sets. The measured cross sections of heavy quark production will be compared with the predictions of Monte Carlo event generator Cascade [24]. Cascade is a full hadron level Monte Carlo event generator for , , and processes, which uses the CCFM evolution equation for the initial state cascade in a backward evolution approach supplemented with offshell matrix elements for the hard scattering^{1}^{1}1A discussion of the phenomenological applications of Cascade can be found in [25].. In this way we will investigate the influence of parton showers in initial and final states for the description of the data. We will study the possible sources of theoretical uncertainties of our predictions (i.e. uncertainties connected with the gluon evolution scheme, heavy quark mass, hard scale of partonic subprocess and the heavy quark fragmentation functions). To investigate the dependence of our predictions on the noncollinear evolution scheme we will apply the unintegrated gluon densities derived from the usual (DGLAPevolved) parton distributions (in the KimberMartinRyskin (KMR) [26] approximation). Our special goal is to study specific kinematic properties of the final heavy quarkantiquark pair which are strongly related to the nonzero initial gluon transverse momentum.
The outline of our paper is following. In Section 2 we recall shortly the basic formulas of the factorization approach with a brief review of calculation steps. In Section 3 we present the numerical results of our calculations and a discussion. Section 4 contains our conclusions.
2 Theoretical framework
The main formulas have been obtained previously in [14, 15, 16]. Here we only recall some of them. The cross section of heavy quark hadroproduction at high energies in the factorization approach is calculated as a convolution of the offshell (i.e. dependent) partonic cross section and the unintegrated gluon distributions in a proton. It can be presented in the following form:
where is the unintegrated gluon distribution in a proton, is the offshell (i.e. depending on the initial gluon virtualities and ) matrix element squared and averaged over initial gluon polarizations and colors, and is the total centerofmass energy. The produced heavy quark and antiquark have the transverse momenta and and the centerofmass rapidities and . The initial offshell gluons have a fraction and of the parent protons longitudinal momenta, nonzero transverse momenta and (, ) and azimuthal angles and . The analytic expression for the can be found, for example, in [14, 17].
The unintegrated gluon distribution in a proton in (1) can be obtained from the analytical or numerical solution of the BFKL or CCFM evolution equations. In the numerical calculations we have tested a few different sets. First of them (A0) was obtained in [23] from the CCFM evolution equation. The initial (starting) distribution has been used in the following form:
where all input parameters have been fitted to describe the proton structure function . An equally good fit can be obtained using different values for the soft cut and a different value for the width of the intrinsic distribution (the CCFM set B0). A reasonable description of the data can be achieved [23] by both these sets.
To evaluate the unintegrated gluon densities in a proton we apply also the KMR approach [26]. The KMR approach is a formalism to construct the unintegrated parton (quark and gluon) distributions from the known conventional parton distributions , where or . In this scheme, the unintegrated gluon distribution is given by the expression
where are the usual unregulated leading order DGLAP splitting functions, and are the conventional quark and gluon densities^{2}^{2}2Numerically, we have used the standard GRV 94 (LO) [27], MSTW 2008 (LO) [28] (in LZ calculations) and MRST 99 [29] (Cascade) sets. and is the Sudakov form factor. The theta function implies the angularordering constraint specifically to the last evolution step to regulate the soft gluon singularities.
In Fig. 1 we plot the CCFM set A0 (as a solid lines), the CCFM set B0 (as a dashed lines) and KMR (as a dotted lines) unintegrated gluon densities at probing scale GeV as a function of for different values of . One can see that the shapes of gluon densities under consideration are very different from each other In the following we will study the possible manifestations in the total and differential heavy flavour cross sections.
In our analytical calculations (LZ) the multidimensional integrations in eq.(1) have been performed by the means of Monte Carlo technique, using the routine vegas [30]. The full C code is available from the authors on request^{3}^{3}3lipatov@theory.sinp.msu.ru.
3 Numerical results
We now are in a position to present our numerical results. First we describe our input and the kinematic conditions. After we fixed the unintegrated gluon distributions, the cross section (1) depends on the renormalization and factorization scales and . In the numerical calculations we set , (where is the transverse momentum of initial offshell gluon pair), GeV, GeV and use the LO formula for the coupling constant with active quark flavours at MeV, such that .
3.1 Inclusive charm and beauty production
We begin the discussion by studying the role of nonzero gluon transverse momentum in the offshell matrix elements involved in (1). In Fig. 2 we plot the differential cross section for and pair production as a function of transverse momentum , rapidity and the azimuthal angle difference between the transverse momenta of produced quarks at GeV. The solid histograms correspond to the results obtained according to the master formula in eq.(1). The dotted histograms are obtained by using the same formula but without virtualities of the incoming gluons in partonic amplitude and with the additional requirement . There are no cuts applied on the phase space of the produced quarks. As it was expected, in the backtoback region both results coincide with each other. However, we find that a sizeable effect appears at low . Therefore the nonzero gluon transverse momentum in the hard matrix element is important for the description of the data at low and mediate . This effect is more significant for charm production due to smaller . We have checked that the predictions between the LZ and Cascade calculations agree well at parton level.
Now we turn to the transverse momentum distributions of charm and beauty production. In the case of inclusive quark production, the transverse momentum distribution has been measured [1] by the D0 collaboration and has been presented in form of integrated cross section (as a function of quark minimal transverse momentum ) at the total energy GeV for (note that there is no cut on the rapidity of ). Our predictions are shown in Fig. 3 and compared to the data. We find good agreement between the LZ and Cascade predictions. We find a good description of the data when using the CCFMevolved (A0) gluon distribution. The results obtained by using the KMR gluon density are rather close to the NLO pQCD ones [31] (not shown) but lie below the data. The difference between the CCFM and KMR predictions comes from the different behaviour of these gluon densities (see Fig. 1) which is due to absence of small resummation in the KMR distributions.
The CDF collaboration has measured [4, 5] the transverse momentum distributions of and several mesons (namely, , , and ) with (where is the meson rapidities in the centerofmass frame) at GeV. Our predictions are shown in Figs. 4 – 6 in comparison with the data [4, 5]. The fragmentation of the charm and beauty quarks into a and mesons is described with the Peterson fragmentation function [32] with and . According to [33, 34], the following branching fractions are used: , , , and . The observed difference between the LZ and Cascade predictions is due to the missing parton shower effects in the LZ calculations. We address this point in more detail in Section 3.3.
Since the predicted meson transverse momentum distributions depend on the quarktohadron fragmentation, we have repeated our calculations for and mesons with the shifted values of the Peterson shape parameter , namely and . These values are also often used in the NLO pQCD calculations. Additionally, we have applied the nonperturbative fragmentation functions which have been proposed in [7, 8, 35] and which have been used in the FONLL calculations. The input parameters were determined [8, 35] by a fit to LEP data. The results of our calculations are shown in Fig. 7. We find that the predicted cross sections (in the considered region) are larger for smaller values of parameter or if the fragmentation function from [7, 8, 35] is used. However, the typical dependence of numerical predictions on the fragmentation scheme is much smaller than the dependence on the unintegrated gluon density.
The D0 experiment has measured muons originating from the semileptonic decays of quarks at GeV [1] for GeV, (for both muons) and GeV, where is the muon pseudorapidity and is the invariant mass of the produced muon pair. In [2] the measurements are extended to the forward muon rapidity region, namely . To produce muons from quarks in the LZ calculations, we first convert quarks into mesons (using the Peterson fragmentation function with default value ) and then simulate their semileptonic decay according to the standard electroweak theory. The branching of as well as the cascade decay are taken into account with the branching fraction taken from [34]. The predictions of the LZ and Cascade calculations are shown in Figs. 8 and 9. We find that our predictions with both CCFMevolved unintegrated gluon densities describe the experimental data for both the transverse momentum and rapidity distributions of muons reasonably well.
The calculated total cross sections of the quarks, and mesons and their decay muons compared to the CDF experimental data [4, 5] are listed in Table 1. In Table 2 the systematic uncertainties of our calculations are summarized. To estimate the uncertainty coming from the renormalization scale , we used the CCFM set A0 and A0 instead of the default density function A0 in Cascade. These two sets represent a variation of the scale used in in the offshell matrix element. The A0 stands for a variation of , while set A0 reflects . In all heavy flavour analyses studied here, we observe a deviation of roughly for set A0. The uncertainty coming from set A0 is generally smaller, but still positive. The dependence on the heavy flavour masses is investigated as well using Cascade. We varied our default values of GeV by GeV and GeV by GeV. All heavy flavour cross sections are observed to be lower, if the masses are enlarged and vice versa.
Source  

CDF data [b]  
A0 (LZ/Cascade)  3.47/2.76  13.34/9.31  4.53/3.45  3.81/3.19  0.49/0.35 
B0 (LZ/Cascade)  2.54/2.02  9.55/6.78  3.27/2.54  2.78/2.32  0.35/0.26 
KMR (LZ/Cascade)  1.33/1.16  6.92/5.06  2.40/1.92  2.13/1.75  0.29/0.22 
KMR (LZ, MSTW 2008)  1.16  6.04  2.17  1.79  0.25 
Source  

CCFM set A0  2.76  9.31  3.45  3.19  0.35 
CCFM set A0  +10%  +6%  +6%  +5%  +9% 
CCFM set A0  +1%  +4%  +4%  +3%  +3% 
GeV, GeV  8%  4%  4%  4%  3% 
GeV, GeV  +11%  +4%  +1%  +4%  +11% 
,  +7%  +21%  +23%  +21%  +25% 
Total 
We turn now to the investigation of the angular correlations between the produced particles. As it was mentioned above, such correlations have been not studied yet in the framework of GMVFNS scheme (which has been applied for transverse momentum distributions of the charm and beauty mesons only [9, 10]). Experimental data on the azimuthal correlations in charm and beauty production come from both the CDF and D0 collaborations. In the case of quark production, CDF data [3] refer to the azimuthal angle distribution measured at , GeV, GeV and GeV. The D0 data [1] refer to the muonmuon correlation measured in the region of GeV, and GeV at the same energy. In the case of charm production, the CDF collaboration has presented the measurement [6] of , and , pair cross section as a function of angle between the two charm mesons for the kinematic range , GeV, GeV and GeV. Our predictions are shown in Figs. 10 and 11 in comparison with the data [1, 3, 6]. We observe that the predicted shapes of azimuthal angle distributions are very different for different unintegrated gluon density functions. This is in a contrast to the cross sections as a function of transverse momenta or rapidities where all unintegrated gluon densities gave a similar behaviour. We can conclude that the cross section as a function of is very sensitive to the details of the noncollinear gluon evolution in a proton. As it was pointed out in [11, 15, 16], such observables can serve as an important and crucial test discriminating the different approaches of the small physics. The CCFMevolved gluon densities overshoot the data at and tends to underestimate them at . We observe that the peak at is not at all reproduced by the CCFM unintegrated gluon densities; we address this point in more detail in Section 3.3 where we discuss the the process . However the KMR gluon density has obviously a very different distribution (see Fig. 1) and therefore provides a better description.
3.2 dijet production
Recently, the experimental study of the dijet production in collisions at GeV has been presented [21, 22]. The total cross section has been measured in the kinematical region defined by , , GeV and GeV. The differential cross sections as a functions of the leading jet transverse energy and of the dijet invariant mass have been measured with the azimuthal angle correlation between the two jets. The jets are reconstructed with the JETCLU cone algorithm with radius . Our predictions are shown in Figs. 12 –14. We observe that at large and at large dijet invariant masses the predictions fall below the measurement. Note, however, that in this kinematical region the quark induced subprocesses (such as ) becomes important but are not taken into account here. At small and moderate and , where gluon induced subprocesses play the leading role, the overall agreement of our predictions and the data is quite good. The measured distribution has a significantly different tail at small compared to the predictions; this will be addressed in more detail in the next section.
3.3 Parton shower effects and additional gluon processes
As shown above, the full set of the data on the heavy flavour production is in general reasonably well described. However the measured cross sections as a function of show significant differences. These distributions are directly sensitive to additional parton radiation treated by parton showers or by other additional processes. In the following we investigate the influence of the details of the parton shower and the process with subsequent branching of one of the outgoing (offshell) on the distribution. For these studies we use the CCFM set A0.
In Cascade, the initial state parton shower is angular ordered and based on the CCFM evolution equations. The final state parton shower is also angular ordered, but based on the DGLAP evolution equations. To investigate the dependence on the parton shower, we calculated the cross sections without parton shower, only initial state, only final state and both initial and final state parton shower, as shown in Fig. 15. Here, dijet production is selected as an example. We observe a very small contribution of the initial state parton shower, since in factorisation the initial state parton shower does not influence the of the gluons (since it is determined from the unintegrated gluon density). However, a significant contribution comes from the final state parton shower. The prediction with full parton shower underestimates the backtoback region of the azimuthal separation of the two jets. This suggests that the final state parton showering generates too much gluon radiation, which causes less correlated jets as shown in Fig. 15 (a). To investigate this further we changed the final state parton shower scale from to . As shown in Figure 15 (b), this leads to a higher correlation of the jets and a very good description of the backtoback region of the jets is achieved.
At lower values of the azimuthal separation of the dijet system, higher order processes are expected to contribute significantly. Therefore, we repeated the dijet analysis for the process . In the matrix element calculation, which was performed in [35], one gluon in the initial state is onshell and the other one is offshell, as shown in Fig. 16. The heavy quark pair is then produced via parton showers in the final state as . As shown in Fig. 17 (a), this process contributes significantly to the tail of the distribution very well, without any additional adjustment of parameters. Part of the contribution from with is already included in the matrix element , so simply adding both contributions would result in double counting [37]. In [22] the influence of multiparton interactions in the collinear frame is investigated, and it was found, that the contribution even in the small region is small, because a large of the jets is required.
In conclusion, the azimuthal cross section measurements can be very well described by adjusting the scale for the final state parton shower to and by including the process with subsequent branching , which in a collinear calculation would contribute only at NNLO.
4 Conclusions
We have studied charm and beauty production at Tevatron energies within the framework of the factorization. Our calculations are based on the CCFMevolved unintegrated gluon densities in a proton. The analysis covers the total and differential cross sections of open charm and beauty quarks, and mesons (or rather muons from their semileptonic decays). The cross sections of dijet production have been studied also. Special attention was put on the specific angular correlations between the finalstate particles. Using full hadronlevel Monte Carlo generator Cascade, we investigated the effects coming from the parton showers in initial and final states. Different sources of theoretical uncertainties have been specially studied.
We obtain good agreement of our calculations for all observables and the recent experimental data taken by the D0 and CDF collaborations. We have demonstrated, that the parton shower plays a significant role in the description of the cross section. We obtain a good description of the correlations once the higher order process with subsequent splitting is included.
5 Acknowledgements
We thank S.P. Baranov for his encouraging interest and useful discussions. We are very grateful o S. Vallecorsa for many discussions on the CDF results. We are also grateful to S. Baranov for pointing out a possible double counting when processes for heavy quark production. The authors are very grateful to DESY Directorate for the support in the framework of Moscow — DESY project on MonteCarlo implementation for HERA — LHC. A.V.L. was supported in part by the Helmholtz — Russia Joint Research Group. Also this research was supported by the FASI of Russian Federation (grant NS4142.2010.2) and FASI state contract 02.740.11.0244.
References
 [1] B. Abbott et al. (D0 Collaboration), Phys. Lett. B 487, 264 (2000).
 [2] B. Abbott et al. (D0 Collaboration), Phys. Rev. Lett. 84, 5478 (2000).
 [3] D. Acosta et al. (CDF Collaboration), Phys. Rev. D 71, 092001 (2005).
 [4] A. Abulencia et al. (CDF Collaboration), Phys. Rev. D 75, 012010 (2007).
 [5] D. Acosta et al. (CDF Collaboration), Phys. Rev. Lett. 91, 241804 (2003).
 [6] J. Rademacker, Proceedings of Charm’07, Ithaca, NY, August 2007.

[7]
M. Cacciari and P. Nason, Phys. Rev. Lett. 89, 122003 (2002);
M. Cacciari, S. Frixione, M.L. Mangano, P. Nason and G. Ridolfi, JHEP 0407, 033 (2004).  [8] M. Cacciari and P. Nason, JHEP 0309, 006 (2003).
 [9] B.A. Kniehl, G. Kramer, I. Schienbein, and H. Spiesberger, Phys. Rev. D 77, 014011 (2008).
 [10] B.A. Kniehl, G. Kramer, I. Schienbein, and H. Spiesberger, DESY 09008.
 [11] S.P. Baranov and M. Smizanska, Phys. Rev. D 62, 014012 (2000).
 [12] Ph. Hägler, R. Kirschner, A. Schäfer, L. Szymanowski and O.V. Teryaev, Phys. Rev. D 62, 071502 (2000).

[13]
M.G. Ryskin, Yu.M. Shabelski and A.G. Shuvaev, Phys. Atom. Nucl. 64, 1995 (2001);
Yu.M. Shabelski and A.G. Shuvaev, Phys. Atom. Nucl. 69, 314 (2006).  [14] N.P. Zotov, A.V. Lipatov and V.A. Saleev, Phys. Atom. Nucl. 66, 755 (2003).
 [15] S.P. Baranov, N.P. Zotov and A.V. Lipatov, Phys. Atom. Nucl. 67, 837 (2004).
 [16] A.V. Lipatov, L. Lönnblad and N.P. Zotov, JHEP 0401, 010 (2004).

[17]
L.V. Gribov, E.M. Levin, and M.G. Ryskin, Phys. Rep. 100, 1 (1983);
E.M. Levin, M.G. Ryskin, Yu.M. Shabelsky and A.G. Shuvaev, Sov. J. Nucl. Phys. 53, 657 (1991);
S. Catani, M. Ciafoloni and F. Hautmann, Nucl. Phys. B 366, 135 (1991);
J.C. Collins and R.K. Ellis, Nucl. Phys. B 360, 3 (1991). 
[18]
E.A. Kuraev, L.N. Lipatov and V.S. Fadin, Sov. Phys. JETP 44, 443 (1976);
E.A. Kuraev, L.N. Lipatov and V.S. Fadin, Sov. Phys. JETP 45, 199 (1977);
I.I. Balitsky and L.N. Lipatov, Sov. J. Nucl. Phys. 28, 822 (1978). 
[19]
M. Ciafaloni, Nucl. Phys. B 296, 49 (1988);
S. Catani, F. Fiorani and G. Marchesini, Phys. Lett. B 234, 339 (1990);
S. Catani, F. Fiorani and G. Marchesini, Nucl. Phys. B 336, 18 (1990);
G. Marchesini, Nucl. Phys. B 445, 49 (1995). 
[20]
B. Andersson et al. (Small Collaboration), Eur. Phys. J. C 25, 77 (2002);
J. Andersen et al. (Small Collaboration), Eur. Phys. J. C 35, 67 (2004);
J. Andersen et al. (Small Collaboration), Eur. Phys. J. C 48, 53 (2006).  [21] T. Aaltonen et al. (CDF Collaboration), CDF note 8939.
 [22] S. Vallecorsa. Measurement of the bbar dijet cross section at CDF. PhD thesis, University Geneve, 2007.
 [23] H. Jung, arXiv:hepph/0411287.

[24]
H. Jung, Comp. Phys. Comm. 143, 100 (2002);
H. Jung et al., DESY 10107.  [25] F. Hautmann and H. Jung, JHEP 10, 113 (2008).

[26]
M.A. Kimber, A.D. Martin and M.G. Ryskin, Phys. Rev. D 63, 114027 (2001);
G. Watt, A.D. Martin and M.G. Ryskin, Eur. Phys. J. C 31, 73 (2003). 
[27]
M. Glück, E. Reya, and A. Vogt, Phys. Rev. D 46, 1973 (1992);
M. Glück, E. Reya, and A. Vogt, Z. Phys. C 67, 433 (1995).  [28] A.D. Martin, W.J. Stirling, R.S. Thorne, and G. Watt, Eur. Phys. J. C 63, 189 (2009).
 [29] A.D. Martin, R.G. Roberts, W.J. Stirling, and R.S. Thorne, Eur. Phys. J. C 14, 133 (2000).
 [30] G.P. Lepage, J. Comput. Phys. 27, 192 (1978).
 [31] M. Mangano, P. Nason and G. Ridolfi, Nucl. Phys. B 373, 295 (1992).
 [32] C. Peterson, D. Schlatter, I. Schmitt and P. Zerwas, Phys. Rev. D 27, 105 (1983).
 [33] R. Barate et al. (ALEPH Collaboration), Eur. Phys. J. C 16, 597 (2000).
 [34] K. Hagiwara et al. (PDG Collaboration), Phys. Rev. D 66, 010001 (2002).
 [35] E. Braaten, K.M. Cheng, S. Fleming and T.C. Yuan, Phys. Rev. D 51, 4819 (1995).
 [36] M. Deak, F. Hautmann, H. Jung, K. Kutak, JHEP 0909, 121 (2009).
 [37] S.P. Baranov, private communication