Exciting Ions: a Systematic Treatment of Ultraperipheral Heavy Ion Collisions with Nuclear Breakup

We present an updated theoretical treatment of ultraperipheral collisions (UPCs) of heavy ions, within the SuperChic Monte Carlo generator. This in particular accounts for mutual ion excitation through additional photon exchanges between the colliding ions. This effect occurs frequently in UPCs, and indeed can be (and has been) measured in data through the use of zero degree calorimeter (ZDC) detectors installed in the far forward region. The theoretical approach presented here accounts for the non-trivial and non-negligible impact such ion dissociation has on the measured cross sections and distributions of the produced particles in the central detectors. This builds on previous work, whereby the survival factor probability of no additional inelastic ion-ion scattering due to the strong interaction, and its kinematic impact, are also accounted for within the same overall framework. We compare to data from ATLAS and CMS at the LHC, and STAR at RHIC, and find in general encouraging agreement for a range of observables and ZDC neutron tags, with some room for further improvement, suggesting the inclusion of higher order QED effects and/or tuning of the the $\gamma A \to A^*$ cross section may be desirable. Overall, this gives confidence in the approach considered here and for applications to new phenomena within and beyond the Standard Model.


Introduction
The high energy collision of heavy ions is well established as a testing ground for the strong interaction in extreme regimes, see e.g.[1][2][3] for reviews.However the ions themselves carry significant electric charges Z and hence can act as intense sources of electromagnetic radiation, or in the language of particle physics, photon-photon collisions.This is in particular the case for so-called ultraperipheral collisions (UPCs), where the impact parameter separation of the ions is significantly larger than the range of QCD, with the ions remaining intact after the collision.In such a case an object of interest can be produced by a purely QED interaction, with no additional particle activity present in the detector.
This production mechanism is a key ingredient in the LHC precision and discovery programme, providing a unique probe of physics within and beyond the SM, see e.g.[4] for a review.One topical example is the case of light-by-light (LbyL) scattering, γγ → γγ, for which the first ever direct evidence was presented by ATLAS in the UPC channel [5] and subsequently by CMS [6] (see [7] for the first observation and [8] for further analyses).This is in general sensitive to a range of new physics phenomena, with the possible s-channel contribution from the decay of an axion-like-particle (ALP) being particularly relevant to UPCs; indeed, the tightest available bounds on such states in the 10-100 GeV mass region have been placed by both the ATLAS [8] and CMS [6] collaborations (see [9] for a recent study in pp interactions).Even within the SM, possible sensitivity to tetraquark states has been discussed in [10].A further topical case is the potential to use the UPC channel as a probe of the anomalous magnetic moment of the τ lepton via γγ → τ τ process [11,12].Recent observations of this process by both ATLAS [13] and CMS [14] have already placed competitive constraints on this, with further data to come.
The theoretical framework required to model such UPC processes has been the focus of much study (see e.g.[15] for an early and [16] for a recent review) and a range of Monte-Carlo event generators are available [17][18][19][20].However, while the basic ingredients are universal to these, the specific implementation (as well processes considered) is not.A particularly important element is the 'survival factor' probability that the colliding ions do not interact strongly, leading to a high multiplicity event more typical of standard central ion-ion collisions, where the QED production mechanism is not so straightforwardly isolated.Now, provided the initiating photons are coherently emitted from the ions, their typical Q 2 is very low, which precisely corresponds to the larger impact parameter separations where QCD interactions are suppressed.However, they are not negligible, and a precise account of this is mandatory.The survival factor in particular depends on the specific process under consideration as well as the kinematics of the produced object.The first and so far only complete treatment of this for the case of UPCs was presented in [18] (see also [21]) and is implemented in the SuperChic Monte Carlo (MC) generator [22].
A complete treatment of the survival factor, as well as all other elements of the theoretical calculation, is crucial to taking advantage of this UPC process, where a key motivation depends on a precise understanding of the underlying photon-initiated production process.A useful testing ground for this is the production of electron and muon pairs in UPCs, which has been measured at the LHC [23][24][25][26] and RHIC [27].The ATLAS data [24,25] are in particular compared to the results of SuperChic, with generally good agreement seen, but with some systematic differences observed.This is discussed in detail in [21], where the potential importance of higher order QED corrections is emphasised.
An important process of this type is due to additional photon exchanges between the colliding ions, whereby one or both ions can be excited into a higher energy state that subsequently decays by emitting a single or multiple neutrons.The dominant such excited state is known as giant dipole resonance (GDR), although higher excitations are also present, see e.g.[16,[28][29][30][31][32].Thus, in the presence of such exchanges the ions no longer remain intact, although as there is no colour exchange between the ions and the impact parameter can remain large, the event signature in the central detector remains the same.Moreover, these exchanges can be assumed to factorize from the underlying γγ → X process, and hence the cross section and distributions of the UPC process are insensitive to this, provided these are inclusive with respect to such ion dissociation.Nonetheless, this is not always true, and one can in particular measure these ion dissociation processes via zero degree calorimeter (ZDC) detectors, which have been used in UPC measurements at ATLAS [24,25], CMS [26] and STAR [27].These are designed to detect neutral particles produced at very forward rapidities and so can tag events with or without additional neutron production.As ion excitation typically results in the emission of relatively low momenta neutrons in the ion rest frame, these are heavily boosted in the lab frame, and hence such ion dissociation can be measured by tagging these neutrons in the ZDCs.In other words, this is an observable effect that should be modelled appropriately.
An important element of modelling this process is, as with the case of the survival factor, to fully account for the appropriate process and kinematic dependence.The treatment of the latter effect has in particular been found to be incomplete in the Starlight MC generator [32] implementation by the CMS data [26] on the impact of ion dissociation on the acoplanarity of the dimuon system in UPCs.Theoretical calculations that do account for this effect [33] on the other hand describe the trend of the CMS data (see e.g.[34][35][36][37][38] for other studies), but without providing a full MC treatment.In this paper, we therefore extend the SuperChic MC generator to account for mutual ion dissociation, including its full process and kinematic dependence, as well as the interplay with the survival factor probability.We compare to data from ATLAS, CMS and STAR, and find overall rather good agreement, as well as general consistency with the theoretical results of e.g.[33].Therefore, this extension will allow predictions for UPC production to be compared more directly with what is measured experimentally, given as we will see the ion dissociation probability is certainly non-negligible and does impact on the corresponding kinematic distributions in the presence of ZDC tags.
The outline of this paper is as follows.In Section 2.1 we present the key ingredients in the model we apply for UPCs.In Section 2.2 we discuss the ion-ion survival factor, and how this is accounted for.In Section 2.3 we describe how this framework can account for the effect of mutual ion dissociation.In Section 3 we compare to a range of data from ATLAS, CMS and STAR.Finally, in Section 4 we conclude.

Key ingredients
The basic formalism follows that described in for example [18].Omitting the survival factor for now, the cross section can be written as where x i and q i⊥ are the photon momentum fractions (see [39] for precise definitions) with respect to the parent ion beams and the photon transverse momenta, respectively.The photons have momenta q 1,2 , with q 2 1,2 = −Q 2 1,2 , and we consider the production of a system of 4-momentum k = q 1 + q 2 = N j=1 k j of N particles, where dΓ = N j=1 d 3 k j /2E j (2π) 3 is the standard phase space volume.β is as defined in [39] and s is the ion-ion squared c.m.s.energy.
In (1), T is the process amplitude, and is given by where V µν is the γ * γ * → X vertex, i.e. the amplitude that in the on-shell case would couple to the photon polarization vectors ϵ.The normalization factors are where F 2 p (Q 2 ) is the squared form factor of the ion and The squared form factor is given in terms of the proton density in the ion, ρ p (r), which is well described by the Woods-Saxon distribution [40] where the skin thickness d ∼ 0.5 − 0.6 fm, depending on the ion, and the radius R ∼ A 1/3 .To be precise, for Pb ions we take the experimentally determined values [41] R p = 6.680 fm , d p = 0.447 fm , R n = (6.67 ± 0.03) fm , d n = (0.55 ± 0.01) fm .
For Au ions we scale these values of R n,p by A 1/3 , while keeping the d p,n fixed.The density ρ 0 is set by requiring that The ion form factor is then simply given by the Fourier transform in the rest frame of the ion; in this case we have q 2 = Q 2 , so that written covariantly this corresponds to the F p (Q 2 ) which appears in (3).
In (3) G E corresponds to the form factor of the protons within the ion; numerically this has a negligible impact, as the ion form factor falls much more steeply, however we include this for completeness.To be precise, G E is proton EM 'Sachs' form factor, given by in the dipole approximation.In this work we do not use the dipole approximation but rather, as in [42], the fit from the A1 collaboration [43].We note that in (2) we can see that the amplitude depends in general on a range of other kinematic variables (the photon momentum fractions x i , the kinematics of the produced final state and so on).These are always implied, but omitted in the arguments of T for brevity.We keep the dependence on the photon transverse momentum q i ⊥ explicit for clarity, in particular when discussing the survival factor in the following section.
In the high energy limit, and neglecting the off-shellness of the initial-state photons in the γγ → X process, the above expression reduces to a well known result of the equivalent photon approximation [44], which is formulated purely at the cross section level.To be precise, we can rewrite (1) as in these limits.The flux n(x i ) is given by where Equivalently the above formulation follows simply from the structure function expression in e.g.[45], once we drop the subleading magnetic form factor of the ion, which is suppressed by a factor of ∼ Z, and take the high energy limit.However, formulating things as in (1) allows us to fully keep track of the non-zero photon off-shellness, while working in terms of the amplitude T allows us to include the survival factor, and ion dissociation effects, in an appropriate way.We consider this in the next section.Finally, we note that, in the EPA limit it was shown in [46] that (2) can be recast via where M ±± corresponds to the γ(±)γ(±) → X helicity amplitudes and J P z the spin-parity of the corresponding γγ configuration.While we do not explicitly make use of this formula in our calculation, it will be useful to consider the above expression later on when discussing the comparison to the STAR data [27].

The Survival Factor
The inclusion of the survival factor closely follows the description in [18], although here we present this in a somewhat different way, in order to facilitate the discussion when we consider mutual ion dissociation.In the high energy limit, and neglecting the off-shellness of the initialstate photons in the γγ → X process kinematics we can write although we emphasise that in all calculations we make use of the full form as per (1).Our discussion of the survival factor then concerns the object The physical interpretation is clearest when we move to impact parameter space, i.e. taking the Fourier transform of the amplitude that appears above to give The survival factor is then accounted for by considering where Γ A 1 A 2 represents the probability that no inelastic scattering occurs at impact parameter b ⊥ = |b 1⊥ + b 2⊥ |, and weights the cross section including the survival factor in the appropriate way.It is typically written in terms of the ion-ion opacity This is given in terms of the opacity due to nucleon-nucleon interactions, Ω nn , which is in turn given by a convolution of the nucleon-nucleon scattering amplitude A nn and the transverse nucleon densities T n , see [18] for a more detailed discussion.To first approximation, we have i.e. it corresponds to a requiring that the ions to do not overlap in impact parameter, although in the full calculation the opacity shows some departure from this, see [18].
At the amplitude level (18) simply corresponds to replacing where we will comment on the (lack of) any complex phase at the end of this section.Moving back to transverse momentum space, we are therefore interested in where the 'S 2 ' subscript indicates that the survival factor has now been appropriately accounted for; we then substitute this in (1) to get the cross section.A convenient form for the above amplitude comes from defining in terms of which we have where q ′ 1⊥ = q 1 ⊥ −k ⊥ and q ′ 2⊥ = q 2⊥ +k ⊥ .We can see from (20) that as it stands (23) involves an integral that extends from b ⊥ ∼ 2R A to infinity, which while formally convergent, is numerically unstable.With this in mind we can also define in terms of which we have While these two formulations are in principle completely equivalent, we can see that (25) now involves an integral from the finite range b ⊥ = 0 to b ⊥ ∼ 2R A , which is numerically stable.We note that it is common to refer to the survival factor, which corresponds to the ratio of the cross section evaluated using T S 2 to that evaluated using T .However, it is clear from (22) that the ratio is dependent on the photon transverse momenta q i ⊥ as well as the particular form of T (q 1 ⊥ , q 2 ⊥ ).
The former dependence implies that the survival factor will modify the predicted kinematic distributions, both those directly dependent on q i ⊥ such as acoplanarity distributions but also through (4) which couples the q i ⊥ dependence to the photon momentum fractions x i (and hence the γγ invariant mass and rapidity).The latter dependence implies that the survival factor is dependent on the particular process under consideration.Further discussion of this can be found in e.g.[21,42] and references therein.We also note that although only UPCs are considered here, one can in principle readily extend the above approach to include collisions at other centralities, that is by suitably modifying Γ AA in (18) to have support for the appropriate b ⊥ range corresponding to the required centrality class.Finally, one might worry in (21) about the uniqueness of this replacement, given one could in principle multiply by an arbitrary b ⊥ dependent complex phase and still give the same integrated cross section (18) evaluated in impact parameter space.This will however in general modify (23) and, one might worry, the corresponding predictions for the differential cross section with respect to the photon transverse momenta.The solution to this comes from noting that the individual photon transverse momenta are not observable, but only their vector sum, i.e. the transverse momentum of the measured centrally produced state.We can then write (22) as where we have defined b ′ ⊥ = b 1⊥ − b 2⊥ , with b ⊥ defined as before, and It then follows straightforwardly from (28) that and hence any b ⊥ dependent phase in (21) will not contribute.This can therefore without loss of generality be simply dropped.

Including Mutual Ion Dissociation
We include ion dissociation following the well-established formalism described in e.g.[17,32,47,48].The basic idea is that the ion dissociation occurs via a single (or multiple) photon exchanges between the ions, such that coherent photon emission process on one ion side leads to ion dissociation via γA → A * on the other ion side.If the ion dissociates, it will emit some number n neutrons as part of its 'de-excitation'.As these are typically emitted with rather low momenta in the ion rest frame, in the lab frame these are highly boosted to forward rapidities, where they can be detected by ZDCs.The dominant nuclear excitation is a so-called 'giant dipole resonance' (GDR) [16,[28][29][30][31][32], which is the lowest energy nuclear excitation and typically leads to the emission of a single neutron.However, higher excitations are possible, for which a larger number X neutrons will be emitted (in addition to other particles such as pions).We therefore consider three possibilities in our calculation, namely the emission of zero, one or any number X > 0 neutrons; while in principle one can provide predictions for any exclusive number of neutron emissions, we will for simplicity not consider these here, beyond the single neutron case.
We also note that we will only consider mutual ion dissociation here.In general, the initial state photons in the γγ initiated process may be emitted inelastically from the ions, however this is subleading, being suppressed by a factor of ∼ Z, so that it is relevant at higher photon transverse momentum (i.e.p ll ⊥ in the dilepton case), where the coherent production mechanism is suppressed by the elastic ion form factor.At higher transverse momentum, the dominant emission process is due to incoherent inelastic emission from the nucleons in the ion.Therefore, to good approximation, one can follow a technique as in [25] to model this process, namely by using the SuperChic distribution from dissociative lepton pair production in pp collisions and setting the normalization in a a data-driven way; this is found to work rather well in the ATLAS analysis.A related excitation mechanism is for both ions to be excited by a single photon exchange in addition to the γγ initiated process, i.e. for the additional photon exchange not to be due to coherent emission from either ion.This is however suppressed for the same reason, being subleading in Z, and indeed as it is relevant at larger exchange virtualities, contributes in the region of smaller ion-ion impact parameters, which are strongly suppressed by the ion-ion survival factor.Other higher order QED processes that we will not consider in detail here are when additional photons are exchanged between the production process (i.e. the dilepton system in the case of lepton pair production) and the ions, as well as the production of additional low momenta electron pairs, i.e. so-called unitarising effects.These are discussed further in [21] and the references therein; while these will dominantly not lead to ion dissociation, they may effect the cross section predictions at the precision level.A more detailed account of these is left to future work.
Returning to the case of mutual ion dissociation we will consider here, the nuclear excitation is assumed to occur independently of the γγ → X production process, and so in impact parameter space we can write where P X 1 X 2 is the breakup probability, such that X i = 0, 1, X corresponds to the the emission of 0 (i.e.no nuclear excitation) 1, or X > 0 neutrons emitted for each ion i = 1, 2. The probability P X 1 X 2 (b ⊥ ) depends on the total ion-ion impact parameter b ⊥ = |b 1⊥ + b 2⊥ | and factorizes into independent breakup probabilities for each ion, i.e.
In the same way as the γγ → X cross section (10) in the EPA, formulated in impact parameter space, the lowest order breakup probability for each ion i = 1, 2 is then given by a convolution of the photon emission flux from the ion j = 2, 1 and the γA → A * cross section: where ω is the photon energy in the A rest frame, i.e. ω = xs/(2m A ) in the s ≫ m 2 A limit, which holds to very good approximation; the dissociation probability therefore depends on s, as well as the ion beam type, but we omit these arguments for brevity.The flux factor | Ñ | 2 is built from exactly the same ingredients that enter the γγ → X cross section (1).To be precise it is given in terms of the Fourier transform where N(x i ) is defined in (12).σ γA→A * (ω) is the photon-ion excitation cross section, which has been measured over a wide range of photon energies from fixed target ion scattering experiments.The corresponding data points, which we use to perform the integral (34) are show in Fig. 1; these are in many cases as implemented in the Starlight MC [17] (see also [29,48]), with some exceptions that we will describe below.The cross section is shown in Fig. 1 (left), while the result weighted by the flux factor as per (34) is shown in Fig. 1 (right) for the, as we will see, representative value of b ⊥ = 3R A .The largest peak in the ω ∼ 10 − 20 MeV region corresponds to the dominant GDR resonance, data for which is taken from [49] for both the single and multiple neutron emission cases, fitted to a Lorentz shape.In the single neutron emission case the upper limit of the integral (34) is set to the upper limit of these data, at ω = 23.4MeV, where indeed the cross section is negligible.Above this energy, photon absorption will very dominantly lead to multiple neutron emission.In the 23.4 < ω < 440 MeV region we use the data from [50,51] on photonuclear scattering.In the 440 < ω < 2000 MeV region we use the data from [52] on photonuclear scattering; this is in contrast to [17] (see also [29,48]) where older data [53,54] on γp and γn scattering were  used, scaled to the nuclear case but with no shadowing applied.In the 2 < ω < 16. 4 GeV region, we use data from [55,56] on γP b and γAu scattering, suitably interleaved to cover the full energy region and with appropriate A rescaling applied.In our MC implementation, for the GDR region data for both Au and Pb beams are implemented, with A rescaling applied in other cases.For energies above this an effective per nucleon cross section is derived as described above (i.e. using a selection of Au and Pb data), which can be scaled to the appropriate ion A. Finally, above 16.4GeV there is some limited direct photon-ion data in the analysis of [57], in particular in the ∼ 40 − 80 GeV region for Pb ions.To cover the entire energy region, including energies beyond this, we use the fact that the γ-nucleon cross section is observed to obey approximate Regge scaling as in [58].The high-energy data for this is limited to a handful of measurements at HERA, and for concreteness we take the ZEUS extraction [59].The γ-nucleon cross section is then suitably scaled by A, but with a nuclear shadowing factor of 0.65 applied, in order to match the direct data from [57] in the Pb case.
The contribution from this region, while naively suppressed by ∼ 1/ω, as evidenced in Fig. 1 (right), is in fact in principle rather significant, due to the large photon energies available.In particular, from the considerations in Section 2.1 the photon x is cutoff at roughly x max ∼ 1/(R A m A ), and hence the photon ω max ∼ s/(m 2 A R A ).At the LHC this is of order 850 TeV, which corresponds to a γn c.m.s.energy that is roughly an order of magnitude larger than that probed at HERA, and many orders of magnitude higher than the highest energy direct data on photon-ion absorption.To give an estimate, we can assume for simplicity that the flux term in (34) is constant up to the cutoff x max , and that the γA cross section is constant with energy, which is roughly consistent with the Regge parameterisation.We then have if we take the measured value from [57] at ω 0 = 80 GeV.The (in theory dominant) contribution from the GDR region can be estimated using the TKR sum rule (see e.g.[60]) for the Pb case.At the LHC, a more precise numerical evaluation gives 260 mb for the integrated single neutron emission cross section (equal to the total GDR contribution to reasonable approximation) and 177 mb for (36) at a representative value of b ⊥ = 3R A , consistent with Figure 2: Lowest order breakup probabilities for PbPb collisions at √ snn = 5.02 TeV, for single and multiple neutron emission, as defined in (34), given by the dashed curves.Also shown in solid is the result applying the unitarising corrections discussed in the text.
these rough expectations.Including the contribution below 80 GeV and above the GDR region, we find that (36) accounts for roughly 25% of the total contribution to (34) at the LHC, again for b ⊥ = 3R A .Given the lack of direct data in this region, we can conservatively assign an uncertainty of this order in the corresponding neutron tagged cross section.We note that a useful way to bypass this uncertainty source is to consider single neutron tag data (i.e.0n1n and 1n1n), where the contribution from this region is negligible.Indeed, at these high energies, the photon-ion interaction cannot be viewed as a purely electromagnetic one.In particular, according to the rather well established vector meson dominance model [61][62][63] we can view the photon as a superposition of light vector mesons, which then undergoes a hadronic (dominantly inelastic) interaction with a nucleon within the ion.Indeed, for sufficiently high photon energies the γn system will become relatively more central in rapidity, while the underlying γn interaction will be a relatively high multiplicity inelastic hadronic event.Hence it is arguably possible that the produced neutrons may not be detected in the ZDCs or, more significantly, whether some of the products of the inelastic γn interaction will be seen centrally and hence fail the experimental veto on additional particle production in the central detector.Given it has been observed in e.g. the ATLAS data on muon pair production [24] in PbPb UPCs that the Starlight [17] predictions (which include this high energy contribution up to the kinematic limit) tend to overshoot the observed 0nXn and XnXn fraction, we will to be precise cut the cross section off at ω = 500 GeV, which corresponds to |y γn | ∼ 5 for PbPb collisions at √ s nn = 5.02 GeV.As particle production will in general occur at rapidities lower than this, we may expect this to fail the experimental veto, although a precise evaluation would require that we account for the particle multiplicity distribution and the particular experimental cuts; for e.g. the ATLAS analyses [24,25] we will consider in the following section, the relevant requirements are for no additional tracks with p ⊥ > 100 -200 MeV out to |η| = 2.5 − 3.75.Even absent this we can view this as an effective cutoff, driven by the comparison to data, and given the uncertainties in the calculation discussed above.For b ⊥ = 3R A this removes roughly 15% of the contribution.As we will see when comparing to the ATLAS measurement of dilepton production in UPCs, it may be that a more stringent cut is required to match the data.
The resulting lowest order breakup probabilities, for single and multiple neutron emission, are shown in Fig. 2 by the dashed lines.The multiple neutron emission probability is moderately larger than the single emission probability, as expected from the discussion above.However, we can also see that at smaller impact parameters, the probability increases, and in both cases rises above unity.This is driven by the flux in (34) which one can readily show scales as for the dominant region of x (i.e.x ≪ 1) that contributes to (34) and for b ⊥ > R A such that the ion can be treated as a point-like charge.That is the flux is peaked towards low b ⊥ , until we reach b ⊥ < R A where the extended nature of the ion comes into play and the flux becomes suppressed by the ion form factor (although the contribution from this region is in any case negligible once one accounts for the ion-ion survival factor).This indicates an inadequacy in the lowest order (in Z 2 α) perturbative approximation for the ion excitation process, in particular given the Z 2 enhanced flux for photon emission from the spectator ion.To account for this, we follow the approach described in [32], namely assuming each ion excitation process happens independently we have that the number of excitations follows a Poissonian probability, such that The impact of this unitarising is shown in Fig. 2, where we can see this by construction gives dissociation probabilities that never exceed unity.While the dominant effect of this is in fact below the b ⊥ ∼ 2R A region, which is in any case removed once the ion-ion survival factor is included, it is not negligible; it reduces the XnXn cross section by a factor of ∼ 1.5 − 2, depending on the precise kinematics.The corresponding dissociation probabilities (33) are shown in Fig. 3.In the left plot the results prior to multiplying by the probability exp(−Ω A 1 A 2 (s, b ⊥ )) of no ion-ion inelastic scattering are shown for reference, and in the right plot this is accounted for.As expected from the discussion above, the probability for ion dissociation is peaked towards b ⊥ ∼ 2R A by the scaling of the flux (38), before being cut off by the no ion-ion inelastic scattering probability, which rapidly approaches zero below b ⊥ ≲ 2R A .Conversely, if no ion dissociation is required, as in the 0n0n case, than the probability become increasingly close to unity as b ⊥ increases.For the mixed 0nXn case the peaking towards b ⊥ ∼ 2R A is again present, although this is less strong than for XnXn; the 0n1n result is not shown for clarity, though it follows a similar trend.
So far, we have worked purely in impact parameter space, however for a full treatment, and in particular to account for the process and kinematic dependence of the ion dissociation probability, we must translate these to transverse momentum space.To do this, we simply replace and then use this as in (22) to get the corresponding amplitude in transverse momentum space and hence cross section.More precisely, from (23) we are interested in the integral where J 0 is a Bessel function of the first kind.Making use of (38) we have where The large b ⊥ limit of ( 41), for which we have Γ A 1 A 2 (b ⊥ ) ∼ 1, is therefore driven by the integrand where in this limit we have such that which is numerically rapidly converging1 .For the same reason, once we appropriately use (25) rather than (23) then the 0n0n results in an integral that has the same convergence.For the remaining cases however we have and the numerical convergence of the integral (though it is certainly finite) is more problematic.
To resolve this, we can consider the integrands The first (square bracketed) term now to good approximation scales as ∼ J 0 (b ⊥ k ⊥ )/b 2 ⊥ and is hence numerically convergent, while the latter can be evaluated as 2πA where we have used while we can see that the integrand of the first term in (49) only has support for b ⊥ ∼ 2R A and hence can readily be evaluated.The 0nXn case can be evaluated in a similar way, or alternatively by rearranging where 'inc.' denotes the inclusive case, with no ion dissociation requirement applied, and for the 0nXn it is always implied that either ion can dissociate (but not both).
Having established the framework to account for ion dissociation in transverse momentum space, in the following sections we will now consider results of this and compare them to data.We note that, as discussed at the end of Section 2.2 for the case of the ion-ion survival factor, the impact of the ion dissociation requirement (i.e. the predicted cross section for ion dissociation) will depend on the kinematics of the central system as well as the corresponding process under consideration.This is clear from Fig. 3, where we can see that the impact parameter dependence of the dissociation probability is different depending on the particular neutron tag.This will therefore lead to a non-trivial dependence in the Fourier conjugate transverse momentum space, as we will see.

Comparison to ATLAS data
In this section we will compare our predictions to ATLAS measurements of muon [24] and electron [25] pair production in ultraperipheral PbPb collisions at √ s nn = 5.02 TeV.We will focus throughout on those observables that are presented explicitly with a neutron tag imposed in the ZDCs, although in [24,25] the inclusive data are compared to SuperChic predictions, with the agreement in general found to be good.We note that by 'tag' it is always implied that neutrons are either required or vetoed on, while by inclusive we mean that both cases are included, i.e. no ZDC requirement is made, although we still require the collisions to be ultraperipheral, with no colour flow between the ions.In the ATLAS case, three event categories are considered, namely the case where no neutrons are registered in the ZDCs, where X > 0 neutrons are registered on one side but none on the other, and where X > 0 neutrons are registered on both sides, denoted by 0n0n, 0nXn and XnXn, respectively.The cross section fractions f injn are then defined relative to the inclusive case, such that We begin by comparing to the data on muon pair production [24].In Fig. 4 we show the 0nXn and XnXn fractions as a function of the dimuon invariant mass and for different dimuon rapidity regions; the remaining 0n0n fraction is then found using (52) and so is not independent of these.In terms of the broad trend, we can see that both of these fractions are predicted to increase with m µµ , which is as we would expect.In particular, from (4) we can see that the average initial-state photon transverse momentum will increase with the photon momentum fraction x i ∝ m µµ .In impact parameter space this corresponds to smaller b ⊥ values, and we can see from Fig. 3 that it is precisely the 0nXn and XnXn cases that are enhanced in this region, with the 0n0n, conversely, suppressed.This trend is clearly observed in the data, even within the relatively large uncertainties at larger m µµ .In more detail, however, we can see that there is a distinct tendency to overshoot the data in the lowest mass bin (and at central rapidities, in the second mass bin).That is, too much ion dissociation is predicted, for all dimuon rapidities.
In Fig. 5 we show a comparison to the same dataset, but now presented differentially in the dimuon rapidity, for different dimuon invariant mass regions.We can see that the 0nXn and XnXn fractions are predicted to increase with rapidity, which is again driven by the changing photon q ⊥ dependence in the production amplitude due to (4), and the dependence of the photon momentum fractions on the dimuon rapidity.In this case, forward rapidity corresponds to an increased photon x i on one side, but a decrease on the other, and so it is not immediately obvious that the overall trend should be for larger average photon transverse momenta, and hence smaller impact parameters b ⊥ .This is however the basic trend predicted by the full numerical treatment, and indeed similar behaviour is predicted in earlier studies of photon-initiated production in heavy ion and pp collisions, see e.g.[18,42].We can see from the plots that this trend is indeed observed in the data, again supporting the theoretical framework presented here.On the other hand, the same overshoot in the lowest m µµ bin is clear.
We next consider a comparison to the ATLAS data on electron pair production [25].In Fig. 6 we show the 0n0n, 0nXn and XnXn fractions as a function of the dielectron invariant mass, for different dielectron rapidity regions.As data for the three cases are provided by ATLAS, we present comparisons for all of these, however as noted above only two of the three are independent, i.e. the sum in any given bin is by construction unity.We note that the event selection, given in the figure captions, is rather similar between this and the muon measurement, with in this case somewhat lower invariant masses being probed.The same overall trend as predicted above is seen with respect to the pair invariant mass, m ee , i.e. for the 0nXn and XnXn fractions to increase, and the 0n0n to decrease, as m ee increases.This is again clearly observed in the data, and in general the level of agreement between data and theory is good.The most visible differences are in the forward rapidity, 1.8 < |y ee | < 2.4, bin, both at low and high masses.In the high mass bin, however, the data errors are rather large and certainly the rather extreme suppression in the 0n0n case is not seen in the other rapidity bins.
Looking more closely, we can see that there is again a general trend to overshoot the 0nXn, XnXn data, and undershoot the 0n0n data, as there was in the dimuon data.Given, as discussed in Section 2.3, there is a reasonable theoretical uncertainty in the predicted γA → A * cross section, it is interesting to investigate how much lower the input cross section would need to be in order to better match the data.Before doing so, we note that the second invariant mass bin, 10 < m ee < 20 GeV, in the dielectron measurement covers the same region as the lowest invariant mass bin in the dimuon measurement, see Figs. 4 and 5, and for the same dilepton rapidity region; the only difference from the point of view of the kinematic cuts is the tighter p ⊥ cut in the dimuon case.We would therefore expect rather similar fractions f in both cases, and indeed that is true to very good approximation in the theoretical predictions.In terms of the data, on the other hand, the 0nXn and XnXn fractions are rather higher in the dimuon case, at the ∼ 2σ level.In Fig. 7 we therefore show comparisons to both the dielectron and  dimuon invariant mass distributions on the same plot.We can see that indeed the predicted 0nXn, XnXn distributions overshoot both sets of data, in particular at lower mass, but that this occurs rather more significantly for the dimuon data.We also show in the dashed histogram the predicted fractions when the default γA → A * cross section is multiplied by a factor of 0.8, and can see that in this case the agreement is rather better.It may therefore be that some amount of tuning is required in the future to better match the data.To enable this, in the public SuperChic 4.2 release we provide a flag (fracsigX) by which the normalization of the γA → A * cross section may be modified.We note however, that in practice a reduction in the γA → A * cross section cannot simply be achieved by e.g.removing the higher energy and less well constrained region, where a Regge theory parameterisation must be used; even removing the entire cross section above ω > 20 GeV (corresponding to |y γn | ∼ 6.5) only reduces the cross section by a further ∼ 10%.This therefore corresponds to a fairly significant reduction.We note that in principle another variable that will impact on the predicted dissociation fractions is the treatment of the survival factor.For example, if the suppression due to this is increased and/or pushed to lower impact parameter values, this will modify the average impact parameter sampled in the cross section.Given the dissociation probabilities have distinct impact parameter dependencies as in Fig. 3, this will then modify these.However, on closer investigation we find that it is only with the rather extreme variations in the survival factor (of the type examined in [21]) that a noticeable reduction in the 0nXn, XnXn fractions occurs, and not necessarily with a better description of both cases simultaneously.A further way to shed light on this issue would be to present data in the 0n1n and 1n1n channels.In these cases the ion dissociation is guaranteed to be dominated by the GDR region, where uncertainties due to the As in Fig. 4 but with the ATLAS data in the dielectron channel [25] for the same mass bins (and with a very similar event selection) shown.Theoretical predictions correspond to the dimuon event selection, but results for the dielectron case (which is very similar) are barely distinguishable, and hence are not shown for clarity.The solid histograms correspond to the default SuperChic 4.2 predictions, while the dashed curves correspond to the result with the γA → A * cross section (34) multiplied by 0.8, for comparison.high energy regime that enter the Xn0n and XnXn cross sections are absent.It would therefore be useful to determine whether the difference between data and theory persists in these single neutron channels, or is absent.
Finally, in Fig. 8 we compare predictions for a range of kinematic observables in the 0n0n case.Namely, the cosine of the electron scattering angle θ * in the dielectron rest frame, the dielectron invariant mass and rapidity, and the average single electron/positron transverse momentum, p e ⊥ .We can see that in general the agreement between data and theory is very good, both at the level of normalization and shape; the p e ⊥ and m ee distributions match well across roughly four orders of magnitude in the cross section.The largest discrepancy is in the cos θ * at forward angles, which may indicate some sensitivity to higher order QED effects in the γγ → e + e − process that are absent in the current theoretical treatment.In addition, we can see that in terms of the normalization there is a trend for the theoretical prediction to lie at the lower end of the data uncertainty region, though they are consistent within errors.

Comparison to CMS data
In this section we compare to CMS measurement [26] of muon pair production in ultraperipheral PbPb collisions at √ s nn = 5.02 TeV.These are presented for a wide range of neutron tags, namely the 0n0n, 0n1n, 0nXn, 1n1n, Xn1n and XnXn final states.The acoplanarity, α = 1 − ∆ϕ µµ /π, distributions are measured over a wide range of α and fitted with a 'core' and 'tail' distribution, where the former is taken to be dominated by the leading order γγ → µ + µ − process, and the latter is sensitive to higher order QED effects such as photon radiation from the dimuon system.The core distributions are then used to extract average acoplanarities for each neutron tag, as well as the average invariant dimuon invariant mass, comparisons to which as shown in Fig. 9.
In the Fig. 9 (left) the average acoplanarity is shown.We can see that the broad trend that is predicted is for this to increase as one requires more neutron emission, moving from the 0n0n to the XnXn cases.This is exactly as expected from the discussion in the previous section, namely that as in Fig. 3 the 0n0n (XnXn) dissociation probability is peaked least (most) strongly at lower ion-ion impact parameters b ⊥ .Given this we expect the, Fourier conjugate, dimuon transverse momentum p µµ ⊥ to be peaked most strongly towards lower (higher) values in the 0n0n (XnXn) case, with the intermediate cases lying between these extremes.This is precisely as observed in Fig. 9 (left), bearing in mind that the dimuon acoplanarity increases directly with p µµ ⊥ .This trend is also clearly seen in the data, validating this prediction.We note that this result only comes after including a full treatment of ion dissociation in transverse momentum space, as described in Section 2.3.Indeed, in the CMS analysis [26] the results are compared to the Starlight MC [17], where this is absent, and a rather flat scaling is predicted, in contradiction to the data.A similar level of agreement to the current case is on the other hand seen when comparing to the prediction of [33], which applies a similar framework to us, albeit without a full MC implementation.
While the overall trend in the data is consistent with our predictions, we can see that there is a systematic offset between data and theory, with the data having a somewhat larger average acoplanarity across all neutron tags.To understand this better, in Fig. 10 we show the acoplanarity distributions for the two extreme, 0n0n and XnXn, cases.The data and theory predictions are defined such that the cross section is normalized in the α < 0.004 region, where we expect higher order QED effects due to e.g.photon FSR to be less significant, see e.g.Fig. 12 of [25].We can that at lower acoplanarity the data are described rather well, but that beyond α ∼ 0.002 − 0.004 a clear excess in the data becomes apparent, which increases with increasing α.It is therefore possible that the average acoplanarity, as derived by the 'core' distribution fits in the CMS analysis (which from Fig. 1 of [26] describe the data well out to α ∼ 0.004), is being  The distributions are defined such that the cross section is normalized in the α < 0.004 region, where higher order QED effects are less significant.The muons are required to have p ⊥,µ > 3.5 GeV, |ηµ| < 2.4, 8 < mµµ < 60 GeV.Data errors correspond to systematic and statistical added in quadrature.
driven to somewhat larger values by such higher order QED effects.As these are absent in the current theoretical treatment, this would then lead the observed excess in Fig. 9 (left).In Fig. 9 (right) the average dimuon invariant mass is shown.The basic trend is for this to increase as one requires more neutron emission, and is again as expected from the discussion in the previous section, and indeed observed in the ATLAS data, for which the 0nXn and XnXn event fractions are enhanced at larger dilepton invariant masses.This is therefore again an encouraging validation of the overall approach.Some excess of data over theory is on the other hand observed in the 0nXn and XnXn cases, albeit within relatively large experimental errors.

Comparison to STAR data
Finally, in this section we compare to the STAR measurement [27] of ultraperipheral electron pair production in AuAu collisions at √ s nn = 200 GeV.These data are taken with a XnXn neutron tag imposed, or more precisely a Y nY n tag with Y = 1, 2, 3 suitably corrected up to the full X > 0 case.A particular observable of interest is the azimuthal angle ∆ϕ, defined in [27] as the angle between the dielectron transverse momentum, p ee ⊥ , and the transverse momenta of one  of the e ± .It is in particular predicted in [34] (see also [64][65][66]) that a modulation of the type should be observed in the case of dilepton UPCs, with A 2∆ϕ being zero up to lepton mass corrections ∼ m 2 l /m 2 ll , and the precise value of A 4∆ϕ depending on the specific kinematics.To be precise, the variable ∆ϕ is in fact defined in [34] to be the angle between the the dielectron transverse momentum and the vector difference l ⊥ = p e + ⊥ − p e − ⊥ , which only coincides with the STAR experimental definition in the p e ⊥ ≫ p ee ⊥ limit.While this is true to reasonable approximation, it is as we will see not exact, and so we will for completeness consider both definitions, denoting that of [34] as ∆ϕ, and the STAR definition as ∆ϕ ′ .
In Fig. 11 we compare our predicted normalized ∆ϕ distribution to the STAR data, while in Table 1 we compare the predicted values of A (2,4)∆ϕ , as well as the root mean squared dielectron transverse momentum, p ee ⊥ , and the fiducial cross section to the measured values.We first observe that the predicted total cross section is in excellent agreement within uncertainties with the data.Next, in Fig. 11 (left) we also show for comparison the predicted distribution in the inclusive (with respect to neutron tag) case, with and without including the ion-ion survival factor.In general, we observe an oscillatory behaviour with respect to ∆ϕ, which arises due to the photon polarization dependence in the γγ → X production amplitude (see (13)).This effect was first discussed in [67] and is appropriately incorporated in SuperChic, but is not always accounted for in public MCs, such as e.g.Starlight [32].
We can see by comparing to the inclusive case that the XnXn tag requirement and the survival factor both have a non-negligible (and in fact counteracting) impact on the predicted distribution.This is due to the differing impact parameter dependence of the ion dissociation probability with respect to the inclusive case, which modifies the corresponding amplitude (22) (after making the replacement (40)) in a non-trivial way.We can then see from (13) that this amplitude couples the initial-state photon transverse momenta to the γγ → e + e − vertex such that the weight of the contributing photon helicity amplitudes will be modified by the ion dissociation probability (as well as the ion-ion survival factor) and its particular impact parameter dependence.Indeed, a discussion of this effect from a different perspective is presented in [34].As an aside, we note that for the 1n1n case the predicted distribution is very similar , as well as the root mean squared dielectron transverse momentum, p ee ⊥ , and the fiducial cross section to the STAR data [27].The corresponding coefficients with respect to ∆ϕ ′ are very similar, and are therefore not given.Uncertainties correspond to sum in quadrature of all quoted sources from the experimental analysis.
(although not identical) to the XnXn case, hence any potential effect from the fact that the ZDC tag does not extend beyond 3 neutrons should be negligible.
For the appropriate XnXn case we can see in Fig. 11 that the predicted distribution matches the data reasonably well.The agreement is in particular excellent in the ∆ϕ ≳ π/2 region, while below this there are some discrepancies.In the right plot, however we also show the result of a direct fit using (53) to the data, and while this achieves a somewhat better agreement at low ∆ϕ, overall the agreement is not significantly improved.We also show a comparison to the experimental definition, ∆ϕ ′ ; we can that the predicted distribution is indeed mildly different, with a somewhat better agreement at lower ∆ϕ, but again overall the agreement is not significantly improved.
Therefore, there is possibly a limit to how well any prediction can match the measurement, perhaps due to fluctuations in the data or some other systematic effect.Nonetheless, in Table 1 the predicted values of A 2,4∆ϕ are compared to the values determined in [27] from a fit to the data and there is fair agreement, at the ∼ 1 − 2σ level.With further more precise data the agreement may of course improve, but a more precise analysis, accounting for e.g. higher order QED effects may improve this.Indeed, in e.g.[34,35] the impact of accounting for photon FSR from the dilepton system is discussed and found to be non-negligible in some cases.Similarly, in Table 1 the root mean squared dielectron transverse momentum is also given, and again found to be in fair agreement with the data, but to lie below the measured value; photon FSR effects will act precisely to increase this.Indeed, in Fig. 12 we compare the predicted transverse momentum distribution with the STAR data and the agreement is seen to be very good for most of the p ⊥,ee region, but with some excess of data over theory at the larger values, which again is as we would expect from photon FSR effects.We also show for comparison the predicted distribution in the inclusive (with respect to neutron tag) case, with and without including the ion-ion survival factor.The impact of a full kinematic account of ion dissociation is again clear, and it is only after doing this that the data are matched well.The fact that the XnXn tag leads to a broadening of the transverse momentum distribution towards higher values is again exactly as expected from the impact parameter dependence of the XnXn dissociation probability, which is more strongly peaked to lower values, i.e. larger p ⊥,ee .
As mentioned above, in [34,65,66]  dielectron transverse momentum in the calculation.More precisely, the predicted distribution (53) in fact corresponds [65] to the differential cross section with respect to p ee ⊥ and the vector difference l ⊥ defined above, whereas the observed cross section is of course integrated over these.The lepton transverse momenta cuts p e ⊥ > p cut ⊥ = 200 MeV correspond in terms of these to which can introduce a dependence on cos ∆ϕ, that is not captured by (53) with A 2∆ϕ = 0.The precise form of this depends on the above cut and its non-trivial interplay with the full kinematic dependence of the production cross section, and hence is not straightforward to predict analytically.It is only by providing a full MC treatment, as we do here, that this can be accounted for.We note that if we remove the p e ⊥ cuts, then the predicted value of A 2∆ϕ is indeed consistent with zero, as it is if we increase the threshold on m ee to e.g ∼ 4 GeV; this effect is therefore rather specific to the STAR kinematics.If we reduce the electron mass arbitrarily, then this makes a negligible difference, confirming that this is not the relevant factor.
In the STAR analysis [27] the data are compared to a 'QED' prediction from [68].Although the language used is sometimes different, the basic approach of this is the same as that applied here, i.e. the impact parameter dependent ion dissociation and ion-ion survival probabilities are accounted for and appropriately translated to transverse momentum space, with the standard LO QED γγ → l + l − amplitude applied and the ion photon flux accounted for via the usual ion EM form factors. Qualitatively we see reasonable agreement with these results, but they are not identical.The reason for this is unclear, and may lie in the precise implementation of the above theoretical ingredients.However, we note that the quoted value of A 2∆ϕ is indeed exactly zero, contrary to the discussion above, and therefore these predictions must rely on the p ee ⊥ ≪ p e ⊥ approximation, which as discussed above is not necessarily valid for the STAR data.Hence, this is a possible reason for the observed difference.
Finally, we end this section by noting that some care is needed in interpreting the comparisons made in the STAR analysis [27].In particular, in the previous versions of SuperChic, as is clearly described in [18] for the case of version 3 where heavy ion UPCs are first considered, ion dissociation had not yet been implemented.That is, only inclusive (with respect to ion dissociation) production could be generated.The STAR data are on the other hand taken with a XnXn tag, that is they are not inclusive with respect to ion dissociation.Given this, and as we have discussed in detail above, we would expect the SuperChic 3 predictions not to match the shape of the data distributions; this is precisely what is observed in Fig. 11 (left) and Fig. 12.However, in [27], comparisons are presented to SuperChic 3, and this difference between what the data correspond to and what is actually generated by the Monte Carlo is not highlighted in the comparisons2 .Indeed, the observed discrepancy between the SuperChic 3 predictions and the data is claimed in [27] to invalidate the theoretical approach described in [18].It should be clear from the discussion here that this is not the case.It is furthermore stated that the treatment of [18] neglects the fact that energy spectrum of the colliding photons depends on the nucleusnucleus impact parameter, and therefore, on the spatial distribution of the electromagnetic fields, with the implication being that this is the reason for the discrepancy.This is not correct: as discussed in detail in Section 2, this dependence (i.e. the correlation between the ion-ion impact parameter and the measured kinematic distributions) is fully accounted for.It is simply that mutual ion-ion dissociation was not yet accounted for in SuperChic 3, and hence it should not have been used for comparison to the STAR data.
The above confusion appears to motivate the discussion in [69], where the approach presented here and in previous studies is incorrectly associated with the equivalent photon approximation (as discussed in Section 2.1 this approximation is not applied in our calculation) and an apparently artificial distinction made between the approach discussed here and the result of a 'QED' calculation [33,68].This appears to be based on the incorrect assumption that the method outlined here cannot account for impact parameter constraints of the sort implied by mutual ion-dissociation conditions (see in particular the discussion in Section 3.3 of [69]).It should be clear from the results presented in this paper that this is not the case.The 'QED' results of [33,68] and those of this work are built from precisely the same underlying ingredients, and account for the correlation between the ion-ion impact parameter and the measured kinematic distributions in the same manner.Indeed, the predictions of [68] with respect to both the STAR data and the CMS measurement presented in Section 3.2 agree rather well with our results.

Summary and Outlook
Ultraperipheral heavy ion collisions (UPCs) provide a source of photon-photon collisions with which we can probe new phenomena within and potentially beyond the Standard Model.However, to leverage this initial-state to its maximum potential requires that the underlying process be modelled with as much precision as possible.In this paper we have discussed a new development of an approach which aims to provide this.We have in particular presented an update to the SuperChic MC generator to account for mutual ion dissociation in UPCs, due to additional photon exchanges between the colliding ions, whereby one or both ions can be excited into a higher energy state that subsequently decays by emitting a single or multiple neutrons.This is released in the SuperChic 4.2 MC, the code and a user manual for which can be found at http://projects.hepforge.org/superchicThis simulates a range of exclusive processes in both heavy ion and proton collisions, see [18,42] for more details.
The theoretical framework presented in this paper can provide a full account of the relevant physics effects in UPCs.In particular, while the survival factor probability of no additional inelastic ion-ion scattering due to the strong interaction, and the mutual ion dissociation probability, are both most straightforwardly formulated in the impact parameter space of the colliding ions, we have discussed how these can consistently be accounted for and translated to momen-tum space.In this way we can systematically account for the modifications these entail for the measured kinematic distributions of the centrally produced state.While the importance of such a systematic treatment in the case of the survival factor has been discussed elsewhere (see e.g.[21]), here we have focused on the case of mutual ion dissociation, which is included for the first time in Superchic.This is of particular relevance given it is possible to measure such ion dissociation processes via zero degree calorimeter (ZDC) detectors at the LHC and RHIC.The previous version 3 of Superchic is inclusive with respect to ion dissociation, and hence can be reliably used when the data are also presented inclusively.However, it is possible with ZDCs to presented data for a range of neutron tags, that is with and without ion dissociation requirements in place.Such processes can now be modelled consistently within the same overall framework.
We have compared to a range of data from ATLAS, CMS and STAR and demonstrated how these ion dissociation requirements (or vetos) in the ZDCs have a non-negligible impact on the predicted kinematic distributions, and observed rather good agreement between the data and theory for electron and muon pair production.On the other hand, some apparent discrepancies exist in certain regions of phase space.Some of these we can understand to be due to higher order QED corrections, in particular due to photon FSR from the leptons.A full account of this, and indeed other higher QED effects due to lepton-ion interactions and unitary corrections relating to additional electron pair production, is therefore well motivated as future work, and is left to that.Further mild differences may indicate the need for further tuning in the input γA → A * cross sections, and we provide an approximate method to achieve this directly in the SuperChic 4.2 code.Further data on e.g. the 1n0n and 1n1n channels, as well as the Xn0n, XnXn channels will shed more light on this.
In summary, in this work we have presented a significant update to the Superchic MC to account for mutual ion dissociation, allowing this to be used for the first time to compare against the range of data that can be (and has been) collected in high energy heavy ion collisions with ZDC tags in place.The agreement between data and theory for electron and muon pair production is observed to be in general good, providing key support for the approach presented here and motivation to use it to interpret data on less well tested final states, such as τ leptons (and the implications for the anomalous magnetic moment of the τ ) and light-by-light scattering.

Figure 1 :
Figure 1: (Left) Photon absorption cross section, σ(γA → A * ), for lead ions, as a function of the photon energy ω in the dissociating ion rest frame.(Right) Photon absorption cross section weighted by corresponding flux factor as in (34) for b ⊥ = 3RA, and √ s nn = 5.02 TeV.

Figure 3 :
Figure 3: (Left) Breakup probabilities for single and multiple neutron emission, as defined in (39) and (33), for PbPb collisions at √ snn = 5.02 TeV.The probability that no inelastic ion-ion scattering occurs, as introduced in (19), is indicated by the dashed black curve.(Right) As in the left figure, but now multiplied by the no inelastic ion-ion scattering probability, i.e. including the survival factor.

f 8 <f nn , 1 . 6 <Figure 4 :
Figure 4: Comparison of SuperChic 4.2 predictions to ATLAS data [24] on ultraperipheral muon pair production in PbPb collisions at √ snn = 5.02 TeV as a function of the dimuon invariant mass and for different dimuon rapidity regions.Results for the ratio of the Xn0n and XnXn cross sections to the inclusive UPC case (with respect ion dissociation) are shown.The muons are required to have p ⊥,µ > 4 GeV, |ηµ| < 2.4, mµµ > 10 GeV and p ⊥,µµ < 2 GeV.Data errors correspond to systematic and statistical added in quadrature.

fffFigure 5 :
Figure 5: As in Fig. 4 but now shown as a function of the dimuon rapidity, and for different dimuon invariant mass regions.

Figure 6 :
Figure 6: Comparison of SuperChic 4.2 predictions to ATLAS data [25] on ultraperipheral electron pair production in PbPb collisions at √ snn = 5.02 TeV as a function of the dimuon invariant mass and for different dimuon rapidity regions.Results for the ratio of the 0n0n, Xn0n and XnXn cross sections to the inclusive UPC case (with respect ion dissociation) are shown.The electrons are required to have p ⊥,e > 2.5 GeV, |ηe| < 2.4, mee > 5 GeV and p ⊥,µµ < 2 GeV.Data errors to systematic and statistical added in quadrature.

Figure 8 :
Figure 8: Comparison of SuperChic 4.2 predictions to ATLAS data [25] on ultraperipheral electron pair production in PbPb collisions at √ snn = 5.02 TeV in the 0n0n channel, and for a range of kinematic variables.The electrons are required to have p ⊥,e > 2.5 GeV, |ηe| < 2.4, mee > 5 GeV and p ⊥,µµ < 2 GeV.Data errors correspond to systematic and statistical added in quadrature, and are shown by the grey band in the data/theory ratios.

Figure 10 :
Figure 10: Comparison of SuperChic 4.2 predictions to CMS data [26] on ultraperipheral muon pair production in PbPb collisions at √ snn = 5.02 TeV as a function of the dimuon acoplanarity, α, for different neutron tags.

Figure 11 :
Figure 11: Comparison of SuperChic 4.2 predictions to STAR data [27] on ultraperipheral electron pair production in AuAu collisions at √ snn = 200 GeV as function of the azimuthal angular separation, ∆ϕ, definedin the text.In the left plot the prediction for the inclusive (with respect to the neutron tag) case, both with and without including the ion-ion survival factor, are shown for comparison.In the right plot the predicted distribution with respect to the alternative variable, ∆ϕ ′ , defined in the text, as well as the fit of the functional form(53) to the data are shown.The electrons are required to have p ⊥,e > 0.2 GeV, |ηe| < 1.0, 0.4 < mee < 2.6 GeV, and p ⊥,ee < 0.1 GeV.Data errors correspond to systematic and statistical added in quadrature.

Table 1 :
Predicted values of A 2,4∆ϕ (53)s argued that the A 2∆ϕ coefficient in(53)should be zero up to ∼ m 2 e /m 2 ee electron mass corrections.However, we can see from Table1that this is not the case for our prediction.The reason for this comes from considering precisely what is assumed in these analyses, namely that the photon transverse momenta are much smaller than the electron transverse momenta, i.e. p ee ⊥ ≪ p e ⊥ .However, the STAR data extend down to a minimum p e ⊥ > 200 MeV, which is indeed larger than the root mean squared dielectron transverse momentum, ⟨(p ee ⊥ ) 2 ⟩ ∼ 40 MeV, but not so large that one can necessarily neglect the Figure 12: Comparison of SuperChic 4.2 predictions to STAR data [27] on ultraperipheral electron pair production in AuAu collisions at √ snn = 200 GeV as function of the dielectron transverse momentum, p ⊥,ee .The electrons are required to have p ⊥,e > 0.2 GeV, |ηe| < 1.0, 0.4 < mee < 2.6 GeV, and p ⊥,ee < 0.1 GeV.Data points are extracted from [27]; as precise p ⊥ binning not publicly available theoretical results are presented as curves only.