Two-particle azimuthal correlations in photonuclear ultraperipheral Pb + Pb collisions at 5 . 02 TeV with ATLAS

Citation for published version (APA): Aad, G., Abbott, B., Abbott, DC., Abud, AA., Abeling, K., Abhayasinghe, D. K., Abidi, S. H., AbouZeid, O. S. A., Abraham, NL., Abramowicz, H., Dam, M., Camplani, A., Ignazzi, R., Hansen, J. B., Hansen, P. H., Hansen, J. D., Wiglesworth, G., Petersen, T. C., Xella, S., & ATLAS Collaboration (2021). Two-particle azimuthal correlations in photonuclear ultraperipheral Pb+Pb collisions at 5.02 TeV with ATLAS. Physical Review C, 104(1), [014903]. https://doi.org/10.1103/PhysRevC.104.014903


I. INTRODUCTION
In ultrarelativistic collisions of lead nuclei at the CERN Large Hadron Collider (LHC), the typical processes studied are those for which the nuclei have an impact parameter less than twice their radius (b 2R A ). Such lead-lead (Pb+Pb) collisions are understood to create a quark-gluon plasma and result in a large number of particles in the final state which participate in collective motion as a consequence of the plasma evolution [1][2][3]. In addition to the particles produced in Pb+Pb collisions, those produced in high-energy proton-proton (pp) and proton-lead (p + Pb) collisions also exhibit a collective behavior which manifests as an event-wide azimuthal variation persisting broadly in pseudorapidity, initially observed as a "ridge" [4][5][6][7]. This behavior is characterized in terms of nonzero single-particle azimuthal anisotropies, given by nth-order Fourier coefficients, and here referred to as flow coefficients v n . Nonzero v n values have also been observed in significantly lower-energy  d + Au collisions at the BNL Relativistic Heavy Ion Collider (RHIC) [8]. A natural question is whether such signatures persist in even smaller collision systems [9] and, if so, how this may influence the interpretation of these signatures in pp or p + Pb collisions.
In the prevailing paradigm, these observed anisotropies arise from the creation of a miniature region of quark-gluon plasma [10,11], in which a hydrodynamiclike expansion of the system converts spatial nonuniformities in the initial state of the system into momentum-space anisotropies of the final-state particles. However, momentum correlations already present in the initial state of the collision may also persist into the final state [12]. The relative importance of these two explanations can be tested in collision systems where one or both of the "beams" has a significantly simpler initial state. Recently, studies were performed in archived e + e − collision data at √ s = 91 GeV from the ALEPH detector [13] and in archived ep collision data at √ s = 316 GeV from the ZEUS detector [14] with a deep-inelastic-scattering selection Q 2 > 1 GeV 2 . A ridge signature was not observed, and experimental upper limits were set on the possible magnitude of v n coefficients.
In addition to the hadronic Pb+Pb interactions described above, the strong electromagnetic (EM) fields of the fully ionized nuclei can induce interactions even when the nuclei have significantly larger impact parameters such that no hadronic interaction occurs (b 2R A ). In the equivalent photon approximation [15][16][17], these strong EM fields correspond to a flux of quasireal, high-energy photons. These photons can be emitted coherently from the entire nucleus, producing a flux enhanced by a factor of Z 2 (Z = 82 for Pb) for photons up to 80 GeV at the LHC. These ultraperipheral collisions (UPCs) [18,19] have appreciable rates and include photonphoton (γ γ ) and photonuclear (γ + A) interactions.
At the LHC, previous measurements of ultraperipheral processes in Pb+Pb collisions include light-by-light scattering (γ γ → γ γ ) [20][21][22][23], exclusive dilepton production (γ γ → e + e − and γ γ → μ + μ − ) [24][25][26], and the photonuclear production of various meson states (γ + A → h + X ) [26][27][28]. In the photonuclear case, the photon may act as a pointlike particle interacting with a parton in the nucleus (the "direct" case). However, the vector-meson dominance picture [18,29] suggests that the photon often fluctuates into a vectormeson state such as a ρ or ω (the "resolved" case). In this case, the interaction proceeds as a meson-nucleus collision at an energy lower than that of the associated nucleon-nucleon collision. Figure 1 illustrates the direct and resolved photonuclear interactions. The photon-nucleon collision energy and the boost of the center of mass relative to the nucleus-nucleus rest frame depend on the photon energy and thus vary event to event. For photons with energies at the upper boundary of the coherence region, E = 80 GeV, the resulting photon-nucleon center-of-mass energy is approximately 900 GeV. Thus photonuclear collisions may be used to probe the dynamics of a system with a novel energy and geometry compared to pp or p + A collisions at the LHC, and to e + e − or ep collisions at the CERN Large Electron-Positron Collider (LEP) and the Hadron-Electron Ring Accelerator (HERA) at DESY. Since photonuclear events are the photoproduction limit of deep inelastic scattering on nuclei, these measurements may also shed light on possible collective signatures at the future Electron Ion Collider [30,31].
This paper presents a measurement of azimuthal anisotropies obtained via two-particle correlations in photonuclear collisions, where such analyses have not previously been undertaken. The data were recorded using a trigger designed to select minimum-bias and high-multiplicity photonuclear events in 1.7 nb −1 of Pb+Pb collisions at 5.02 TeV per nucleon pair delivered by the LHC in 2018. Photonuclear event candidates are selected, and are distinguished from peripheral hadronic Pb+Pb events and other background events, using the topology of the distribution of particles in the event as measured in the zero-degree, forward, and barrel calorimeters, as well as the tracking systems through the reconstruction of pseudorapidity gaps [32,33]. Properties of the selected events are compared with the expectations from Monte Carlo (MC) simulations of photonuclear processes. Two-particle correlations as a function of relative separation in azimuth ( φ) and pseudorapidity 1 ( η) are constructed for different selections of event charged-particle multiplicity and charged-particle kinematics. The nonflow contributions to the two-particle correlations (for example, jet correlations which give rise to a φ correlation structure) are suppressed by studying the correlations at large η, and the residual nonflow contribution is subtracted via the template-fitting method used extensively in prior ATLAS measurements [4,5,34]. In the template method, the correlation in high-multiplicity events is described as a combination of the correlation in lower-multiplicity events plus a component modulated by cos(2 φ) (and higher order) Fourier terms. The template method makes particular assumptions about how the nonflow component evolves with multiplicity. Although there are differences between the system explored in this measurement and those in previous two-particle correlation measurements, the sensitivity to the assumptions of the template method can be tested within the standard approach. A test of the template method in simulated photonuclear events which do not include flow or initial-state correlation mechanisms is performed in Sect. VII. The resulting magnitudes of the two-particle correlations are interpreted as arising from the product of global v 2 and v 3 values for individual particles, and are reported as a function of the reconstructed charged-particle multiplicity (N rec ch ) and transverse momentum (p T ). The results are compared with other small collision systems at the LHC, and theoretical expectations from initialand final-state physics mechanisms are discussed. 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector and the z-axis along the beam pipe. The x axis points from the IP to the center of the LHC ring, and the y axis points upwards. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2).

II. ATLAS DETECTOR AND DATA SAMPLE
The ATLAS detector [35] covers nearly the entire solid angle around the collision point. The detector systems used for the measurements presented in this paper include the inner tracking detector, the electromagnetic (EM) and hadronic calorimeters, the zero-degree calorimeter (ZDC), and the trigger and data acquisition systems. The detector halves at positive and negative z values are referred to as the A and C sides, respectively.
The inner-detector system is immersed in a 2 T axial magnetic field and provides charged-particle tracking in the range |η| < 2.5. The high-granularity silicon pixel detector covers the vertex region and typically provides four measurements per track. The innermost layer, the insertable B layer [36][37][38], has been operating as a part of the silicon pixel detector since 2015.
The pixel detector is followed by the silicon microstrip tracker (SCT) which usually provides four measurement points per track. These silicon detectors are complemented by the transition radiation tracker, which enables radially extended track reconstruction up to |η| = 2.0.
The calorimeter system covers the pseudorapidity range |η| < 4.9. Within the region |η| < 3.2, electromagnetic calorimetry is provided by barrel and endcap high-granularity lead/liquid-argon (LAr) calorimeters, with an additional thin LAr presampler covering |η| < 1.8. Hadronic calorimetry is provided by a steel/scintillator-tile calorimeter, segmented into three barrel structures within |η| < 1.7, and two copper/LAr hadronic endcap calorimeters. The solid angle coverage is extended with the forward calorimeter (FCal), composed of copper/LAr and tungsten/LAr modules optimized for EM and hadronic measurements, respectively. The FCal detectors cover the regions 3.2 < |η| < 4.9.
The minimum-bias trigger scintillator (MBTS) detects charged particles over 2.07 < |η| < 3.86 using two hodoscopes of 12 counters positioned at z = ±3.6 m. The ZDC consists of EM and hadronic sections and plays a key role in identifying UPC events in heavy-ion collisions by primarily detecting neutrons resulting from the breakup of one or both nuclei. The ZDC modules are located at z = ±140 m from the IP. They measure neutral particles at pseudorapidities |η| 8.3 and consist of layers of alternating quartz rods and tungsten plates.
A two-level trigger system [39] is used to select events. The first-level trigger (L1) is implemented in hardware and uses a subset of the detector information to restrict the accepted rate to at most 100 kHz. This is followed by a software-based highlevel trigger (HLT) stage that reduces the accepted event rate to 1-4 kHz depending on the data-taking conditions during 2018 Pb+Pb operations.
The measurements presented in this paper were performed using the √ s NN = 5.02 TeV Pb+Pb dataset collected with a variety of triggers in 2018, with a total integrated luminosity of 1.7 nb −1 . Photonuclear candidate events were first selected by the trigger by requiring one ZDC side (referred to as the Pb-going side) to have a minimum amount of energy at L1, E > 1 TeV, consistent with the presence of one or more neutrons. The other side (referred to as the photon-going side) was required to have an energy below a maximum-energy cutoff, E < 1 TeV, consistent with no neutrons. The cut value of E = 1 TeV is several multiples of the energy resolution away from the single-neutron peak at E = 2.5 TeV [24]. Thus, the selected topology is referred to as "0nXn" in the figures in this paper. Events were also required to satisfy an upper bound of 200 GeV on the total transverse energy deposited in the calorimeter at L1, for further rejection of hadronic Pb+Pb events.
After these requirements at L1, events had to pass either a minimum-bias (MB) trigger or one of several highmultiplicity triggers (HMT) with further requirements at the HLT stage. The MB trigger was defined by requiring at least one reconstructed online track with p T > 0.4 GeV. The HMTs were defined by requiring a larger number of tracks originating from the same vertex. The primary vertex, defined as the one with the largest scalar sum of p 2 T values of associated tracks, was used. An HMT with a threshold of 15 tracks was defined by also requiring a hit in the MBTS at L1. Two HMTs with thresholds of 25 and 35 tracks were defined by requiring a minimum amount of energy in the calorimeter at L1. Finally, all HMTs had an additional HLT requirement of less than 5 GeV of transverse energy in the photon-going FCal, which rejected a large fraction of the 0nXn peripheral Pb+Pb background events. The 35-track HMT sampled the full luminosity of 1.7 nb −1 during data taking. The 25-track HMT, 15-track HMT, and the MB trigger were configured with increasingly larger prescale factors, sampling 1.6 nb −1 , 0.13 nb −1 , and 1.0 μ b −1 of data respectively.

III. MONTE CARLO SIMULATION
Monte Carlo (MC) simulations of the relevant physics processes were used to understand the performance of the detector, and provide distributions to be compared with the data. For all the generators and configurations listed below, the simulated events were passed through a GEANT4 simulation [40,41] of the detector and reconstructed under the same conditions as the data.
A sample of one million peripheral Pb+Pb events was generated using HIJING v1.383 [42] with impact parameters in the range 10 < b < 20 fm. These events were used to determine the charged-particle tracking efficiencies for the measurement, and to model the event properties of low-multiplicity, hadronic Pb+Pb events.
Several generators were used to simulate photonuclear interactions. For each of the cases below, events were generated with different minimum requirements on N rec ch to provide good statistical coverage over the N rec ch range accessed in data. First, the distribution of photon flux for 208 Pb beams at the LHC was calculated using STARLIGHT [43]. The flux distribution was passed to a multipurpose generator based on the dual parton model (DPM) and referred to as DPMJET-III [44,45], which simulated direct and resolved photon-lead (γ + Pb) interactions at the generator level, followed by a full GEANT4 simulation of the ATLAS detector. Thirteen million γ + Pb events were simulated in this way. Additionally, the flux from STARLIGHT was used to simulate two million γ + p events using DPMJET-III, where the energy of the proton was set  equal to the energy per nucleon in the Pb beams. Finally, for an alternative description of these processes, PYTHIA 8.240 [46] configured with the NNPDF23LO parton distribution functions [47] and A14 set of tuned parameters [48] was used to generate twelve million γ + p events. The photon flux in PYTHIA was reweighted at the event level to match that calculated by STARLIGHT, and the simulation was configured to include both the direct and resolved photon interactions.

A. Reconstruction and event selection
Given the low particle multiplicities of UPC events, the charged-particle track and calorimeter energy-cluster reconstruction procedures follow those optimized for pp data taking [49,50]. Reconstructed charged-particle tracks are used in the analysis if they satisfy quality criteria as outlined in Ref. [51]. These include criteria for the minimum number of hits in the pixel and SCT detectors. Tracks are further required to have p T > 0.4 GeV, |η| < 2.5, and distances of closest approach to the primary reconstructed vertex in the longitudinal and transverse directions of less than 1.5 mm each. Clusters in the range |η| < 4.9 are constructed from topologically connected groups of calorimeter cells [52]. They are required to have p T > 0.2 GeV and to meet the significance criteria for the measured energy as outlined in Ref.
[32] to suppress the contribution from electronic noise fluctuations.
Events are required to have a reconstructed primary vertex within |z| < 90 mm. Events compatible with multiple in-time interactions are tagged by the presence of multiple reconstructed vertices. These comprise less than 1% of the sample and are rejected. Events must pass a version of the ZDC-based energy selection described in Sec. II which uses the same threshold but an improved calibration for the ZDC energies determined after data taking. Events are further characterized by their charged-particle multiplicity, N rec ch , which is defined as in previous ATLAS measurements of correlations in small systems [4,5,7] as the total number of reconstructed tracks, without efficiency correction, with p T > 0.4 GeV and |η| < 2.5. For the correlation analysis described in Sec. V, events with a given N rec ch were included in the analysis if they were selected by any of the HMTs which was more than 80% efficient for events with this N rec ch value. For the study of event properties in Sec. IV B, each N rec ch range was only populated by the highest-sampled-luminosity HMT which was more than 99% efficient over the entire N rec ch range. In either case, events with N rec ch < 15 were populated only by the MB trigger.
Reconstructed pseudorapidity gap quantities, constructed using the tracks and clusters in each event, are used to distinguish between different physics processes such as photonuclear collisions, low-activity (peripheral) hadronic Pb+Pb collisions, and γ γ → X processes. The requirement of a rapidity gap above a minimum value in the photon-going direction can efficiently remove peripheral Pb+Pb events. Rather than the traditional pseudorapidity gap quantity [32], which determines the pseudorapidity difference between the edge of the detector and closest particle, an alternative "sum-of-gaps" definition is used, which adds together contiguous gaps separated by particle production concentrated in a narrow pseudorapidity regions. This alternative definition is used to retain a large selection efficiency for resolved photon events which may break up a large gap with a hadronic fragment localized in pseudorapidity (see Fig. 1). The quantities γ η and A η correspond to the sum-of-gaps values calculated in the photon-going and nucleus-going halves of the detector, respectively. They are calculated by first sorting the tracks and clusters in η. The differences in η between adjacent particles, η, are included in the sum if they are larger than 0.5. The value of 0.5 was observed in simulation to retain good efficiency for resolved photon events. The position at η = 0 is treated as if it were a particle, which ends the sum-of-gaps calculation. Thus, γ η and A η range from 0 to 4.9.
In addition to the sum-of-gaps quantities defined above, the traditional edge gap quantity, constructed using clusters, η cluster edge , is defined as the pseudorapidity difference between the edge of the detector on the photon-going side, at η = +4.9, and the first reconstructed cluster. In the mixed-event construction discussed in Sec. V, events are classified by their η cluster edge as a way to characterize their overall topology without using charged-particle tracks.

B. Event properties
Figure 2 (left) shows the distribution of γ η and N rec ch values for events recorded with the MB trigger. The γ η can have a value of zero or fall into the range 0.5 γ η 4.9, with the lower value determined by the minimum gap size included in the sum. Events resulting from hadronic Pb+Pb collisions are at small γ η and have a broad multiplicity distribution, whereas events from photon-induced processes are at large γ η and have an N rec ch distribution that falls steeply. Figure 2 (right) shows the distribution of γ η and A η values. Most hadronic Pb+Pb collisions have a broad pseudorapidity distribution of particles and hence have both small γ η and small A η. In contrast, the photonuclear events have large γ η but smaller values of A η. There is no significant yield of events with a large sum of gaps on both sides, which could signal the presence of photonphoton hadronic process backgrounds. The signal events for photonuclear collisions are defined by γ η > 2.5. No event selection is made on A η. Figure 3 summarizes the multiplicity distribution observed in data after the event selection described above. The N rec ch distribution for events passing the photonuclear event selection is shown in the left panel of Fig. 3. The sawtooth pattern arises from the inclusion of the HMTs at different N rec ch thresholds as described above. In each N rec ch range, the selected photonuclear events have a steeply falling multiplicity distribution. The N rec ch distribution for photonuclear events, fully corrected for the different luminosities sampled by the triggers and the N rec chdependent trigger efficiency, is shown using black circles in the right panel of Fig. 3. It is compared with the distribution from events with γ η < 1, which mostly selects Pb+Pb hadronic events, which are backgrounds in this analysis. The selected photonuclear events have a significantly more steeply falling multiplicity distribution.  shows the total N rec ch distribution, corrected for the trigger efficiency and the different integrated luminosities sampled by the HMTs. This is compared with three distributions from MC generators, with the same gap-based selection requirements as the data. The ZDC selection is applied to data but not to the MC samples as the generators do not model the nuclear-breakup processes relevant for forward neutron spectators. Including a ZDC requirement in MC events may thus impact the distributions if particle production in the detector is correlated with the nuclear fragmentation in the forward region. Despite the limitations in the modeling of forward neutrons, a generator-level check requiring neutrons in the ZDC acceptance was performed and found not to impact the level of agreement between data and DPMJET-III. The distributions in PYTHIA and DPMJET-III γ + p are normalized to have the same integral as the data over the full N rec ch range. The models show good agreement with the data at low N rec ch , but systematically predict too low a relative yield at higher N rec ch . DPMJET-III γ + Pb does not describe the full distribution, but has been normalized to have the same integral as the data in N rec ch > 35 to highlight its good agreement in this region over many orders of magnitude. With this normalization, DPMJET-III γ + Pb systematically predicts too low a relative yield for N rec ch < 15. This comparison either suggests the presence of other, nonphotonuclear processes in data at such low N rec ch < 10 values, or points to the need for improved modeling of this region in the simulation.
The right panel of Fig. 4 shows the reconstructed γ η distributions in data and simulation for selected events with N rec ch > 10, without the γ η > 2.5 requirement. Structures in the distributions correspond to transitions between detector subsystems and the change in the detector response as a function of η. At large γ η values 2.5, the shape FIG. 3. Left: Distribution of event charged-particle multiplicity, N rec ch , for the photonuclear event selection from multiple triggers, without prescale correction. Right: Distribution of event charged-particle multiplicity, N rec ch , for the photonuclear event selection (black) with corrections for trigger efficiency and prescale factors, and an alternative selection intended to select hadronic Pb+Pb events (red, see text). of the distribution in data is qualitatively similar to that in DPMJET-III γ + Pb and Pythia γ + p simulation. However, the distributions in the simulated photonuclear events decrease at smaller γ η values, while the distribution in data rises. At low γ η, the shape in data is qualitatively similar to that in peripheral HIJING Pb+Pb events. This comparison suggests that the trigger-selected events contain a mixture of peripheral Pb+Pb events and genuine photonuclear events, with the latter dominant at γ η > 2.5. The possible impact of residual peripheral Pb+Pb events in the set of selected events is discussed in Sec. VI. Figure 5 compares the charged-particle pseudorapidity distribution, dN ch /dη, in data and simulation. The left panel shows the dN ch /dη in data, for charged particles with 0.4 < p T < 5 GeV, for multiple N rec ch selections in photonuclear events. The distributions are corrected for tracking efficiency on a per-track basis, which ranges from 0.7-0.9 depending on track η and p T . To compare the relative shapes between N rec ch selections, the distributions are each normalized to have an integral of one. In all cases, the pseudorapidity distributions are strongly asymmetric, peaking at η = −2.5 (the nucleus-going direction) and then monotonically decreasing until η = +2.5 FIG. 5. Left: Charged-particle pseudorapidity distribution, dN ch /dη, in selected N rec ch ranges. The distributions are normalized to the same integral and are shown in arbitrary units. Here, positive and negative η denote the photon-going and nucleus-going directions, respectively. Right: dN ch /dη distribution in data for N rec ch > 10 (black points), normalized per event, and compared with that in DPMJET-III γ + Pb (dotdashed green histogram), PYTHIA γ + p (dashed blue histogram), peripheral HIJING Pb+Pb (solid magenta histogram), and DPMJET-III γ + p (dotted red histogram) with the same reconstruction-level selection as the data. All distributions have been normalized to have the same value as DPMJET-III γ + Pb at η = 0.
(the photon-going direction). This overall shape is even more asymmetric than that observed in, for example, p + Pb collisions at the LHC [53]. The dN ch /dη shape remains similar over a wide range of multiplicities (N rec ch > 10). The selection with N rec ch < 5 results in an even more strongly asymmetric dN ch /dη shape, suggesting that it may include photonuclear events at significantly lower energies or other processes. Additionally, as discussed above, the region N rec ch < 10 is not well described by simulation, and may contain a significant contribution from background events. On the basis of these observations, the two-particle correlation analysis is performed using events with N rec ch 10. Additionally, an upper bound of N rec ch = 60 is used to remove the region with very few events.
The right panel in Fig. 5 compares the pseudorapidity distribution of particles with 0.4 < p T < 5 GeV in data for N ch rec 10 with the four sets of simulated events, which have the same reconstruction-level selections as the data. While all three generators considered here systematically fail to predict a steep enough slope for the dN ch /dη distribution, the distribution for DPMJET-III γ + Pb events is the most similar to that observed in data.
In general, no single generator is able to quantitatively describe in detail all of the features observed in data. The results here may be used to benchmark future developments in the modeling of photonuclear collisions. The two-particle correlation analysis described in Sec. V, due its data-driven nature, does not rely on a detailed description of photonuclear events in simulation and therefore its results are not sensitive to this mismodeling.
Since the DPMJET-III γ + Pb simulation and the data agree qualitatively in the N rec ch and γ η distributions, the simulated events provide an estimate of what ranges of photon energy in the laboratory frame, E γ , and photon-nucleon center-of-mass energy, W γ N , are selected for analysis by choosing various N rec ch or γ η values. These ranges are shown in Fig. 6 as a function of N rec ch . In DPMJET-III, photonuclear events with a larger N rec ch have a larger average E γ and W γ N . At N rec ch = 10, the mean values are E γ = 26.8 GeV and W γ N = 519 GeV, rising to E γ = 123 GeV and W γ N = 1.11 TeV at N rec ch = 60. However, even a narrow range of N rec ch values selects events with a broad distribution of E γ and W γ N . At lower values of N rec ch < 10, the mean photon energy decreases rapidly and the range of photon energies increases, supporting the data-driven minimum value of N rec ch determined above.

V. TWO-PARTICLE CORRELATIONS
This section describes the two-particle correlation and template fit procedures used to extract the azimuthal anisotropies, v n . These procedures have been used in previous analyses of two-particle correlations in small systems [4,5,34,54,55], and are summarized here.
Correlations between pairs of charged particles are reported as a function of the pair's relative separation in pseudorapidity, η = η a − η b , and azimuth, φ = φ a − φ b . The labels a and b denote kinematic selections on the first and second particle, which are generally different. If two particles FIG. 6. Correlation between photon energy E γ and event multiplicity N rec ch in simulated DPMJET-III γ + Pb events. The right-hand axis shows the center-of-mass energy W for the photon-nucleon collision system. The markers and vertical errors bars show the mean and standard deviation, respectively, of E γ as a function of N rec ch .
each meet selections a and b, the pair is tabulated twice, once for each possible assignment. By convention, the distributions in η and φ, Y ( φ, η), are averaged over events by normalizing them by the total yield of particles in selection a, N a , To correct for the charged-particle reconstruction efficiency, entries in the distribution are weighted by the product of the inverse efficiencies evaluated for the kinematics of the two particles. These efficiencies are derived from the HIJING simulation described in Sec. III. In this analysis, all chargedparticle tracks within |η| < 2.5 are considered, except the tracks which lie in the region η > 4.9 − γ η (with positive η defined as the photon-going direction), if the sum of gaps extends into the tracker. There are two reasons to exclude these tracks. First, in simulations of photonuclear events, such particles are typically fragments of resolved photons, and not the bulk particles produced in the γ + Pb interaction which are being measured for possible collective signatures. Second, since these tracks are typically near the edge of the inner detector's acceptance, the | η| restriction described below would result in a small number of tracks being paired many times, thus making a disproportionately large contribution to the correlation function.
The one-dimensional φ normalized pair yields, Y ( φ), are constructed by integrating over a restricted | η| range, typically 2.0 to 5.0, designed to suppress the contributions to the correlation from nonflow effects (such as jets, particle decays, and resonances), as FIG. 7. Two-dimensional normalized particle pair distributions in photonuclear events, corrected for acceptance effects with the mixedevent distribution, and presented as a function of η and φ. The peak at ( φ, η) = (0, 0) is truncated to better show the structure of the correlation function. Each panel represents a different charged-particle multiplicity N rec ch range for the selection 0.
In principle, due to finite coverage and detection efficiency introduced by a real detector, Y may contain acceptance effects that do not reflect physical particle pair correlations. To account for this, a mixed-event technique may be applied [56], in which particles a and b are taken from different events that have similar overall characteristics. The two events are required to have a ZDC signature on the same side (i.e., same directions of the photon and nucleus), and a similar charged-particle multiplicity, primary vertex position, and cluster edge gap η cluster edge . The mixed-event distribution d 2 N mixed /(d φ d η) contains the effects of the pair acceptance but not the physical correlations.
The acceptance effects can be corrected for by dividing the same-event pair distribution by the mixed-event ones. The two-dimensional (2D) correlation function corrected for acceptance effects is where N a is the total yield of particles with selection a, N mixed is the yield in mixed events, and N mixed pair is the total integral of the mixed-event distribution d 2 N mixed /(d φ d η). Without the mixing correction, the 2D Y ( φ, η) distributions have an artificial triangular shape in η, reflecting the convolution of the single-particle acceptances |η a , η b | < 2.5. This effect is removed by the mixing correction, and the η dependence of structures in the 2D C( φ, η) distributions become more evident.
Correlation functions are constructed for various selections on N rec ch and p a T , with the p T of particle b always in the range 0.4 < p b T < 2 GeV. All events are given an equal statistical weight, and no corrections for the trigger efficiency or differences in sampled luminosity (for N rec ch ranges that span contributions from multiple triggers) are performed. Although the mixed-event correction is necessary to properly construct the two-dimensional correlation functions presented below, it was found to ultimately have only a minor effect on the projected, one-dimensional Y ( φ) distributions, producing compatible results but with increased statistical uncertainty. Thus the mixing correction is not applied for the nominal v 2 and v 3 results. Instead, it is used in the determination of systematic uncertainties associated with the pair acceptance of the detectors as described in Sec. VI.
Examples of two-dimensional C( φ, η) correlation functions are presented in Figs. 7 and 8. In Fig. 7, particles a and b are required to have 0.4 < p a T < 2 GeV, and two example N rec ch selections are shown. In Fig. 8, correlation functions are presented for 20 < N rec ch < 60, with two example p T selections for particle a shown. The two-dimensional correlation functions have features which are broadly similar to those observed in pp collisions. There is a localized "nearside" peak at ( φ, η) ≈ (0, 0) from correlations between jet fragments, and an extended "away-side" ridge at φ ≈ π which extends over a large η range from correlations between fragments of azimuthally opposite jets. A quantitative analysis and nonflow subtraction are necessary to discern if there are additional ridge structures, as detailed below.

A. Nonflow subtraction
To remove the contribution to the correlation function from nonflow effects that persist as long-range contributions on the away side ( φ ∼ π ), a template fit procedure is employed. To perform the fit procedure, samples corresponding to events with low multiplicity (LM) and high multiplicity (HM) are selected. The shape of the nonflow contribution is assumed in the template procedure to be the same in the LM and HM samples. In this analysis, the LM events are chosen to have 15 N rec ch 20. The Y ( φ) in HM events is parametrized as the sum of an azimuthally modulated pedestal (which expresses the azimuthal anisotropy) and a nonflow component, as follows: Above, Y LM ( φ) is the correlation function in LM events, which is parameterized as a truncated Fourier series up to the fourth order (with coefficients c n for n = 0, . . . , 4). The values of F and G, the three v n,n values, and the five parameters FIG. 8. Two-dimensional normalized particle pair distributions in photonuclear events, corrected for acceptance effects with the mixedevent distribution, and presented as a function of η and φ. The peak at ( φ, η) = (0, 0) is truncated to better show the structure of the correlation function. Each panel represents a different p a T range for the selection 20 < N rec ch < 60 and 0.4 < p b T < 2.0 GeV. c n , describing the LM reference, are free parameters in the fit, but F and G are constrained such that the integrals of both sides of Eq. (1) are the same. Modulation terms up to fourth order (v 2,2 , v 3,3 , and v 4,4 ) are considered in the fit, in order to best describe the HM data. By fitting Y LM ( φ) and Y HM ( φ) simultaneously, the extracted uncertainty in F , G, and v n,n correctly accounts for the statistical uncertainty of both the LM and HM samples. An example of the simultaneous fit of the HM selection and LM reference is shown in Fig. 9. Examples of the template fit in additional HM and p a T selections are shown in Fig. 10. In the bottom panels of Figs. 9 and 10, the p values are defined following the procedure described below.
The template fit is performed by minimizing the standard χ 2 between the data points and the functional form. However, the data points within the correlation functions contain nontrivial point-to-point correlations, since a single particle b may be used in combination with multiple particles of type a. The minimum of the χ 2 statistic, when calculated in the traditional way, is found at the appropriate values of the fit parameters. However, the p value and the uncertainty in the parameter values, if also determined in the standard way, would be inaccurate. In order to properly account for these correlations and determine the parameter value uncertainties, a bootstrapping procedure was applied. Pseudoexperiments were generated by giving a random Poisson weight (with a mean of one) to  individual events. These pseudoexperiments are independent sets of events drawn from the observed set of events with the same statistics and correlations as the real data. The parameters values extracted in the pseudoexperiments agreed with those in data and followed Gaussian distributions. The standard deviations of the parameters values extracted from the pseudoexperiments were used to set the statistical uncertainties in the final results, which are 10-30% larger than if they were estimated in the naive way without taking the point-to-point correlations into account. The p values from the bootstrapping experiment are presented in the template fit plots that follow.
The full set of v 2,2 and v 3,3 values are shown in Fig. 11. The filled points show the nonflow subtracted results according to the template fit procedure in Eq. (1). The open points show the v n,n values that would be obtained if instead a direct Fourier decomposition of Y HM were performed, without any attempt to account for the non-flow contribution. For most FIG. 11. Summary of v 2,2 and v 3,3 results as obtained from a direct Fourier fit (open markers) and those obtained after the template-fitting subtraction (filled markers). Results are shown differential in N rec ch for fixed p a T (left) and differential in p a T for fixed N rec ch (right).
of the selections considered here, the nonflow subtraction has a significant effect on the extracted v 2,2 and v 3,3 values. The resulting v n,n values are positive in all selections, with one exception: in the p a T -dependent results with a single HM selection, the v 2,2 value for 3 < p T < 5 GeV is negative. Additionally, the v 2,2 value for 2 < p T < 3 GeV is significantly lower than that for 1.2 < p T < 2 GeV. In these selections, the v 3,3 values also rise significantly. The template fits to these selections are shown in Fig. 12, and are discussed further below.

B. Factorization test
In the flow paradigm, a two-particle azimuthal modulation characterized by a v n,n value arises from the product of nonzero azimuthal anisotropies, v n , for each particle. These are related via v n,n (p a T , if a and b are selected from the identical particle p T range. Thus, a single-particle flow coefficient v n (p a T ) may be determined from twoparticle v n,n values through v n (p a T ) = v n,n (p a T , p b T )/v n (p b T ) = for a given selection on reference particle b. To test whether the v n,n values in data are compatible with this picture, a factorization test can be performed in which v n values for particle a are compared for different particle b selections. The results of this test for the v 2 values as a function of N rec ch are shown in Fig. 13. The test demonstrates that, while the v 2,2 values for different p b T selections may be different, the v 2 values obtained for particle a as a function of N rec ch are independent of the selection for particle b. This is a necessary condition of the paradigm that the observed v 2,2 values arise from the product of event-wide, single-particle v 2 azimuthal modulations. Figure 14 shows the results of a factorization test as a function of particle a p T . Although this test is more statistically limited than that for the N rec ch -dependent test above, the results for p a T < 2 GeV are compatible with the factorization assumption. Since nonflow effects are not necessarily expected to factorize, this validates an important assumption of the template procedure. For the selection with p a T > 3 GeV, the v 2,2 values obtained in Fig. 14 are significantly negative, and are accompanied by large v 3,3 values. A negative v 2,2 value violates the expected factorization for particles with independent overall modulation, and thus cannot be interpreted as arising from hydrodynamic flow.
As can be seen from the growing Fourier v 2,2 values in Fig. 11 and the large away-side peaks in Fig. 12, the nonflow contribution in these high-p a T selections is significantly larger than the magnitude of any possible modulation. Thus, the extraction of any genuine cos(n φ) modulation is highly sensitive to assumptions in the template method about the scaling of the nonflow contribution between multiplicity classes. If violated, this could manifest as a reduced v 2,2 and increased v 3,3 . This assumption is difficult to test given the limited FIG. 14. Summary of the factorization test for v 2 values as a function of p a T . Comparison of v 2,2 (left) and derived v 2 values (right) using the same particle a selection, but with different selections on particle b. Only statistical uncertainties are shown. The points are displaced horizontally for visibility.

C. Physics backgrounds
The measured v n values are potentially biased by the presence of background events, primarily from peripheral hadronic Pb+Pb events, in the selected event sample. The purity of photonuclear events is estimated by fitting templates derived from simulated photonuclear and peripheral Pb+Pb events to distributions measured in data, and determining the fraction of selected events which are compatible with photonuclear interactions.
The purity was first estimated by fitting the γ η distribution in data to the combination of distributions from DPMJET-III γ + Pb events and peripheral HIJING Pb+Pb events over the range 0.5 < γ η < 4.9, with the relative contributions from each as free parameters. An example of these two contributions and the distribution in data is shown in the right panel of Fig. 4. This fit was performed separately in each of the N rec ch selections used in the measurement. The purity for the γ η > 2.5 selection is estimated to be greater than 98% in all N rec ch selections used in the analysis. Thus, for the nominal results, no correction for a possible impurity in the selected events is performed.

VI. SYSTEMATIC UNCERTAINTIES
The primary sources of systematic uncertainty are grouped by category and discussed below. The systematic uncertainties are evaluated by varying the procedure applied to the data, and are not directly derived from the simulation samples described in Sec. IV B. Thus, the uncertainties do not depend on an accurate modeling of the photonuclear processes in simulation. The uncertainties from each source are ultimately combined in quadrature to obtain the total systematic uncertainty. For some individual sources below where only one variation of the analysis procedure was performed, the resulting uncer-tainty was symmetrized by adding an uncertainty of equal magnitude in the other direction. However, in general the final systematic uncertainties are not symmetric. The typical magnitudes of the systematic uncertainties are summarized in Table I.
Pair-pseudorapidity difference. The two-particle correlation function is measured for pairs with 2.0 < | η| < 5.0 as a way to strongly reduce the nonflow contribution, which varies much more strongly with η than the genuine long-range azimuthal modulation. To test the sensitivity of the analysis to a potential residual contribution, the lower bound on η is varied to 1.8 and to 2.2, and the impact of the resulting variations on the data is included as an uncertainty.
Low-multiplicity reference. Photonuclear events with 15 N rec ch 20 are used as the LM reference in the template fit to the other N rec ch selections. In the template method, the shape of the nonflow contribution is assumed to be the same in the LM and HM selections. To test the possible variation of the nonflow contribution with multiplicity, the analysis may be repeated using an alternative reference for LM events. Additionally, a different multiplicity range may select a different sample of photonuclear events (i.e., distribution of photon energies, or mix of resolved and direct processes [18]), to which the results may be sensitive. An alternative selection, 10 N rec ch < 15, composed from the MB trigger data, is considered for the LM reference. The analysis is repeated and the difference in the results is included in the uncertainties.
φ binning. The sensitivity to the particular choice of φ binning was evaluated by repeating the fitting procedure after halving the number of bins. The resulting difference is included in the total uncertainty.
Acceptance effects in correlation function. For the nominal analysis, the φ-dependent correlation functions are formed by projecting the 2D correlation functions in the region 2 < | η| < 5. To test the possible impact of acceptance effects that may result in a distortion of the φ distribution, the correlation functions are corrected by the projection of the mixed-event 2D correlation function onto φ, as described in Sec. V. The differences in the results are assigned as a systematic uncertainty.
Exclusion of tracks from γ η. The default analysis excludes reconstructed tracks that are within the sum-of-gaps range from the two-particle correlations, i.e., those with η > 4.9 − γ η. As a systematic variation, these tracks are included in constructing the correlation functions. The resulting differences are included as a systematic uncertainty.
Photonuclear event purity. The estimate of the purity of photonuclear events in data relies on the description of γ + Pb event properties in simulation. To test the sensitivity to this description, the fit described in Sec. V C was repeated using PYTHIA γ + p events as the model for photonuclear events and, alternatively, by fitting to the N rec ch distribution, separately in different slices of γ η. In the latter case, the purity decreases to approximately 90% in events with N rec ch = 15−20 and to 80% in those with N rec ch = 10−15, which are the ranges used to define the nominal and alternative LM selections for the two-particle correlation analysis.
In the default analysis, no correction is made for a possible presence of backgrounds, primarily peripheral Pb+Pb events, in the selected sample of high-multiplicity γ + Pb events. As a systematic variation, the φ distributions in data are corrected by subtracting the contribution from these background events, scaled to match the purity in each N rec ch range. Since the potential impact of the purity is much larger for events with low N rec ch than those with high N rec ch , the effect on the analysis is evaluated by accounting for the possible presence of backgrounds only in the LM reference. Additionally, since the purity was observed to be lower in the alternative LM reference events, N rec ch = 10−15, than in the nominal LM selection, the difference in the results using the alternative LM reference, with and without the purity correction, was taken as the uncertainty for this source. The model φ distributions for the background events are taken from pp collision data at √ s = 5 TeV [57], with the same reconstruction, definition of N rec ch , and gap-based event selection as this analysis. The resulting differences are included as a systematic uncertainty.
Other uncertainty sources and checks. Several other potential sources were investigated, and were found to have a negligible effect on the results and thus are not included in the total uncertainty. These include: comparing the results using only events where the photon is headed towards the A side of the detector with those where it is headed towards the C side, the number of Fourier terms used in parametrizing the LM reference distribution, the number of additional Fourier terms beyond n = 4 used in the template fit, the impact of a possible misalignment of the tracker on the charged-particle measurement, the impact of the correction for the finite tracking efficiency, and the impact of correcting for the trigger efficiency and different luminosities of the triggered event samples as a function of N rec ch on a per-event basis. While the γ η > 2.5 requirement is treated as an operational event definition which is not a source of uncertainty, the results using alternative definitions γ η > 2.0 and γ η > 3.0 were not significantly different from the nominal results. This demonstrates that the sample of photonuclear events in data is robust against small changes in their selection criteria and against detector effects in the measurement of γ η. Finally, the template method carries an implicit assumption that the modulation in the HM selection is the same as that in the LM selection. If this is not the case, the extracted v n values in the HM selection will be biased to be systematically closer to those in the LM selection. In recent measurements of azimuthal anisotropies in small systems [58], a "one-step" correction was employed using two adjacent LM reference bins to account for the possibly changing v n magnitudes with multiplicity. Since the measured v n values were observed to have no N rec ch dependence within uncertainties, this alternative method was not utilized here.
As can be seen in Table I, the systematic uncertainty is not dominated by any single source, and typically those from the pair-pseudorapidity difference, the choice of LM reference, and the event purity are codominant. The systematic uncertainties tend to be larger for the p T -dependent results than for the N rec ch -dependent results, and are typically larger for the v 3 results than for the v 2 results. For the v 2 and v 3 results as a function of N rec ch , the statistical uncertainties are typically larger than the systematic uncertainties. However, for the results as a function of p ch T , the systematic uncertainties are generally larger. This reflects, for example, the greater sensitivity to the choice of LM reference, which must be divided into different p ch T selections.

VII. RESULTS AND DISCUSSION
This section presents the measurement of the second-and third-order flow coefficients in photonuclear collisions, in both cases as a function of N rec ch and charged-particle p T . These results are compared with previous measurements of v 2 and v 3 for inclusive charged hadrons by ATLAS in inelastic, minimum-bias 13 TeV pp collisions and in 5.02 TeV p + Pb collisions [5]. Figure 15 presents the v 2 and v 3 values as a function of the event charged-particle multiplicity N rec ch . Significant, nonzero values for v 2 and v 3 are observed and they are compatible with no N rec ch dependence within uncertainties. The results are compared with previous measurements in pp and p + Pb collisions, presented here as a function of N rec ch for a slightly different p T selection, 0.5 < p T < 5 GeV. Given the shape of the p T distribution in high-multiplicity pp collisions [59] and the p T dependence of the v 2 values, the v 2 in pp collisions for 0.4 < p T < 2 GeV particles (the selection used in this measurement) may be at most 10% lower than the values presented for 0.5 < p T < 5 GeV. In photonuclear events, the v 2 values as a function of N rec ch are systematically below those in pp and p + Pb events, even after accounting for the different p T ranges in the datasets. The v 3 values are compatible with those in pp and p + Pb events within significantly larger statistical and systematic uncertainties.
The v 2 and v 3 results as a function of charged-particle p T are presented in Fig. 16 for events with 20 < N rec ch 60. The v 2 results have central values similar to those plotted as a function of N rec ch in Fig. 15 but have larger systematic and statistical uncertainties, due to the more limited size of the available data sample for the p T -dependent LM references for these measurements. For v 2 , the results trend towards negative values for p T 2.0 GeV, although with very large systematic uncertainties. It is notable that there is a statistically significant confirmation of factorization, as detailed in Sec. V B, only for v 2 with the selection 0.4 p T 2.0 GeV. Hence the values of v 2 and v 3 for p T 2 GeV have an additional caveat, arising from the fact that the magnitude of the flow signal relative to that of the nonflow contribution is significantly smaller FIG. 15. Flow coefficients v 2 and v 3 for charged particles with 0.5 < p T < 2.0 GeV in photonuclear events, reported as a function of charged-particle multiplicity N rec ch . The vertical error bars and colored boxes represent the statistical and total systematic uncertainties, respectively. The photonuclear data points are positioned at the average N rec ch value in each interval. The measurements in photonuclear events (solid symbols) are compared with those in pp collisions at 13 TeV and p + Pb collisions at 5.02 TeV [5] (open symbols), integrated over 0.5 < p T < 5.0 GeV. than at lower p T . In particular, the trend towards negative v 2 values and rising v 3 values suggests that the factorization assumption could be violated. Figure 17 shows the same data as Fig. 16, but zoomed in on the vertical axis to allow a better comparison with the analogous p T -dependent values in the pp and p + Pb measurements described above, with the selection N rec ch 60. In the region 0.4 < p T < 2 GeV the central values of the v 2 are smaller than those in pp and p + Pb collisions, similar to that observed in the N rec ch -dependent results in Fig. 15. However, due to the larger uncertainties in the p T -dependent case, the v 2 values for photonuclear and pp collisions are compatible within the uncertainties of the former in the range p T < 2 GeV. The v 3 values are compatible between systems within large uncertainties.
There are currently no published theoretical predictions for flow coefficients in photonuclear collisions within a hydrodynamic or parton transport framework. In such frameworks, the elliptic and triangular flow coefficients scale with the initial geometry eccentricities, ε 2 and ε 3 respectively, and the charged-particle multiplicity dN ch /dη. In the vector-meson dominance picture, photon-hadron interactions arise through fluctuations of the photon into hadronic states with the same quantum numbers as vector mesons, which have a nontrivial initial transverse geometry. This geometry is determined by the spectrum of these fluctuations, and while models of this spectrum exist [60], they have not yet been adapted to provide quantitative models. In the absence of a complete model, the magnitude of the eccentricities can be estimated by noting that fluctuations of the photon into light vector-meson states such as the ρ give the largest contribution to the cross section. The initial geometries for ρ + Pb collisions can be computed with a Monte Carlo Glauber calculation [61] which treats the ρ meson as having two constituent quarks. The resulting mean values of the second-and third-order spatial eccentricities, ε 2 and ε 3 , are nearly identical to those in the p + Pb case. Also, when comparing p + Pb and photonuclear events with the same N rec ch , in fact the relevant dN ch /dη is larger in the photonuclear events since the particles are distributed over a smaller pseudorapidity region. Thus, one might naively expect the flow coefficients to be similar in photonuclear events and p + Pb collisions. However, in order to compare any such calculation with data, a full modeling of the photon The v 2 data are also compared with a CGC-based theory calculation from Ref. [31]. These photonuclear data are the same as in Fig. 16 but with different y-axes ranges to allow comparison with additional data and theoretical predictions. fluctuations in the selected γ + Pb events needs to be carried out. In addition, correctly accounting for the boosted kinematics and limited acceptance using a fully three-dimensional simulation may be important.
An alternative interpretation of two-particle correlations in small collision systems involves interactions at the earliest time between gluon fields in the color glass condensate (CGC) framework [62]. Recently such calculations have described heavy-flavor hadron and quarkonia azimuthal anisotropies in p + Pb collisions [63,64], although calculations in the CGC framework fail to describe other aspects of the data, such as the charged-hadron flow coefficients in p + Pb at the LHC and small-systems collisions at RHIC [65,66]. The authors have extended these calculations to consider a color dipole interacting with a Pb nucleus either at the future Electron Ion Collider or in photonuclear collisions at the LHC [31]. The CGC calculation for photonuclear collisions is shown in Fig. 17 and is in reasonable agreement with the v 2 data within uncertainties. In these calculations, the Pb nucleus is described with a saturation scale Q 2 s = 5 GeV 2 and typical parton transverse momentum = 0.5 GeV, as used in calculations of v 2 for heavy-flavor mesons and quarkonia [63,64]. However, in the calculation for the photonuclear case, the parameter B p = 25 GeV −2 , which controls the transverse area of the interaction and thus the number of color domains from the Pb nucleus taking part in the interaction, is significantly larger FIG. 18. Comparison of results for raw Fourier coefficients v 2,2 and v 3,3 (left, without nonflow subtraction) and for nonflow subtracted coefficients v 2,2 and v 3,3 (right, with nonflow subtraction using the template method), shown in data (open points) and in DPMJET-III (filled points). The results in data and DPMJET-III are presented as functions of N rec ch and N truth ch , respectively. than in the previous calculations [63,64]. A larger transverse area generally contains more color domains and results in a smaller value of v 2 . Finally, it is interesting to compare the data with results from DPMJET-III, which has neither CGC-related correlations nor hydrodynamic or parton scattering effects. This comparison sets a baseline and provides a check of the possible impact of issues such as the LM template choice (which accounts for nonflow) selecting events with photon energies that differ from the HM sample. The photon energy is expected to be a function of N rec ch , which is true within DPMJET-III as shown in Fig. 6. If the nonflow contribution changes with photon energy, it may be that the nonflow correlation shape would not be the same in the LM and HM selections, which is an assumption of the non-flow subtraction procedure. In order to test if this effect might mimic the v 2 and v 3 signals observed, the DPMJET-III output has been run through the full analysis procedure, including the LM and HM selection and nonflow template fit. The resulting v n,n values are shown in Fig. 18 for both data and DPMJET-III. The left panel shows the values that would result from a Fourier fit to the HM correlation function, without accounting for nonflow, while the right panel shows the template-subtracted results. In DPMJET-III, the template-subtracted v n,n values are both negative, in contradistinction to the opposite sign for these values observed in data. Since they give negative values, which are impossible in a flow factorization interpretation, one cannot extract v 2 and v 3 values to compare with the data. Nevertheless, this demonstrates that, at least within DPMJET-III, there is no trivial effect in the non-flow subtraction issue which could result in a false signal.

VIII. CONCLUSION
This paper reports a measurement of long-range twoparticle correlations in high-energy photonuclear collisions. Events are selected from 1.7 nb −1 of Pb+Pb collision data at 5.02 TeV recorded by the ATLAS detector at the LHC in 2018 using a combination of criteria for the energy deposited exclusively in one side of the zero-degree calorimeter, and requirements on the sum of gaps in the photon-going direction ( γ η > 2.5). Azimuthal correlation functions are reported for pairs of charged particles separated by 2 < | η| < 5 units of pseudorapidity. A template-fitting method is used to subtract the nonflow contribution. The single-particle flow coefficients v 2 and v 3 are reported as a function of charged-particle multiplicity and single-particle p T . Significant, nonzero v 2 and v 3 values are observed in photonuclear events, indicating that particles produced in these events participate in azimuthally dependent, collective motion. The v 2 values are smaller than those reported in pp and p + Pb collisions at similar particle multiplicities. The results presented here provide new information about long-range collective behavior in an exotic collision system with different initial conditions than those in ordinary hadronic collisions. As such, they provide a new testing ground for different pictures of the physical origin of these correlations, and a starting point for understanding the systems which will be studied at the future Electron Ion Collider.

ACKNOWLEDGMENTS
We thank CERN for the very successful operation of the LHC, as well as the support staff from our institutions without whom ATLAS could not be operated efficiently. We acknowledge the support of ANPCyT