Introduction

The primary objective of relativistic heavy-ion collisions is to create and study quark-gluon plasma (QGP), the deconfined state of strongly interacting matter believed to have existed microseconds after the Big Bang [1-3]. Over decades of experimental and theoretical efforts, the formation of QGP under laboratory conditions has been firmly established, marking a pivotal achievement in high-energy nuclear physics. Current research focuses on characterizing QGP properties, such as its transport coefficients, equation of state, and response to extreme electromagnetic fields, and mapping the phase diagram of quantum chromodynamics (QCD) matter [4-8].

The central focus of these investigations is the probing of collective phenomena in heavy-ion collisions, including anisotropic flow, global spin polarization, and chiral magnetic effects [9-14]. These observables, measured extensively across collision energies and systems (e.g., by the STAR, ALICE, and CMS collaborations), reveal that the QGP behaves as a near-perfect fluid with remarkable vorticity, intense electromagnetic fields, and signatures of chiral symmetry restoration [15, 16]. However, a critical challenge persists: the initial collision geometry (e.g., the reaction plane and participant eccentricity) cannot be directly accessed in experiments [17, 18]. Current methods infer the geometry indirectly via final-state momentum anisotropies, inherently conflating the initial-state properties with medium-induced effects, non-flow correlations, and event-by-event fluctuations [19-21]. This introduces systematic biases that obscure the quantitative links between the QGP properties and initial conditions, underscoring the need for direct probes of the collision geometry.

Recent advances have proposed photon-induced processes in hadronic heavy-ion collisions (HHICs) as a novel pathway to access the initial geometry [22]. When relativistic nuclei collide, their strong electromagnetic fields generate quasireal photons with polarization vectors oriented perpendicular to the motion of the colliding nucleus [23-25]. These polarized photons initiate coherent processes through both QED and QCD mechanisms:

The photon polarization direction is geometrically encoded by the initial collision configuration, which enables polarization-based probes of the reaction plane. Pioneering studies have validated this hypothesis. Xiao et al. [38] first predicted quadrupole modulation in the azimuthal angle of dileptons (defined by the relative momentum direction ) from γγ processes, directly reflecting the photon polarization anisotropy. Subsequently, Wu et al. [22] proposed leveraging vector meson photoproduction to reconstruct the reaction plane using polarization-induced decay asymmetries. Recently, STAR measurements of coherently photoproduced J/ψ mesons [39] in HHICs confirmed the alignment between the J/ψ decay anisotropy and the reaction plane, cementing photon polarization as a direct probe of the initial geometry.

In this study, we employed our established QED framework [40-42] to quantify dilepton production in HHICs, focusing on their dual sensitivity to both azimuthal emission angles () and momentum orientation () relative to the reaction plane. Unlike prior studies that primarily utilized Δϕ-dependent observables, we demonstrate for the first time, that the dilepton pair’s total momentum orientation (Φ_pair) exhibits significantly enhanced sensitivity to geometric information. Our calculations at RHIC and LHC energies reveal that the geometric constraints from Φ_pair exceed those of Δϕ-based analyses in precision, while both observables provide orthogonal perspectives on the reaction plane, bypassing the systematic uncertainties inherent to traditional flow analyses. This establishes the QED-calibrated dilepton framework, simultaneously leveraging Δϕ and Φ_pair, as a multi-dimensional probe for the model-independent extraction of the initial geometry.

Methodological Framework

The Equivalent Photon Approximation (EPA) provides a computationally efficient framework for calculating the total cross sections in heavy-ion collisions through the convolution of photon fluxes with the elementary Breit-Wheeler process [23, 43]: (1)where n(ω) represents the photon flux density and σγγ is the elementary Breit-Wheeler cross-section. While the EPA accurately describes integrated cross sections, it becomes progressively inadequate for differential observables owing to its inherent neglect of photon transverse momentum correlations and polarization interference effects between scalar (σs) and pseudoscalar (σps) interaction channels.

To overcome these limitations, we employ a lowest-order QED formulation based on the external field approximation. Following Ref. [44], the electromagnetic potentials of colliding nuclei in the Lorentz gauge are (2)where q_1,2 are the equivalent photon four-momenta, are four-velocities in the collider frame,γ and β are the Lorentz factor and velocity of colliding nuclei, b is the impact parameter of the collision, and f(-q²) denotes the nuclear charge form factor obtained via Fourier transform of the charge density distribution. The charge density of a symmetrical nucleus is characterized by the Woods-Saxon profile as follows: (3)where R_WS represents the nuclear radius and d denotes the surface diffuseness parameter, both of which are determined from electron scattering measurements [45, 46]. The normalization constant ρ₀ ensures .

With the direct and cross Feynman diagrams of the lowest-order two-photon interaction for lepton pair creation, the matrix element can be expressed as [47] (4)where (5)with . The four-momenta of the produced leptons are denoted by p₊ and p_-, while the longitudinal components of the quasi-real photon four-momenta q₁ are constrained through the relations: (6) (7)where is the lepton energy, and m is the lepton mass. The probability of the lowest order pair production can then be expressed as: (8)with (9)where F(Ni) represents the photon propagators that describe the interaction with the Coulomb field of both nuclei [42, 48, 49]: (10)and Δq is the four-momentum difference between q₁ and , with propagator assignments organized as

We construct angular correlations as follows: (11) (12)where ϕb denotes the azimuthal angle of impact parameter b. The differential distribution is obtained by integrating over: (13)where maps the momentum combinations to the angular differences derived from Eqs. 11 and 12, contains the Jacobian factors from the coordinate transformations, and dΩ_reduce denotes integration over the non-angular momentum components. Similarly, can be obtained by following the procedure outlined in Eq. 13. High-dimensional phase space integration was performed using the Vegas adaptive Monte Carlo algorithm [50] within the ROOT framework [51].

Although the QED formulation provides precise calculations of angular correlations, the dynamical mechanism through which photon polarization induces angular asymmetry (with respect to the impact parameter) remains implicit. To unveil the physical origins of this connection, we traced how the field geometry (governing polarization directions) dynamically “locks” onto the momentum-space characteristics of the produced particles. Consider the electromagnetic field configuration in peripheral heavy-ion collisions: quasi-real photons carry linear polarization determined by classical electromagnetic fields. Specifically, the electric field vector at a transverse position (defined with respect to the nucleus center as the origin) satisfies the radial constraint:this indicates that the electric field direction is radially oriented away from the nucleus center, as illustrated by the field line distribution in Fig. 1. This spatial polarization pattern becomes imprinted on the photon momentum distribution through the Fourier relation , leading to the momentum-polarization locking (14)In the dilepton production process , the scattering amplitude depends critically on the relative orientation of the photon polarization: (15)The polarization overlap reaches maximum when the two photons’ transverse momenta are collinear (). As illustrated in Fig. 1, this collinear condition forces the pair orientation Φ_pair to align with the impact parameter. The observed Δϕ asymmetry stems from the geometric suppression of the polarization correlations owing to angular misalignment. Specifically, as shown in Fig. 1, when the transverse momenta and acquire non-collinear components (induced by a finite impact parameter b), their polarization vectors and become mismatched. These misalignments introduce a characteristic modulation in the angular distribution, dominated by and terms with respect to the impact parameter axis, which are direct consequences of the spin-1 nature of photons and spin-1/2 nature of leptons in the production process [25].

Fig. 1Full-screen View

(Color online) Electromagnetic field configuration in transverse plane with collision impact parameter b

Results

Figure 2 illustrates the azimuthal modulations of dimuon pairs with respect to the reaction plane in peripheral Au+Au collisions at GeV (RHIC) with an impact parameter b = 11 fm. The detailed fiducial cuts applied in the calculations are presented in Table 1. The Δϕ distribution exhibits a dominant pattern owing to the spin of the interacting photons. A secondary modulation, originating from lepton mass effects, becomes visible as the finite muon mass breaks the helicity conservation. In contrast, the Φ_pair distribution shows a pronounced modulation, which directly reflects the alignment of the dilepton total momentum with the reaction plane. This alignment arises from the momentum-polarization locking effect discussed in Sect. 2, where the global polarization geometry encoded in the electromagnetic fields is imprinted onto the angular correlation of the dilepton pairs.

Fig. 2Full-screen View

The periodic oscillation pattern in the alignment of the and related to the reaction plane ϕb for the dimuon pair production in Au+Au peripheral collisions at GeV with impact parameter b=11 fm

The invariant mass dependence of anisotropic coefficients and is shown in Fig. 3. At low mass (Mμμ < 0.8 GeV/c²), the component is visible owing to the violation of helicity conservation. This coefficient diminishes as Mμμ increases beyond 0.8 GeV/c², and is suppressed by the restoration of helicity conservation. The coefficient displays a non-monotonic trend, initially decreasing to a local minimum at GeV/c² before rising again at higher mass. This behavior reflects competing contributions from two physical mechanisms: geometric depolarization effects dominate at intermediate mass, while polarization coherence strengthens at larger Mμμ as dileptons are preferentially produced near the nuclear periphery. In contrast, the coefficient retains a consistently positive non-zero value and increases monotonically with the dimuon invariant mass. This distinct mass dependence highlights its enhanced sensitivity as a global geometric observable to local kinematic configurations compared to .

Fig. 3Full-screen View

Second- and fourth-order modulation coefficients for Δϕ and Φ_pair as a function of dimuon invariant mass in Au+Au peripheral collisions at GeV with impact parameter b=11 fm

Figure 4 shows the rapidity dependence of these coefficients for 0.4 < Mμμ < 2.6 GeV/c². The magnitude of increases with |y|, driven by the enhanced polarization alignment in the forward/backward regions, where photoproduction occurs closer to the nuclear surfaces. This spatial localization amplifies the correlation between the photon polarization vectors and the impact parameter orientation. However, within the rapidity range |y| < 0.8, the coefficient remains approximately constant and exhibits a negligible rapidity dependence. This distinct behavior arises because the total momentum orientation of dilepton pairs depends primarily on the vector sum of photon momenta, which is a global property governed by the reaction plane geometry rather than local rapidity-dependent dynamics.

Fig. 4Full-screen View

Second- and fourth-order modulation coefficients for Δϕ and Φ_pair as a function of dimuon pair rapidity in Au+Au peripheral collisions at GeV with impact parameter b=11 fm

The impact parameter dependence of the modulation coefficients calculated for both the RHIC and LHC collision systems is presented in Fig. 5. Central collisions () yield vanishing anisotropies as the polarization geometry becomes azimuthally asymmetrical. The magnitudes of the coefficients grow with increasing b, reaching maxima near (RA is the nuclear radius). Beyond this point, and stabilize and slightly decrease owing to the depletion of photon flux coherence at ultra-peripheral separations. The impact parameter dependence of behaves similarly to Xiao’s calculations [38], despite the distinct pair p_T ranges. Remarkably, the values consistently exceed those of across the entire b range, demonstrating the superior sensitivity of pair momentum orientation to collision geometry. This systematic trend persists in Pb+Pb collisions at TeV (LHC), confirming the universality of QED-dominated polarization effects across various collision energies.

Fig. 5Full-screen View

Second- and fourth-order modulation coefficients for Δϕ and Φ_pair as a function of impact parameter in Au+Au collisions at GeV (RHIC) and Pb+Pb collisions at TeV (LHC)

Collectively, these results establish Φ_pair correlations as a precision tool for the initial geometry determination. The distinct dependence of on mass and rapidity can effectively mitigate the systematic bias induced by final-state interactions, while its strong b-dependence enables direct mapping between observable angular modulations and collision geometry. Complementary information from Δϕ harmonics provides cross-checks on polarization purity, which is crucial for separating two-photon processes from competing hadronic backgrounds. The combined analysis of both observables opens a new pathway for the model-independent extraction of the reaction plane in heavy-ion collisions, free from assumptions about medium evolution or hadronization dynamics.

Summary

In summary, we employed a QED-based approach to systematically investigate how photon-induced processes can constrain the initial collision geometry in heavy-ion collisions at RHIC and LHC energies. This work introduces a novel observable, the total momentum orientation of the dilepton pair Φ_pair, which demonstrates significantly enhanced sensitivity to geometric features, exceeding Δϕ-based analyses in precision. By combining the complementary sensitivities of the photo-produced dilepton to primordial geometry via both Δϕ and Φ_pair, our calculations establish the QED-calibrated dilepton framework as a robust, multi-dimensional probe that enables model-independent reconstruction of the collision geometry, offers orthogonal constraints on reaction plane determination, and bypasses systematic uncertainties inherent to traditional flow methods. Future measurements of these observables are expected to provide direct and unambiguous experimental constraints on primordial geometry in relativistic heavy-ion collisions.