Disentangling coherent and incoherent quasielastic $J/\psi$ photoproduction on nuclei by neutron tagging in ultraperipheral ion collisions at the LHC

We consider $J/\psi$ photoproduction in ion--ion ultraperipheral collisions (UPCs) at the LHC and RHIC in the coherent and incoherent quasielastic channels with and without accompanying forward neutron emission and analyze the role of nuclear gluon shadowing at small $x$, $x=10^{-4}-10^{-2}$, in these processes. We find that despite the good agreement between large nuclear gluon shadowing and the ALICE data in the coherent channel, in the incoherent channel, the leading twist approximation predicts the amount of nuclear suppression which is by approximately a factor of $1.5$ exceeds that seen in the data. We hypothesize that part of the discrepancy can be accounted for by the incoherent inelastic process of $J/\psi$ photoproduction with nucleon dissociation. To separate the high-photon-energy and low-photon-energy contributions to the $d \sigma_{AA\to AAJ/\psi}(y)/dy$ cross section, we consider ion--ion UPCs accompanied by neutron emission due to electromagnetic excitation of one or both colliding nuclei. We describe the corresponding PHENIX data and make predictions for the LHC kinematics. In addition, in the incoherent quasielastic case, we show that the separation between the low-photon-energy and high-photon-energy contributions can be efficiently performed by measuring the correlation between the directions of $J/\psi$ and the emitted neutrons.


Introduction
Recently coherent and incoherent photoproduction of J/ψ in ultraperipheral collisions (UPCs) of nuclei was measured by the ALICE collaboration at the LHC [1,2]. In the coherent channel, a large reduction of the coherent cross sectionapproximately by a factor of three-as compared to the impulse approximation has been reported. Such a magnitude of the suppression was found to be in the reasonable agreement with the expectations of the approaches predicting significant nuclear gluon shadowing at x ≈ 10 −3 (x is the fraction of the nucleus momentum carried by gluons), notably with predictions of the leading twist approach to nuclear shadowing [3,4] and with the results of the EPS09 global QCD fit to nuclear parton distributions [5]. Thus, charmonium photoproduction on nuclear targets is a useful tool to study nuclear gluon shadowing at small x.
The aim of this paper is twofold. First, we extend application of the leading twist approach to nuclear shadowing [6] to incoherent quasielastic photoproduction of J/ψ on nuclei and show that the suppression of both coherent and incoherent J/ψ photoproduction in ion-ion UPCs can be described in the same framework. A comparison of the resulting theoretical prediction for the cross section of incoherent J/ψ photoproduction in Pb-Pb UPCs at the LHC to the ALICE data, which is also characterized to correspond to an incoherent quasielastic process [1], shows that the expected suppression due to nuclear shadowing is larger than that seen in the data. We argue that this does not only place additional constraints on models of nuclear shadowing down to x ≈ 10 −4 but also indicates that additional processes can contribute to the ALICE data. In particular, on top of incoherent J/ψ photoproduction on nuclei resulting from the target nucleus excitation, the γ+A → J/ψ + Y + (A − 1) * process driven by the γ + N → J/ψ + Y nucleon dissociation (Y denotes products of the nucleon dissociation) accompanied by the nucleus breakup into the (A − 1) * system consisting of nucleus debris or nucleons also leads to the inelastic final state. The calculation of the γ+A → J/ψ+Y +(A−1) * contribution is rather involved reflecting different mechanisms of the elementary reaction at small and large |t| considered in [7] and will be addressed in a separate publication.
Second, we discuss specifics of and make predictions for coherent and incoherent charmonium production in nucleus-nucleus UPCs accompanied by forward neutron emission which can be studied at the LHC with the ALICE, CMS and AT-LAS detectors equipped by zero degree calorimeters (ZDC). The following channels can be studied: (i) one of the nuclei emits at least one neutron while its partner does not -(0nXn); (ii) both nuclei emit neutrons in opposite directions -(XnXn), (iii) neither of the nuclei emits neutrons -(0n0n). We show that selection of a specific channel can strongly influence the ratio of the cross sections of incoherent to coherent J/ψ photoproduction in Pb-Pb UPCs at the LHC. In particular, we argue that the study of incoherent production of charmonium in ion-ion UPCs with the nucleus breakup allows one to separate the low-photonenergy and high-photon-energy contributions to nuclear J/ψ photoproduction and, hence, to provide additional information on the dynamics of nuclear shadowing of the gluon distribution in nuclei which is complementary to that obtained from coherent onium production. This paper is organized as follows. In Sect. 2, we discuss the suppression of the coherent and incoherent nuclear J/ψ photoproduction cross sections due to nuclear shadowing. We briefly recapitulate main results of the vector meson dominance and the color dipole models for these processes and present the derivation of the coherent σ γA→J/ψA and the incoherent σ γA→J/ψA ′ cross sections in the leading twist approximation. In Sect. 3, using the results of Sect. 2, we make predictions for the coherent and incoherent cross section of J/ψ photoproduction in ion-ion UPCs without and with neutron emission and analyze the obtained results. Section 4 presents a brief summary of the obtained results.
2 Nuclear gluon shadowing in coherent and incoherent J/ψ photoproduction on nuclei 2.1 The coherent nuclear J/ψ photoproduction cross section The earliest model for production of vector mesons off nuclei is the vector meson dominance model based on hadronic degrees of freedom [8]. In the high-energy optical limit of the Glauber model, the standard expression for the cross section of coherent J/ψ photoproduction on a nuclear target reads (see, e.g., [9]): where we assumed that the multiple interactions leading to the nuclear shadowing effect do not distort the shape of the transverse momentum distribution of the vector meson. In Eq.
is the nucleus form factor (its normalization is F A (0) = A) whose square takes into account that the nucleus remains intact in the coherent process, t is the four-momentum transfer squared and t min = −M 4 J/ψ m 2 N /W 4 γp is its minimal value (M J/ψ and m N are the masses of J/ψ and the nucleon, respectively); Wγp is the photon-nucleus center of mass energy per nucleon; σ tot J/ψA and σ tot J/ψN are the total J/ψ-nucleus and J/ψ-nucleon cross sections, respectively. Note that in the first and second terms in Eq. (1), the dependence on t min has been safely neglected compared to the Φ A (t min ) term.
In the optical limit of the Glauber model, the σ tot J/ψA cross section is: where is the nuclear density. Equation (2) describes successive multiple interactions of J/ψ with target nucleons, whose destructive interference results in the nuclear attenuation (shadowing) [10] of the J/ψ-nucleus cross section, σ tot J/ψA < Aσ tot J/ψN . The main issue with Eqs. (1) and (2) is what value of the elementary σ tot J/ψN cross section to use. It has been well known for a long time that if one tries to determine σ tot J/ψN using the vector meson dominance model and the data on the elementary γp → J/ψp process, the obtained value of σ tot J/ψN is small, namely, σ tot J/ψN (Wγp = 5 GeV) ≈ 1 mb and σ tot J/ψN (Wγp = 100 GeV) ≈ 3 mb, see, e.g., [11]. As a result, the effect of nuclear shadowing in the σ VMD γA→J/ψA cross section predicted using Eqs. (1) and (2) turns out to be small in the small-x region, which contradicts the ALICE data. Also, the smallness of σ tot J/ψN serves as an indication that in the strong interaction, J/ψ reveals properties of a small-size dipole built from a heavy quark-antiquark pair.
The simple space-time picture of heavy onium production in the vector meson dominance model is superseded by the one in high-energy QCD, where the process of charmonium photoproduction involves three stages: (i) the photon conversion into a qq component (dipole) long before the target, (ii) the interaction of the dipole with the target, and (iii) the conversion of the qq component into the final state vector meson. This space-time picture is properly taken into account/realized in the framework of the QCD dipole approximation [12,13], where the large mass of the c-quark, mc, sets the hard scale in diffractive charmonium photoproduction, , and thus, leads to the dominance of the dipoles of the small transverse size d t , d 2 t ∼ 1/µ 2 . The corresponding qq dipole-target cross section reads [14]: where αs is the running strong coupling constant; G T (x,Q 2 ) is the gluon distribution in the target T , which depends on the momentum fraction . In general, in the dipole approach, the J/ψ photoproduction cross section involves dipoles of different transverse sizes d t , which corresponds to different scales Q 2 in Eq. (3). However, using the leading logarithmic approximation and making simplifying assumptions about the gluon transverse momentumk t (k 2 t ≪ µ 2 ) and the J/ψ wave function (assuming that the transverse momenta of c-quarks, k t , are small (k 2 t ≪ m 2 c ), while the longitudinal ones are equally shared between the c-quark and antiquark), it was shown that the imaginary part of the J/ψ photoproduction amplitude involving the dipole cross section of Eq. (3) is expressed through the gluon density of the target at the scale of µ 2 = M 2 J/ψ /4 = 2.4 GeV 2 and one obtains [15]: where C(µ 2 ) = (1 + η 2 )R 2 g F 2 (µ 2 )M 3 J/ψ Γeeπ 3 αs 2 (µ 2 )/(48αe.m.µ 8 ); Γee is the width of the J/ψ electronic decay and αe.m. is the fine-structure constant. The factors of η and Rg correct Eq. (4) for the real part and the skewness of the γT → J/ψT scattering amplitude, respectively; the factor of F 2 (µ 2 ) absorbs all effects not included in the used approximation (F 2 (µ 2 ) = 1 in the non-relativistic limit for the J/ψ wave function).
In a more general case [16], (i.e., beyond the k 2 t ≪ m 2 c limit for the J/ψ wave function), there exists a theoretical uncertainty in the value of µ 2 in Eq. (4) which means that one could use a reasonable range of values, e.g., µ 2 = 2.4 − 3.4 GeV 2 . For example, the suitable value of µ 2 can be determined phenomenologically [4] comparing predictions of Eq. (4) for the proton target with the data.
Working in the framework of the color dipole model, one can calculate the γA → J/ψA scattering amplitude by summing multiple rescatterings of a dipole of the fixed size d t on the target nucleons essentially using Eq. (2) and then integrating over d t with the weight given by the photon and J/ψ wave functions [17], see   (3) is small, the resulting nuclear shadowing is also small [18,19], which contradicts the ALICE data, see the discussion in [4]. In this respect, the situation is similar to the VMD case considered above. In other words, in the dipole formalism, when only qq-dipoles are included, nuclear shadowing in the γA → J/ψA scattering amplitude is a higher twist effect [20]. At the same time, the dipole approach describes reasonably well the RHIC UPC data corresponding to lower energies (larger values of x), where the effect of nuclear shadowing is rather small, see Fig. 5 and its discussion in Sect. 3.2.
In contrast to the dipole formalism, one can use the leading twist framework of QCD factorization theorems, which enables one to apply Eq. (4) directly to nuclear targets. Applying Eq. (4) to the nucleus and proton targets, we obtain the following expression for the t-integrated cross section of coherent J/ψ photoproduction on nuclei at high energy: where G A (x, µ 2 ) and G N (x, µ 2 ) are the gluon distributions in a nucleus and the free proton, respectively; dσ pQCD γp→J/ψp (Wγp, t = 0)/dt is the perturbative QCD (pQCD) cross section on the proton calculated using Eq. (4). Thus, the nuclear modification (suppression) of σ LTA γA→J/ψA (Wγp) is given by the factor of R = G A (x, µ 2 )/[AG N (x, µ 2 )] < 1 quantifying the amount of nuclear gluon shadowing at small x.
The leading twist theory of nuclear shadowing [6] is based on the space-time picture of the strong interaction at high energies, the generalization of the Gribov-Glauber theory of nuclear shadowing in soft hadron-nucleus scattering [10,21] to hard processes with nuclei, and the QCD collinear factorization theorems for the total and diffractive cross sections of deep inelastic scattering (DIS). The approach allows one to make predictions for the leading twist shadowing correction to nuclear parton distributions (nPDFs), structure functions and cross sections, which are given as a series in the number of simultaneous interactions with the target nucleons (the multiple scattering series). The structure of each term in the series is unambiguously given by the Gribov-Glauber theory supplemented by Abramovsky-Gribov-Kancheli (AGK) cutting rules [22] and the QCD factorization theorems.
In the graphic form, the multiple scattering series for the γA → J/ψA scattering amplitude in the leading twist theory of nuclear shadowing is shown in Fig. 2 graph a is the impulse approximation, graph b corresponds to double scattering (the simultaneous interaction of the probe with two nucleons of the target), and graph c corresponds to the interaction with three nucleons of the target. The multiple scattering series of Fig. 2 can be summed as follows. The Gribov result on the inelastic shadowing correction in hadron-nucleus scattering can be conveniently implemented using the formalism of cross section fluctuations [23]. In this approach, the interaction of a high-energy projectile with a nucleus is a two-step process. First, long before the target, the projectile fluctuates into different configurations interacting with a hadronic target with different cross sections σ characterized by the distribution over cross sections P (σ). Second, these fluctuations interact with the nucleus. The corresponding cross section is calculated separately for each fluctuation (for individual σ) using the Glauber method and then averaged with P (σ), for details and references, see [6]. In particular, for the γA → J/ψA scattering amplitude, we obtain: where σ N = dσP (σ)σ N . The factor of κ contains the factors associated with the overlap of the photon and J/ψ wave functions; its value is determined by the elementary γp → J/ψp cross section: dσ pQCD γp→J/ψp (t = 0)/dt = κ 2 σ 2 /(16π). The first term in Eq. (6) describes photoproduction of J/ψ on a single nucleon and, hence, is proportional to the number of nucleons A; it is the impulse approximation corresponding to graph a in Fig. 2.
The second term in Eq. (6) corresponds to the simultaneous interaction of the hard probe with two nucleons of the target nucleus and gives the leading contribution to the shadowing correction; this term corresponds to graph b in Fig. 2. According to the Gribov-Glauber theory of nuclear shadowing supplemented by the collinear factorization theorem for hard diffraction in deep inelastic scattering (DIS) [24], this contribution is unambiguously expressed in terms of elementary diffraction, notably, in terms of the diffractive gluon distribution of the proton G D(3) N [25]. The corresponding interaction cross section is σ 2 (x, µ 2 ): where B diff ≈ 6 GeV −2 is the slope of the t dependence of the diffractive cross section; η ≈ 0.17 is the ratio of the real to the imaginary parts of the diffractive scattering amplitude; the diffractive parton distribution G is the nucleon momentum fraction carried by the diffractive exchange presented by a zigzag line in Fig. 2 (M X is the invariant mass of the intermediate diffractive state) and β = x/x I P is the diffractive exchange ("Pomeron") momentum fraction carried by the active parton.
The structure of the interaction with three and more nucleons of the target (graph c in Fig. 2 and higher terms that we do not show) presents extension of that of graph b: in the interaction with N nucleons of the target, the hard probe diffractively scatters off two nucleons of the target and the produced diffractive state rescatters on the remaining N − 2 nucleons, which leads to its attenuation (absorption). In particular, the third term in Eq. (6) corresponds to the simultaneous interaction of the hard probe with three nucleons of the target; its contribution corresponds to graph c in Fig. 2. This contribution cannot in general be expressed only in terms of diffractive distributions of the proton and needs to be modeled. Since the cross section of hard diffraction in ep DIS exhibits the Wγp dependence typical for soft processes, it appears plausible to model the rescattering cross section responsible for the interaction with N ≥ 3 nucleons (the solid circle in graph c in Fig. 2) using the formalism of cross section fluctuations. Exactly this was assumed in Eq. (6); the corresponding effective cross section is where we suppressed the x and µ 2 dependence of σ 3 for brevity. In practice, the σ 3 cross section is calculated using the distribution P (σ) modeled using the dipole formalism or Pπ(σ) of the pion, see details in [6]. This is reasonable at µ 2 ∼ few GeV 2 , where soft physics dominates. For larger values of µ 2 (e.g., in the case of Υ photoproduction), one can use Eq. (8) only as input at the low initial scale for the subsequent Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution to the desired value of µ 2 .
For the interaction with N ≥ 4 nucleons (not shown in Fig. 2), we assume that the effect of cross section fluctuations is the same as for the N = 3 term, i.e., With this input, the multiple scattering series in Eq. (6) can be summed and the result presented in the following compact form: ) is the total hadron-nucleus cross section in the case when the total hadron-nucleon cross section is σ 3 . Note that in Eq. (9) we did not take into account the small real part of the soft scattering amplitude corresponding to the σ 3 cross section, whose numerical effect is small.
Expressing the γA → J/ψA differential cross section in terms of the amplitude of Eq. (9), we obtain: Note that the expression in the square brackets in nothing but the nuclear gluon (10) and (5) are consistent with each other.
It is important to note that unlike the case of the color dipole formalism, the shadowing correction in Eq. (10) (and also in R) is a leading twist quantity determined by the elementary hard diffraction in lepton-proton DIS. In the low nuclear density limit, when the interaction with N ≥ 3 nucleons can be neglected, the shadowing correction is driven by the leading twist σ 2 cross section. At the low values of x, the N ≥ 3 terms also become important; their contributions are also leading twist quantities, which can be summed using the σ 3 cross section.
It is convenient to characterize the suppression of σ LTA γA→J/ψA (Wγp) due to nuclear gluon shadowing in terms of the S LTA coh (Wγp) ratio: where σ IA γA→J/ψA (Wγp) is the γA → J/ψA cross section in the impulse approximation (IA) neglecting all nuclear effects except for coherence: Theoretical predictions for the suppression factor of S LTA coh of Eq. (11) calculated using the leading twist theory of nuclear shadowing [6] and the EPS09 global QCD fits [5] agree well with the nuclear suppression factor obtained in the recent analysis of the ALICE data on coherent J/ψ photoproduction in Pb-Pb UPCs at x ≈ 10 −2 and x ≈ 10 −3 [3,4] (see also Fig. 4).

The incoherent nuclear J/ψ photoproduction cross section
Similarly to the case of coherent J/ψ photoproduction on nuclei considered in Sect. 2.1, the VMD model, the color dipole formalism and the leading twist approximation can be used to calculate the cross section of incoherent nuclear J/ψ photoproduction, σ γA→J/ψA ′ , where A ′ denotes the final nuclear state containing products of the nuclear disintegration (A ′ = A). Using completeness of the A ′ states, in the high-energy optical limit of the Glauber model, the VMD model gives [9]: is the inelastic J/ψ-nucleon cross section; B J/ψ is the slope of the t dependence of the γp → J/ψp scattering amplitude. Note that Eq. (13) is valid at not too small |t| = 0.
Equation (13) has the straightforward and well-known interpretation: the probability of incoherent (quasielastic) photoproduction of a vector meson on a nucleus is given by the product of the probability of elastic scattering on a single nucleon times the probability for the produced vector meson to survive the passage through the nucleus on its way out. Since σ in J/ψN ≈ σ J/ψN is small, the nuclear suppression of σ VMD γA→J/ψA ′ due to nuclear shadowing is also small. In the color dipole formalism, the coherent and incoherent nuclear J/ψ photoproduction cross sections can be calculated on the same footing. The only difference is the order of averaging over the dipole sizes d t : in the coherent case, one first averages the γA → J/ψA amplitude over d t and then squares the result, while in the incoherent case, one first squares the appropriate scattering amplitude and then averages the result over d t , see, e.g., Ref. [19]. Since the relevant dipole cross section is small, similarly to the coherent case considered above, nuclear suppression of incoherent nuclear J/ψ photoproduction cross sections is small [19]. At small x typical for the LHC kinematics, x ∼ 10 −3 and below, the dipole formalism predictions are subject to rather significant theoretical uncertainties due to the choice of the model for the dipole cross section and for the J/ψ wave function. Nevertheless, one can still make the observation that the shadowing suppression of the incoherent cross section of J/ψ photoproduction on nuclei appears to be larger than that for the coherent case. At the same time, in the RHIC kinematics corresponding to much larger values of x, x ≈ 10 −2 , where the effect of nuclear shadowing is small, both the dipole framework and the leading twist approach provide the good description of the PHENIX UPC data [26]. Note also that at x ∼ 10 −2 , numerous versions of the dipole model correspond to the similar color dipole cross section because the models have been fitted to the same data, see, e.g., Refs. [27,28].
In the leading twist formalism, the cross section of incoherent J/ψ photoproduction on nuclei can be readily calculated using the input employed in the coherent case (10). Generalizing the standard expression for the incoherent (quasielastic) nuclear cross section [9] to include cross section fluctuations, we obtain in the high-energy limit: Note that in Eq. (14), averaging over cross section fluctuations should be performed at the amplitude level which explains presence of two integrals over σ and σ ′ . Note also that in Eq. (14) we assumed that the slopes of the t dependence of the soft scattering amplitudes corresponding to the σ 2 and σ 3 cross sections are equal; the numerical effect associated with the inequality of the slopes is negligibly small. Recalling that dσ pQCD γp→J/ψp (t = 0)/dt = κ 2 σ 2 /(16π) and Eqs. (7) and (8), after some algebra, Eq. (14) can be written in the following compact form in terms of the σ 2 and σ 3 cross sections: In the last line of Eq. (15) we used σ in 3 ≈ σ 3 ; we checked that this approximation works with high accuracy at the level of a few percent.
The physical interpretation of Eq. (15) is similar to that of Eqs. (9) and (11): the nuclear suppression factor in the square brackets arises from multiple interactions of the produced diffractive state with nucleons of the target, which are driven by the σ 2 and σ 3 cross sections.
A comparison of Eqs. (15) and (10) shows that in the leading twist approximation (LTA), nuclear suppression in both coherent and incoherent photoproduction is determined by the same quantities: σ 2 and σ 3 [see Eqs. (11) and (16)]. The σ 2 cross section is a model-independent quantity whose magnitude and x dependence are fixed by the experimentally measured diffractive parton distributions, inclusive gluon distributions and DGLAP evolution equations, see Eq. (7). The σ 3 cross section is a model-dependent quantity of the LTA approach, whose value is constrained using the formalism of cross section fluctuations. In general, σ 3 ≥ σ 2 [see Eq. (8)]; the lower limit on the value of σ 3 , σ 3 = σ 2 , corresponds to the upper limit on the predicted nuclear shadowing.
Equation (15) defines the shadowing suppression factor for incoherent nuclear J/ψ photoproduction, S incoh : Note that Eqs. (15) and (16) are valid at not too small |t| = 0. Figure 3 presents our predictions for S LTA incoh as a function of x = M 2 J/ψ /W 2 γp for 208 Pb at µ 2 = 3 GeV 2 . The shaded band represents the theoretical uncertainty of our predictions associated with the range of possible values of the σ 3 cross section [6].
One should note that since both suppression factors of S LTA coh (11) and S LTA incoh (16) are determined by the essentially soft physics, we expect them to be numerically of a a similar magnitude, with S LTA incoh being somewhat smaller than (S LTA coh ) 2 . Indeed, at x = 10 −3 corresponding to y ≈ 0 for Pb-Pb UPCs at

Coherent and incoherent cases
A high energy nucleus-nucleus ultraperipheral collision takes place when the colliding ions pass each other at the distance |b| in the transverse plane (impact parameter) exceeding the sum of the nucleus radii, |b| > (2 − 3)R A , where R A is the nuclear radius (the UPC physics is reviewed in [29]). In this case, the strong interaction between the nuclei is suppressed and they interact electromagnetically via emission of quasi-real photons. Thus, nucleus-nucleus UPCs offer a possibility to probe very high energy photon-nucleus scattering and, in particular, photoproduction of J/ψ on nuclei. The corresponding cross section has the following form: dσ AA→AA ′ J/ψ (y) dy = N γ/A (y)σ γA→J/ψA ′ (y) + N γ/A (−y)σ γA→J/ψA ′ (−y) , (17) where N γ/A (y) = ωdN γ/A (ω)/dω is the photon flux; y = ln(2ω/M J/ψ ) is the J/ψ rapidity, where ω is the photon energy and M J/ψ is the mass of J/ψ; σ γA→J/ψA ′ is the nuclear J/ψ photoproduction cross section (see Sect. 2). Note that Eq. (17) includes both the case of coherent scattering without the nuclear breakup (A ′ = A) and the case of incoherent (quasielastic) scattering when the final nucleus dissociates (A ′ = A).
The photon flux at large impact parameters b = |b| > 2R A emitted by a fastmoving nucleus, N γ/A (ω), can be approximated very well by the following simple semi-classical expression for the flux of equivalent photons produced by a point-like particle with the electric charge Z: where αe.m. is the fine-structure constant; K 0 (X) and K 1 (X) are modified Bessel functions; X = bω/γ L , where γ L is the nucleus Lorentz factor. The strong interactions between the colliding nuclei is suppressed by the requirement that b > 2R A . Experimentally this corresponds to the selection of events with only two leptons from the J/ψ decay and otherwise no charged particles in the whole rapidity range covered by the detector. The presence of two terms in Eq. (17) reflects the fact that each nucleus can radiate a photon as well as can serve as a target. In the case of symmetric UPCs (e.g., in the case of Pb-Pb UPCs at the LHC and Au-Au UPCs at RHIC), at a fixed value of the J/ψ rapidity y = 0, Eq. (17) contains two contributions: one corresponding to the interaction of high-energy photons with a nucleus and another corresponding to the interaction of low-energy photons with a nucleus.
In the coherent case, the separation of these overlapping contributions is not an easy problem since a priory one cannot say in the interaction with which of the two nuclei J/ψ was produced. As a result, the photoproduction cross section σ γA→J/ψA can be unambiguously extracted from the measured dσ AA→AAJ/ψ (y)/dy cross section only in the following two cases: (i) at y = 0, where the two photon energies are equal and, hence, both terms in Eq. (17) contribute equally, and (ii) in the rapidity range where one of the contributions strongly dominates. Exactly this situation is realized in the ALICE experiment [1,2]: the dσ AA→AAJ/ψ (y)/dy cross section was measured (i) for −1 < y < 1, which allowed one to extract σ γA→J/ψA (Wγp) at Wγp ≈ 100 GeV corresponding to the gluon momentum fraction of x ≈ 10 −3 (Wγp is the γ-nucleus center-of-mass energy per nucleon) and (ii) for −4 < y < −2.5, where the low-energy photon contribution dominates (more than 95%), which probes the nuclear gluon density at x ≈ 10 −2 .
As we already explained in the Introduction, the nuclear suppression factor for coherent nuclear J/ψ photoproduction determined from the corresponding UPC cross section measured by the ALICE collaboration [3,4] compares favorably with the theoretical models predicting large nuclear gluon shadowing, notably, with the leading twist approximation (LTA) [6] and with the EPS09 [5] result. This is illustrated in Fig. 4, where the ALICE data on the coherent dσ AA→AAJ/ψ (y)/dy cross section at the central and forward values of the rapidity |y| are compared to the LTA predictions combined with the CTEQ6L1 gluon parameterization [30] at µ 2 = 3 GeV 2 . One can see from from Fig. 4 that the theoretical calculations, which are made using Eqs. (10) and (17), describe the data well. The coherent dσ AA→AAJ/ψ (y)/dy and incoherent dσ AA→AA ′ J/ψ (y)/dy cross sections as functions of the J/ψ rapidity y at √ s = 2.76 GeV. The ALICE data [1,2] is compared to the LTA theoretical predictions; the bands span the uncertainty of the theoretical predictions.
In the same figure, the LTA predictions for the incoherent dσ AA→AA ′ J/ψ (y)/dy cross section made using Eqs. (15) and (17) are compared to the ALICE data point at |y| ≈ 0 [1]. One can see from the comparison that the LTA predicts the magnitude of suppression due to nuclear gluon shadowing exceeding the one seen in the data by approximately a factor of 1.5.
The shaded bands in Fig. 4 represent the dominant theoretical uncertainty of the LTA predictions associated with the uncertainty in the value of the σ 3 cross section, which in turn results in the uncertainty of the LTA predictions for nuclear parton distributions [6]. The uncertainty associated with the value of the hard scale µ 2 , which was studied in [4], is much smaller and has been safely neglected.
Note that in our calculations, we consider quasielastic scattering and do not take into account the incoherent contribution associated with the nucleon dissociation γ + N → J/ψ + Y [31]. We explained in the Introduction that this process could potentially contribute to the inelastic final state and, thus, affect the ALICE extraction of the incoherent dσ AA→AA ′ J/ψ (y)/dy cross section [1] due to the fact that the ALICE detector does not cover the full range of the rapidity y. While the calculation of the γ +A → J/ψ+Y +(A−1) * contribution requires a a separate publication, one can still make several qualitative observations. First, this contribution is expected to have approximately the same A dependence as that in Eq. (15) (it is proportional to A in the impulse approximation). Second, the magnitude of this contribution is expected to be sizable: (dσ γN →J/ψY /dt)/(dσ γN →J/ψN /dt) ≈ 0.15 at t ≈ 0 and increases with an increase of |t| when the dσ γN →J/ψY /dt cross section becomes progressively more important and eventually exceeds that of the elastic γ + N → J/ψ + N process; σ γN →J/ψY /σ γN →J/ψN ≈ 0.8 for the t-integrated cross sections and for M Y < 10 GeV (M Y is the invariant mass of the proton-dissociative system Y ) [32]. It would be desirable to perform an additional analysis of the AL-ICE data [1] by assuming that the γ + N → J/ψ + N and γ + N → J/ψ + Y contributions to incoherent nuclear J/ψ photoproduction have different slopes of the t dependence, which would enable one to experimentally estimate the contribution of the nucleon dissociation and, thus, will enable a direct comparison of the data with predictions of Eq. (15). In addition, it is likely that due to the interaction of the system Y with the nucleus, nucleon dissociation will lead to a larger number of neutrons originating from the nucleus dissociation. Finally, the study of neutron production in the quasielastic γA → J/ψA ′ process at |t| ≥ 1 GeV 2 , where the γ + N → J/ψ + Y contribution dominates, may be interesting for understanding of the formation time in diffraction.

UPCs accompanied by neutron emission
Besides ALICE, the ATLAS and CMS detectors at the LHC are capable to measure UPC production of J/ψ in the −2.5 < y < 2.5 range of rapidity. While for central rapidities, the interpretation of the corresponding measurements is unambiguous, it is difficult to disentangle the high-photon-energy and low-photon-energy contributions to dσ AA→AAJ/ψ (y)/dy for non-central values of y and, thus, to access the small-x region that we are interested in. In particular, according to the estimates of [3,4], for 1.5 < |y| < 2.5 in Pb-Pb UPCs at 2.76 TeV, the dσ AA→AAJ/ψ (y)/dy rapidity distribution in the coherent case is exceedingly dominated (by the factor of four) by the low-photon-energy contribution corresponding to 10 −2 > x > 5×10 −3 of the probed nuclear gluons. The high-photon-energy contribution is suppressed by the much lower photon flux and the larger nuclear gluon shadowing. Hence, it is rather difficult to extract the high energy nuclear coherent J/ψ photoproduction cross section from the UPC data and, hence, to probe the nuclear gluon distribution around x ≈ 10 −4 .
The method to overcome this difficulty was suggested in [33]. It is based on the observation of [34] that coherent photoproduction of vector mesons in heavy ion UPCs can be accompanied by additional photon exchanges which lead to electromagnetic excitation of one or both nuclei with the subsequent neutron emission. These neutrons will have the energy close to that of the colliding beams and can be detected by zero degree calorimeters placed at large distances on both sides of the detectors. With the additional requirement to have in the final state only two muons from the J/ψ decay (in addition to the neutrons) and the large rapidity gap, i.e., by requiring the absence of any other charged particle in the whole range of y covered by the detector system, the strong interaction of the colliding nuclei should be suppressed.
To calculate the cross section of photoproduction of vector mesons in UPCs accompanied by the additional electromagnetic excitation of the colliding nuclei followed by their subsequent neutron emission, we use the model developed in [34]. This approach is justified by the success of the calculations of [35] describing very well the ALICE data on electromagnetic dissociation in Pb-Pb UPCs [36]. The model is based on the assumption that an additional photon exchange does not destroy coherence of the photoproduction process but influences the impact parameter of the ultraperipheral collision. The latter is taken into account by the modification of the flux of the photons participating in photoproduction: where the impact parameter dependent factor of P i (b) takes into account different channels of the nuclear decay by the neutron emission (i = 0n0n, 0nXn, XnXn, . . .). In particular, the 0n0n-channel corresponds to the selection of events without additional electromagnetic dissociation with the nucleus neutron decay; the 0nXn-channel corresponds to one-side excitation with the nucleus neutron decay of only one of the colliding nuclei; the XnXn-channel corresponds to mutual electromagnetic dissociation with both excited nuclei decaying by neutron emission.
To obtain a rough estimate of the size of the effect, each additional photon exchange in Pb-Pb UPCs leads to the suppression of the cross section by the factor of Z 2 α 2 e.m. ≈ 0.3 − 0.4. Based on the assumption that an additional photon exchange influences only the flux of photons, see Eq. (19), one can try to separate the low-energy (ω 1 ) and high-energy (ω 2 ) J/ψ photoproduction contributions in the rapidity spectra measured in heavy ion UPCs. To this end, it is necessary [33] to measure rapidity distributions of vector meson production in any two channels, e.g., 0nXn and XnXn, and using the calculated photon fluxes for these channels, to solve simple equations at a fixed value of the rapidity y. For the coherent case, one then obtains: Since the photon fluxes N 0nXn γ and N XnXn γ can be calculated with good accuracy, measurements of different channels will allow one to study nuclear gluon shadowing in a wide range of x. According to our calculations presented and discussed below, the relative contributions of low-photon-energy and high-photon-energy J/ψ photoproduction in these channels are strongly different: while in the 0nXnchannel they are almost equal, the high-energy photoproduction dominates in the XnXn-channel.
Production of forward neutrons in quasielastic incoherent photoproduction of J/ψ in heavy ion UPCs with the nuclear breakup has been considered in [37]. Since in this case the momentum transfer in elastic J/ψ photoproduction on the nucleon can be as large as |t| = 1 GeV 2 , this target nucleon escaping from the nucleus participates in additional quasielastic rescattering. The average excitation energy of a heavy nucleus in the one-nucleon removal process is about 20−25 MeV, which is much higher than the separation energy of 7 − 8 MeV of one neutron. It was shown in [37] that in incoherent J/ψ photoproduction in heavy ion UPCs, the residual nucleus will decay emitting one or more neutrons with the probability of about 85%. Therefore, imposing the constraint that no neutrons are emitted, i.e., considering the 0n0n-channel, one can almost completely (at the level of 10 − 15%) suppress the incoherent contribution.
As we mentioned in the end of Sect. 3.1, the contribution of nucleon dissociation becomes important/dominant with an increase of the transverse momentum of J/ψ. This process should lead to at least as many neutrons as the quasielastic process. Therefore, the γ + N → J/ψ + N and γ + N → J/ψ + Y contributions should be either separated experimentally or the latter should be included in the theoretical calculation of the γA → J/ψA ′ cross section. The procedure for the extraction of the high-photon-energy contribution that we discuss below involves the use of the different p t dependences of the γ + N → J/ψ + N and γ + N → J/ψ + Y cross sections, which allows one to separate their contributions. We also note in passing that the study of the neutron multiplicity at p t ≥ 0.8 GeV, where the process of nucleon dissociation dominates, would produce for the first time information about the space-time formation of hadrons in the diffractive processes like γ + N → J/ψ + Y .
The recent PHENIX data on J/ψ photoproduction in ultraperipheral Au-Au collisions at 200 GeV accompanied by neutrons detected in both ZDCs [26,38] gives an opportunity to check main assumptions about the sources of forward neutrons in coherent and incoherent J/ψ photoproduction in heavy ion UPCs at high energies. Indeed, since the nuclear gluon shadowing in these processes is small at RHIC energies (the shadowing effect leads to an approximately 20% reduction at the cross section level), the coherent and incoherent photoproduction cross sections can be calculated using either the dipole model or the leading twist approximation; at x ∼ 10 −2 corresponding to the RHIC kinematics, the dipole model and LTA predictions largely converge. In this work, we use a simple version of the dipole model which employs the standard Glauber expressions [Eqs. (1) and (13)] with the elementary σ J/ψN cross section approximated by the phenomenological Golec-Biernat-Wusthoff dipole cross section [39]. The resulting predictions are labeled "GBW+Glauber" in Fig. 5.
It is worth noting here that for the discussed kinematics, the results for the dipole-nucleon cross section obtained in different dipole models are rather close since they are well constrained by the DIS data for these energies [27,28]. The theoretical uncertainty is much smaller than the PHENIX experimental errors and, hence, it is not shown in Fig. 5. Note also that in the discussed model, the nuclear shadowing effect is driven by the σ ccN dipole cross section and, hence, shadowing is suppressed (it is a higher twist effect) for the dipoles of such a small size, see the discussion above.
Selection of events with two-side neutron detection means that in the XnXnchannel, coherent production should take into account the mutual electromagnetic dissociation that requires at least two additional photon exchanges. In contrast, according to the predictions of [37], incoherent production is associated with excitation and neutron decay of only the target nucleus. Hence, in this case, detection of the XnXn-channel requires only one photon exchange to excite the nucleus, which serves as a source of the photon flux in the process of incoherent production. Figure 5 presents a comparison of the PHENIX results for the XnXn-channel with our theoretical calculations using the simple dipole described above. Despite large experimental errors, it seems that the agreement between our calculations and the PHENIX data demonstrated in Fig. 5 justifies our approach to neutron production in UPCs. Coupling this with the good description of coherent J/ψ photoproduction at the LHC makes our approach reasonable for the prediction of coherent and incoherent UPCs accompanied by neutron emission at the LHC. Figure 6 presents the results of our calculation of the rapidity distributions in J/ψ photoproduction in Pb-Pb UPCs with neutron emission in the 0nXn-channel (upper row of panels) and the XnXn-channel (lower row of panels) as functions of the rapidity y at the LHC energy of √ s N N = 2.76 GeV; the two panels on the left correspond to the coherent case and the panels on the right are for the incoherent case. Our theoretical predictions are made using the leading twist approximation (LTA); the shaded areas represent the dominant LTA theoretical uncertainty associated with the uncertainty in the predicted nuclear gluon distributions [6]. The uncertainty associated with the value of the hard scale µ 2 is much smaller and, hence, has been neglected. Each panel contains two sets of curves: the upper shaded area is a sum of the two terms in Eq. (17) and the lower shaded area represents the contribution of J/ψ production by the photon emitted by the nucleus moving with the momentum in the direction of positive J/ψ rapidity (the dashed curves labeled "one side").
We can draw several conclusions from Fig. 6. First, one can see from the two left panels that using the data on the 0nXn and XnXn-channels, one can try to extract the high-photon-energy contribution using Eq. (20). In the range of rapidities of 1.5 < y < 2.5, this will allow one to determine nuclear gluon shadowing in Pb down to x ≈ 10 −4 and at the scale of µ 2 ≈ 3 GeV 2 from coherent J/ψ photoproduction in UPCs.
Second, a comparison of the corresponding upper shaded areas in Fig. 6 shows that we predict that at central rapidities, (i) in the 0nXn-channel, the coherent and incoherent contributions will be practically of the same magnitude and (ii) in the XnXn-channel, the coherent contribution exceeds the incoherent one by approximately a factor of two.
Third, besides the coherent channel, the incoherent cross section can also be used to study nuclear gluon shadowing at small x. Based on the dominance of the strong interaction mechanism of neutron production in incoherent photoproduction, we predict that there is a good opportunity to separate the low-photon-energy and high-photon-energy contributions in the 0nXn-channel (see the upper right panel in Fig. 6). Indeed, since in this case neutrons are emitted by the target nucleus, there should by a correlation between the direction of the produced J/ψ and the direction of neutrons in UPC events. In the kinematics of UPCs, the direction of charmonium produced by a high-energy photon and, hence, low-x gluons from the target, is opposite to the direction of the target nucleus and, correspondingly, to the direction of neutrons. Conversely, in the low-photon-energy production, the direction of charmonium coincides with the direction of the target and neutrons. The standard procedure of the separation of coherent and incoherent events consists in the analysis of momentum transfer distributions. In coherent photoproduction, this distribution is dictated by the nuclear form factor squared and presents several distinct peaks at small p t < 200 MeV/c, where p t is the transverse momentum of produced J/ψ and t = −p 2 t . One can see such first three diffractive peaks in Fig. 7. The t dependence of the incoherent nuclear cross section is the same as in the elementary process, i.e., exp[−B J/ψ |t|]. Our predictions for the sum of the coherent and incoherent cross sections of J/ψ photoproduction in Pb-Pb UPCs accompanied by neutron emission in the 0nXn-channel as a function of p t in the rapidity range of 1.5 < y < 2.5 is shown in the upper panel of Fig. 7. As in Fig. 6, our theoretical predictions are based on the LTA; the shaded bands represent the LTA uncertainty in the predicted nuclear gluon distribution. In the figure, the peaks at small p t correspond to the coherent signal. Therefore, the p t < 200 MeV/c cut will effectively reject incoherent events with good accuracy.
As we already discussed above, the incoherent cross section can also probe the small-x nuclear gluon distribution. In particular, in the 0nXn-channel, the directions of J/ψ and the emitted neutrons can be unambiguously correlated with either the low-photon-energy or with the high-photon-energy contributions. At the same time, in the coherent channel, due to the assumption of independence of the processes of coherent photoproduction and electromagnetic dissociation followed by the nucleus neutron decay, it is impossible to assign emitting neutrons to the source or the target, and, hence, the probabilities to detect neutrons in the same and in the opposite direction with respect to the direction of J/ψ should be equal.
These predictions are demonstrated in the middle panel of Fig. 7, where the upper band corresponds to the situation when the direction of J/ψ coincides with that of the neutrons and the lower band is for the opposite directions of J/ψ and the neutrons. The bottom panel shows the ratio of the two curves from the middle panel. Since we identify the events, where the directions of J/ψ and neutrons are opposite, with high-energy photoproduction, there could be two sources of suppression of the corresponding cross section: (i) the falloff of the flux of highenergy photons emitted by the nucleus and (ii) stronger nuclear gluon shadowing of the small-x gluons at x ≈ 10 −4 , which we are mainly interested in.

Conclusions
We considered J/ψ photoproduction in ion-ion UPCs at the LHC and RHIC in the coherent and incoherent channels with and without accompanying forward neutron emission and analyzed the role of nuclear gluon shadowing at small x, x = 10 −4 − 10 −2 , in these processes. We extended the formalism of leading twist nuclear shadowing characterized by large nuclear gluon shadowing to the incoherent σ γA→J/ψA ′ cross section. We found that despite a good agreement between the approaches predicting large nuclear gluon shadowing at x ≈ 10 −3 and the large nuclear suppression factor extracted from the ALICE data on coherent J/ψ photoproduction in Pb-Pb UPCs at the LHC [1,2], in the incoherent channel, the leading twist approximation predicts the amount of nuclear suppression due to gluon shadowing which exceeds that seen in the data by approximately a factor of 1.5 (Fig. 4). We hypothesize that one source of the discrepancy could be the contribution of incoherent nucleon dissociation, γN → J/ψY , which could potentially contribute to the ALICE data [1] and which is not taken into account in our theoretical analysis.
In coherent J/ψ photoproduction in ion-ion UPCs, it is problematic (except for y ≈ 0 and large |y|) to separate the contributions of high-energy and low-energy photons to the dσ AA→AAJ/ψ (y)/dy cross section, which reduces the range of x for the studies of small-x nuclear gluon shadowing. This problem can be circumvented by considering J/ψ photoproduction in ion-ion UPCs accompanied by neutron emission due to electromagnetic excitation of one or both colliding nuclei.
Using the leading twist approximation for nuclear gluon shadowing, we made predictions for coherent and incoherent nuclear J/ψ photoproduction in Pb-Pb UPCs accompanied by neutron emission in various channels at the LHC (Fig. 6). In particular, we discuss the strategy allowing one to separate the low-photonenergy and the high-photon-energy contributions to coherent J/ψ photoproduction performing a joint analysis of the data in the 0nXn and XnXn-channels. This gives an opportunity to shift the study of nuclear gluon shadowing to the lower x region of x ≈ 10 −4 .
In addition, in the incoherent case accompanied by neutron emission, we show that the separation between the low-photon-energy and high-photon-energy contributions can be efficiently performed by measuring the correlation between the directions of J/ψ and the emitted neutrons (Fig. 7).
In the kinematics where nuclear shadowing is small, we showed (Fig. 5) that theoretical predictions based on the dipole model agree with the PHENIX data on J/ψ photoproduction in Au-Au UPCs at