Exclusive heavy quark-pair production in ultraperipheral collisions

In this article, we study the fully differential observables of exclusive production of heavy (charm and bottom) quark pairs in high-energy ultraperipheral $pA$ and $AA$ collisions. In these processes, the nucleus $A$ serves as an efficient source of the photon flux, while the QCD interaction of the produced heavy-quark pair with the target ($p$ or $A$) proceeds via an exchange of gluons in a color singlet state, described by the gluon Wigner distribution. The corresponding predictions for differential cross sections were obtained by using the dipole $S$-matrix in the McLerran--Venugopalan saturation model with impact parameter dependence for the nucleus target, and its recent generalization, for the proton target. Prospects of experimental constraints on the gluon Wigner distribution in this class of reactions are discussed.


I. INTRODUCTION
In QCD, the hadron structure is encoded in the so-called Wigner distributions [1][2][3]. These distributions are known to provide the most detailed information on the parton multidimensional imaging (tomography) in the target. The 5D Wigner distribution W (x, k ⊥ , b ⊥ ) depends on both the transverse momentum k ⊥ of an exchanged parton and its impact parameter b ⊥ . While the Wigner distribution is in impact-parameter b ⊥ representation, their Fourier transform as b ⊥ → ∆ ⊥ is known as the generalized transverse momentum distribution (GTMD) [4][5][6][7] in momentum representation.
These distributions are therefore sensitive to the angular correlation between b ⊥ and k ⊥ whose magnitude is determined by the elliptic Wigner distribution [8][9][10]. It was shown earlier that the angular dependence of the Wigner distribution is particularly responsible for an elliptic flow in pA collisions [11,12], the angular correlations in deeply virtual compton scattering [13] and in quasi-elastic photon-nucleus scattering [10], etc. For a comprehensive review on the fundamental role of these distributions, see also Refs. [14,15] and references therein.
In Refs. [8,16] considering an important example of electron-ion collisions in the highenergy limit it was demonstrated that the low-x GTMD is directly related to a Fourier transform of an impact parameter dependent forward dipole amplitude (or dipole S-matrix), S( r, b) providing an important connection with the gluon saturation phenomena at low-x (for a detailed review of the saturation effects and the Color Glass Condensate, see e.g. Ref. [17]). Just like for lower-dimensional descendents, the collinear parton densities, the QCD perturbation theory cannot predict the key characteristics of the partial dipole amplitude S( k ⊥ , ∆ ⊥ ) and hence the gluon Wigner distributions, so potential possibilities for experimental measurements of such distributions or for setting constraints on them directly from the data gain large importance and have to be studied in detail [8]. The basic difficulty on the extraction of the Wigner distributions (or GTMDs) is typically associated with the fact that the differential cross section is not proportional to a GTMD itself but is given by its convolution integral with the light-cone wave function for a given projectile Fock state scattering off a target. Such integral is originated as a remnant of the loop integral in the exclusive production amplitude formed by the two exchanged gluons with the target (in the color-singlet state), and it is in general not analytically invertible.
A particular relevant class of scattering processes at hadron colliders, the high-energy ultraperipheral collisions (UPCs) with fully exclusive final states, provide essential means for accessing the hadron structure at relatively low momentum transfers and at low-x due to both clean environment and complete reconstruction of kinematics of the exchanged gluons with a target. In UPCs, the relativistic systems scatter at typically large impact parameters by means of quasi-real Weiszäcker-Williams (WW) photon exchange [18,19]. The WW flux is enhanced for heavy nucleus as the square of its charge making it an efficient source of photons. As was demonstrated recently in Ref. [20], such gluon Wigner distribution can be constrained, or even directly extracted, from the data on exclusive light-quark di-jet photoproduction in the UPCs. Besides the largest contribution to the di-jet photoproduction signal, the use of light quarks in the final state as a source of di-jets is particularly beneficial as this channel enables to directly extract the gluon Wigner distribution the data on fully differential exclusive di-jet cross section. Indeed, the loop integral is analytically invertible in this particular case such that the components of the gluon Wigner distribution can be found as integrals of the components of the differential cross section.
However, a relevant non-trivial structure of the gluon Wigner distribution, in particular, its elliptic component, emerges when the transverse momenta of the produced q andq are relatively low and do not significantly exceed the saturation scale of the process. This further prompts valid questions about the applicability of the QCD perturbation theory for reliable computation of the γ + (gg) → qq matrix element in such a region of predominantly soft or semi-soft kinematics. As the only hard scale of this process is associated with the transverse momentum of a produced jet (in dominant nearly back-to-back di-jet configurations), going to low (below a few GeV) jet transverse momenta is severely restricted by potentially large higher order effects, thus, limiting the capability of this method for a reliable extraction of the gluon Wigner distribution in the domain of its maximal enhancement and focussing on the tail high-P ⊥ regions only.
Another technical challenge is related to experimental capabilities for exclusive di-jet photoproduction measurements at high energies. Such a measurement requires reconstruction of full di-jet kinematics simultaneously triggering on large rapidity gap events only in order to suppress backgrounds causing the leakage of energy and transverse momentum into unreconstructed hadronic activity originating from the break up of the target, of the exchanged Pomeron and/or of the quasi-real photon. For this purpose, the forward proton (or nucleus) reconstruction is needed to ensure exclusivity of the corresponding diffractive reaction. While ATLAS Forward Physics program offers certain possibilities for such a measurement, the acceptance in jet transverse momenta is highly limited to high-P ⊥ kinematics only, with the lowest cut-off hardly going below 20 GeV or so.
FIG. 1: Feynman diagram for a quark-pair photoproduction in UPCs. The projectile nucleus denoted by a shaded blob is a source of quasi-real WW photons. The photon fluctuates into a qq dipole, which interacts by means of two-gluon exchange in a color singlet state with the target proton or nucleus.
Such situation indicates that the use of heavy (c and b) quarks, with measurements of exclusive open heavy flavor (DD and BB meson pairs) photoproduction, may help in resolving both theoretical and experimental issues with probing the gluon Wigner distribution in the gluon saturation regime. First of all, it naturally provides a hard scale associated with the heavy quark mass, thus, enabling the use of QCD perturbation theory, with less degree of uncertainty, even for a vanishing quark transverse momentum. Second, reconstruction of heavy-flavored mesons can be performed at much lower transverse moments than for jets and, despite of a smaller cross section, with a much better control of QCD backgrounds. As a price to pay for such an improvement, for heavy quarks the convolution integral in the diffractive amplitude is no longer analytically invertible, so for now we can only make predictions for the corresponding observables in the framework of a given model for the gluon Wigner distribution.
These arguments motivate us to perform the first detailed study of exclusive cc and bb pairs photoproduction in UPCs in fully resolved kinematics when the target survives the interaction and is detected in a forward region as shown in Fig. 1. Such process is driven by at least the exchange of two gluons in a color singlet state, to the leading order in QCD perturbation theory. This study is performed for a large nucleus target in the framework of the recently upgraded McLerran-Venugopalan (MV) incorporating the gluon saturation phenomenon with impact parameter dependence [11,21]. For the case of proton target, we have employed a recent generalization of the MV model presented in Ref. [11]. We analyze the corresponding observables particularly sensitive to the elliptic gluon Wigner distribution in the physically relevant regions of exchanged gluon transverse momenta close to the saturation scale. While a direct reconstruction of the Wigner distribution from the data remains a challenging task, we identify certain non-trivial behavior in the differential cross section directly related to the features of the elliptic distribution.
The paper is organized as follows. In Section II, the light cone approach for heavy quark exclusive photoproduction is formulated in terms of the dipole S-matrix, and that is a new analytic results as far as we know. In Section III, we provide a short description of the MV model for both the large nucleus target and the proton target. Section IV contains numerical results for the structure functions for the massive and massless cases as well as the production cross sections for cc and bb in the case of lead and protons targets. Finally, concluding remarks and a summary are given in Section V.

A. Kinematics
We work in the nucleus-target center of mass frame. The photon is collinear to the z-axis direction, and carries energy ω. Since we are working with quasi-real photons (q 2 ≈ 0), its momentum can be written as q = ( √ 2ω, 0, 0 ⊥ ) in Sudakov (light-cone) variables and the longitudinal polarization contribution to the cross section is negligible.
The incoming gluon transverse momenta are denoted as k ⊥ − ∆ ⊥ /2 and − k ⊥ − ∆ ⊥ /2, where k ⊥ is integrated as the loop momentum. Their energy and longitudinal momenta are neglected in the limit of x → 0. In what follows, we do not consider QCD evolution in x variable (or rapidities). Instead, we perform our analysis in the forward kinematics, such that the gluonic contribution to the quark longitudinal momentum is small, and we justify neglecting it. In this approximation the target contributes only with a total transverse momentum − ∆ ⊥ to the final qq pair.
The final-state heavy quarks studied here are charm and bottom quarks, with masses m c = 1.4 GeV and m b = 4.7 GeV, respectively. The quark will carry momentum fraction z of the projectile photon and has transverse momentum − P ⊥ + k ⊥ , coming from the photon splitting, while the antiquark will carry momentum fraction (1 − z) from the photon and have transverse momentum P ⊥ − k ⊥ . After the di-gluon exchange, the quark will acquire the following light-cone momentum components and, analogously, for the antiquark The k − i (where i = 1, 2) components are determined from the condition that the final states are on mass shell. If the quark rapidities y i =ln( √ 2k + i / k 2 i⊥ + m 2 Q ) and transverse momenta are measured, then the quark momentum fraction z and the photon energy ω are fixed.

B. Exclusive qq photoproduction cross section
The cross section for a UPC between a projectile nucleus A and a target, which can be either a nucleus or a proton, i.e. j = A, p, can be written as dσ Aj where ωdN γ /dω is the photon number density. These quasi-real photons coming from the projectile heavy ion are modeled by the WW photon distribution [18,19,22] In the above expression, Z is the atomic number of the projectile, α is the fine structure constant, and ξ jA = ω(R j + R A )/γ is defined in terms of the Lorentz factor γ = √ s jA /2M p , the target and nucleus radii, R j and R A , respectively, and the jA center-of-mass energy, √ s jA . The R j + R A dependence guarantees that the photons can only interact with the target when there is no overlap between the projectile and the target in impact parameter space. For a review on peripheral collisions and photon fluxes, see Ref. [23]. The partonic cross section was calculated using the light-cone Feynman rules [24,25]. The photon-splitting into a qq dipole can be calculated analytically via the photon wave-function, while the di-gluon-dipole interaction is encoded in the transition matrix T . The latter is defined as T = 1 − S in terms of the dipole S-matrix describing the quark-antiquark dipole scattering off the target and will be discussed in the following sections. For the parton-level cross section, we have the following expression where the phase space element is given by dPS = dy 1 dy 2 d 2 P ⊥ d 2 ∆ ⊥ . The functions M 0 and M 1 are expressed in terms of the T -matrix as follows and Formally, in the above expression, it is S that should be in place of T . However, we took out the non-interaction term 1 from the dipole S-matrix, so for vanishing P ⊥ we will have M 0 = M 1 = 0 (see e.g. Ref. [20]).

C. T matrix
It has been shown that azimuthal angular correlations between k ⊥ and ∆ ⊥ are important in the T -matrix when working in the saturation regime [16,26,27]. In the limit where P ⊥ ∆ ⊥ , these correlations can be taken into account by performing an expansion in Fourier harmonics The main contribution to the cross section comes from the isotropic part T 0 , and we have a sub-leading contribution from the elliptic part T [8]. The latter is the subject of our further analysis, while the higher order harmonics are suppressed.
Since the T matrix contains information on the strong interaction in the saturation regime, it is unfeasible to use pQCD to calculate it, and we must utilize a physically reasonable model for it that should incorporate the gluon saturation effects. Typically, such models are formulated in impact parameter space, where we define r ⊥ to be the dipole size and b ⊥ -the impact parameter in the transverse plane. To connect the model with our representation in momentum space, we define the Fourier transform as This integral can be related to the isotropic and elliptic contributions by expanding the Fourier exponentials in Bessel functions of first kind according to the following relation The isotropic part is related to n = 0, while the elliptic term appears for n = ±2. The odd terms do not contribute to UPCs, but it should be noted that they can be of importance in other processes where the nucleus exhibits an inhomogeneity in the radial direction [6,28]. As a result, we obtain (2.10) Here, the Gaussian-type exponentials act as dumping terms to improve the convergence of the integrals, which are highly oscillatory on the periphery. Physically, they provide a physical cut-off accounting for confinement effects, therefore, the parameters are inversely related to the typical size of the bound systems (nucleon and/or nucleus) and thus should be small compared to the hard scale of the process. We expect the b ⊥ parameter to be smaller than the target size, such that b = 1/R 2 j . We also set r = (0.5 fm) −2 as the photon splitting into the quark-antiquark pair will be suppressed when the transverse separation within the pair is larger than a typical hadron length-scale.

D. The structure functions in the massive quark case
Since the angular dependence is explicit, we can calculate analytically the azimuthal integrals in Eqs. (2.4) and (2.5). The first expression is then reduced to: which is now written in terms of the following two structure functions These generalize the results of Ref. [20] to the massive case. On top of that one should consider two additional structure functions C and D entering the amplitude in Eq. (2.5), in proportion to the quark mass, such that where (2.16) The C and D are defined so that all structure functions have the same mass dimension Starting from the parton-level cross section in Eq. (2.3), the above results enable us to represent the hadron-level cross section as follows (2.17)

III. PHENOMENOLOGICAL MCLERRAN-VENUGOPALAN MODEL FOR THE DIPOLE S MATRIX
In the MV model [21,29], a heavy nucleus can be treated as a semi-classical color field, where the gluons have a high occupation number which is controlled by the saturation scale Q s . In a recent proposal [11], Iancu and Rezaeian generalize the MV formalism for the target in order to incorporate a non-trivial impact-parameter dependence following the other saturation models such as IP-Sat and bCGG [30,31]. Such a generalization enables to analyze azimuthal asymmetries in heavy-ion collisions as a consequence of the collective phenomena in the initial state. In what follows, we employ this model in studies of heavyquark photoproduction in UPCs. We will work in the picture where the dipole experiences multiple soft scatterings off the target. The T -matrix is, in the Glauber approximation, where N ( b ⊥ , r ⊥ ) is the single dipole scattering amplitude. We use A = 1 for the proton target. The angular correlations in configuration space are generated by expanding the single scattering amplitude in Fourier harmonics, as done in the previous section For a nucleus target, one gets [11] and In the above equations, the parameter m g is an IR regulator, which acts as an effective gluon mass. We set it to 0.25 GeV. We define the saturation scale at zero impact parameter as Q 0,s = 1/R, where R is a scale related to the width of the proton color-charge distribution. We use R = 2 GeV −1 , based on best fit values obtained in different saturation models to the HERA and NMC data [32,33]. The proton radius is fixed as R p = 0.8 fm, whilst that of the nucleus is given by R A = (1.12 fm)A 1/3 . The T A is the thickness function of the nucleus target found as follows where ρ A (r) is the Woods-Saxon nuclear density distribution, ρ A ( r) = N A [1 + exp ((r − R A )/δ)] −1 , and N A is determined by the normalization condition d 3 rρ A ( r) = 1. In numerical analysis, we use δ = 0.54 fm. The T -matrix for a lead (Pb, A = 208) nucleus is plotted in Fig. 2. The isotropic part falls very fast with ∆ ⊥ and has a larger contribution from soft gluon momentum (k ⊥ 2 GeV). The elliptic part falls slower, becoming more important at higher ∆ ⊥ . Note that its peak is kept at an almost constant k ⊥ , as expected, since this is controlled by the saturation scale.
For a proton target, the single scattering amplitudes are given by an integral over the relative transverse momentum between the soft gluons p ⊥ , and we have evaluated analytically one of the integrals given in Ref. [11]. Namely, The Gaussian-like exponentials in impact parameter, b ⊥ , and relative momentum, p ⊥ , are introduced in the model due to the color distribution in the transverse plane of the proton. The first term of Eq. (3.5) is the usual one presented in the original MV model [21]. We present the T -matrix for the proton in Fig. 3. While the isotropic component falls much slower then in the case of a nucleus target, the elliptic one rises up to ∆ ⊥ ≈ 0.7 GeV, and then starts to decrease.

A. Structure functions
As one of the main new results of our analysis, we show the structure functions A, C, B, D in Figs. 4, 5, 6, 7, respectively, for heavy (c, b) quarks exclusively produced off a lead (Pb) nucleus target as functions of the quark-antiquark relative transverse momentum, P ⊥ . Since the A, B structure functions are present also for massless quark case, previously studied in Ref. [20], these are shown by dashed-dotted curves in addition to the corresponding ones for heavy quarks. These results are shown for two distinct values of ∆ ⊥ = 0.1 and 0.2 GeV, on left and right panels, respectively. Due to the presence of the additional structure functions, the information that can be obtained by probing the nucleus structure by means of heavy quarks is clearly richer than for the case of massless quarks. When analyzing absolute values, one notices that the structure function C is comparable to A, both determined in terms of the isotropic Wigner distribution, as is D to B, given in terms the elliptic Wigner distribution. Therefore, the introduction of quark masses is not a matter of only correcting the standard A and B functions but surely the new C and D structure function considered here for the first time must be included.
Figs. 4-7 show that, as quark mass increases, absolute values of the peaks of the structure functions get smaller. This behavior is similar to a rescaling of P ⊥ accompanying a "stretch" of the shape of the distributions along the horizontal axis such that the decreasing (due to a mass suppression) peaks of the structure functions move to larger P ⊥ . More specifically, in the large quark mass limit, the position of the peaks approximately scales with the quark mass as P peak ⊥ ∼ m Q , while this dependence is violated for small quark masses. The effect of going from ∆ ⊥ = 0.1 to 0.2 GeV is basically a change of the sign and a reduction of the absolute value. However, for larger ∆ ⊥ , the elliptic structure functions B and D become more relevant when compared to A and C, as the latter have a stronger reduction in absolute values. These features are a direct consequence of oscillations and the general behavior of the T 0 and T functions in the MV model, as it is shown in Fig. 2. It is important to notice that the P ⊥ values of the peaks are kept almost constant for the different ∆ ⊥ . This is the same as for the T -matrix: the peaks are controlled by the saturation scale.
FIG. 6: Same as Fig. 4, but for the B structure function. Analogously, we show the structure functions for the proton target in Figs. 8 and 9 for heavy (c, b) and massless quarks as functions of the relative quark-antiquark momentum P ⊥ for ∆ ⊥ = 0.2 GeV. These are also new results, but now regarding the proton structure. The proton has structure functions with smaller absolute values than the nucleus as our results are not normalized by the number of nucleons.
The absolute value of structure functions C and A are again comparable, while in here D is even larger than B. The behaviour w.r.t. the quark masses is analogous to the nucleus case.
In contrast to the nucleus, the MV model for the proton target presents a smaller difference in the A and C structure functions for different ∆ ⊥ , i.e. just a small decrease in the absolute value and no sign flip within the kinematic range considered here. That is why in this case we chose not to show any figures with different ∆ ⊥ , as such a behavior can be directly inferred from Fig. 3.
Also for larger ∆ ⊥ , the elliptic structure functions B and D become even more relevant than in the nucleus case. Comparing Figs. 2 and 3, we see that they will increase with ∆ ⊥ in the proton case, while in the nucleus case they decrease. Of course, this depends on the  ∆ ⊥ range we work with. If very large, the proton structure functions may show a similar behavior. Our choices of ∆ ⊥ are reasonable within the detector constraints available at the LHC.

B. Exclusive quark-pair photoproduction cross sections
Now we present our cross section results. We calculate the hadron cross section integrated in angle with exact kinematics. However, it is instructive to understand what happens in the limit k 1,2⊥ → P ⊥ . In this limit, azimuthal integration produces terms proportional to 2A 2 + B 2 or 2C 2 + D 2 . We numerically investigated this approximation and found that it has negligible impact on the final result for our choice of small ∆ ⊥ . For all plots we chose the quarks rapidities y 1,2 equal to 1. In this way, we work in the forward region accessible by ATLAS and CMS (which also justifies our choices of c.o.m. energies). This choice is justified also by the fact that the MV model is fitted for low x. Therefore, the forward moving nucleus will provide the photons while the target (backward moving) will provide the small x gluons. First, we present our cross section results for bb production in lead-lead UPCs, where A = 208 and Z = 82. The hadron cross section is depicted in Fig. 10 as a function of P ⊥ . Two plots are shown for ∆ ⊥ = 0.1 and 0.2 GeV, on left and right panels, respectively. The behaviour at large P ⊥ resembles an exponential decay and is expected because MV model does not take into account QCD evolution.
When looking closer at the ∆ ⊥ variation, it is seen that there is no substantial difference in the shape of the cross section other than it is roughly divided by a factor of the order of 100. In the previous section, it was shown that in going from ∆ ⊥ = 0.1 to 0.2 GeV the sign of all the structure functions flips and their magnitudes reduces by a factor of the order of 10. As the cross section is dominated by A 2 , the sign is irrelevant. The observation that there is a larger cross section for smaller ∆ ⊥ is just a result of the fact that the quark pair will more likely be created in a back-to-back configuration.
In Fig. 10 we also show what happens if one sets to zero some of the structure functions. For instance, with C = D = 0, the elliptic part B plays an insignificant role compared to the A contribution. The same can be said about the other elliptic structure function D. However, when C is turned on, there is a small but significant difference at P ⊥ ≈ 10 GeV and a relatively large contribution at small P ⊥ . Therefore, the contribution proportional to C 2 can be directly measured with an appropriate choice of kinematical cuts. In Fig. 11 charm production in Pb+Pb collisions for ∆ ⊥ = 0.2 GeV is depicted. Overall, it has the same features as for bottom production, but it is a larger cross section, especially at small P ⊥ when it becomes a factor 10 3 larger. Another point is that at P ⊥ ≈ 10 GeV the relative size of the C structure function contribution is larger when compared to the bottom quark case. So it appears that charm quark pair production would be a better observable than bottom quark-pair for our purposes. It has a much higher cross section and it discerns the C structure function better. However, there is an important issue: quark-meson fragmentation. The detectors will be able to measure transverse momentum of D and B mesons only, i.e., they do not access quark level variables. In spite of that, for a large heavy quark mass, there is a well known effect: the produced meson will have most of the momentum of the heavy quark. Therefore, the bottom quark observables have an advantage that fragmentation does not washes away too much the momentum distributions [34]. For the charm quark case, the c → D fragmentation may be relevant for a comparison with future measurements. But since the corresponding calculations become very difficult to perform numerically, and their theoretical interpretation becomes less transparent, we leave this point for future studies. The next Fig. 12 shows the same bb and cc production cross sections but as functions of the target final transverse momentum ∆ ⊥ . The variable P ⊥ is fixed at 10 GeV. We remark that the dips (minima) are not affected by changing the c.m. energy. They are a direct result of the oscillations in MV model, whose scale is related to the saturation scale. It is still important to state that they also do not depend on whether bottom or charm quarks are being produced, and thus exhibit important probes for the proton or nucleus structure.
At this point we turn our attention to a different observable. As seen above, the angularintegrated cross sections discussed above are not very convenient for getting any physics information about the elliptic part. Therefore, instead we would like use the cosine-weighted angle average determined as follows: Roughly speaking, the more positive this observable is, the more P ⊥ and ∆ ⊥ are parallel (or antiparallel); the negative case correlates with perpendicular vectors. An azimuthal angle distribution is relatively easy to measure, besides not being influenced by fragmentation. Also, Eq. 4.1 is a ratio of the cross sections implying largely reduced theoretical and experimental uncertainties. For instance, the luminosity as well as the overall normalization are canceled out in such a ratio. We conclude that this is a very reliable observable. Again, we use exact kinematics, but if we set the quark momenta to be equal P ⊥ , the integration over the differential cross section multiplied by cos 2(φ P − φ ∆ ) lets only terms A · B or C · D survive. As such, the study of this observable is relevant to determine the size of the elliptic contribution, and to determine which kinematic region is more interesting to obtain phenomenological information on the Wigner distribution. The ratio of the cosine-weighted angle integrated over the angular-integrated cross section as a function of P ⊥ is shown in Fig. 13. We show both the bb and cc production cross sections, with fixed ∆ ⊥ = 0.2 GeV. We see that with appropriate choices of P ⊥ the information can be extracted about both B and D structure functions.
It is expected from Figs. 4 and 6 that the contribution from the elliptic term rises as ∆ ⊥ increases. This can be seen better by fixing P ⊥ and varying ∆ ⊥ , as in Fig. 14. Besides, we see an oscillation as expected. What would change if the target were a proton? In Figs. 15-18, we show this case as well. Of course the cross sections are smaller since in the heavy ion case we did not divide by A. The dependence on P ⊥ is pretty much the same as in the nuclear case. However, the dependence on ∆ ⊥ does not show oscillations for the ranges studied. For instance, the cosine-weighted average increases steadily with ∆ ⊥ as is seen in Fig. 18.

V. CONCLUSIONS
In this work we have analyzed exclusive heavy quark photoproduction in the forward region in pA and AA UPCs with quasi-real WW photons as a mean to constrain the gluon Wigner distribution. The study of heavy-quark di-jets is important since measurements can be done at lower transverse momenta then for jets emerging from the light quarks. Besides, these observables are less affected by fragmentation effects and have a cleaner QCD background.
We introduced two new structure functions C(P ⊥ , ∆ ⊥ ) and D(P ⊥ , ∆ ⊥ ), besides adding mass corrections to the already known A(P ⊥ , ∆ ⊥ ) and B(P ⊥ , ∆ ⊥ ). By using the generalized MV model for lead and proton targets, we studied the non-trivial features of these functions, especially the ones related to the elliptic component of the Wigner distribution, B and D. These are of definite importance to understand the angular correlations between the transverse momenta k ⊥ and ∆ ⊥ in the GTMD, and can be related to elliptic flow in hadron and/or nuclei collisions.
We numerically calculated the bottom and charm pair production cross sections in a fully differential form using the above structure functions, both for proton and nucleus target. We provided results in the forward region, which is accessible by ATLAS and CMS. Regarding the azimuthal-angle integrated cross section, we showed that the new C structure function has a big impact in the low-P ⊥ region (P ⊥ 4 GeV), while A dominates the differential cross section in the range 4 P ⊥ 7 GeV. The C function starts to be relevant again for P ⊥ 7 GeV.
We then defined the cosine-weighted angular average of the differential cross section in order to access the elliptic part of the hadron structure. This observable can help to constrain the elliptic distribution, since it is directly connected to the products A · B and C · D. As neither of the products dominates in the considered kinematic regions, they can be be probed simultaneously by a measurement.
We hope that the observables of this paper will be taken into account by the forward physics experimental groups when planning for new measurements.