Nonperturbative transverse-momentum-dependent effects in dihadron and direct photon-hadron angular correlations in $p+p$ collisions at $\sqrt{s}=200$ GeV

Dihadron and isolated direct photon-hadron angular correlations are measured in $p$$+$$p$ collisions at $\sqrt{s}=200$ GeV. The correlations are sensitive to nonperturbative initial-state and final-state transverse momentum $k_T$ and $j_T$ in the azimuthal nearly back-to-back region $\Delta\phi\sim\pi$. In this region, transverse-momentum-dependent evolution can be studied when several different hard scales are measured. To have sensitivity to small transverse momentum scales, nonperturbative momentum widths of $p_{\rm out}$, the out-of-plane transverse momentum component perpendicular to the trigger particle, are measured. These widths are used to investigate possible effects from transverse-momentum-dependent factorization breaking. When accounting for the longitudinal momentum fraction of the away-side hadron with respect to the near-side trigger particle, the widths are found to increase with the hard scale; this is qualitatively similar to the observed behavior in Drell-Yan and semi-inclusive deep-inelastic scattering interactions. The momentum widths are also studied as a function of center-of-mass energy by comparing to previous measurements at $\sqrt{s}=510$ GeV. The nonperturbative jet widths also appear to increase with $\sqrt{s}$ at a similar $x_T$, which is qualitatively consistent to similar measurements in Drell-Yan interactions. To quantify the magnitude of any transverse-momentum-dependent factorization breaking effects, calculations will need to be performed to compare to these measurements.

widths of pout, the out-of-plane transverse momentum component perpendicular to the trigger particle, are measured. These widths are used to investigate possible effects from transverse-momentumdependent factorization breaking. When accounting for the longitudinal momentum fraction of the away-side hadron with respect to the near-side trigger particle, the widths are found to increase with the hard scale; this is qualitatively similar to the observed behavior in Drell-Yan and semi-inclusive deep-inelastic scattering interactions. The momentum widths are also studied as a function of centerof-mass energy by comparing to previous measurements at √ s = 510 GeV. The nonperturbative jet widths also appear to increase with √ s at a similar xT , which is qualitatively consistent to similar measurements in Drell-Yan interactions. To quantify the magnitude of any transverse-momentumdependent factorization breaking effects, calculations will need to be performed to compare to these measurements.

I. INTRODUCTION
Quantum chromodynamics (QCD) research has entered a period in which the focus of nucleon structure has shifted from a one-dimensional to a multidimensional picture. To take into account additional degrees of freedom besides the longitudinal momentum of partons within hadrons, the transverse-momentum-dependent (TMD) framework has been developed which also accounts for partons' transverse momentum. Rather than parton distribution functions (PDFs) and fragmentation functions (FFs) being integrated over transverse momentum degrees of freedom, TMD PDFs and TMD FFs contain explicit dependence on the nonperturbative transverse momentum of the parton within the nucleon. This dependence offers a way to describe the three dimensional momentum distribution of unpolarized partons within unpolarized hadrons, in addition to a variety of spin-spin and spin-momentum correlations when the parton and nucleon spin states are considered.
In the past decade, significant effort has been placed on measuring asymmetries that are sensitive to PDFs and FFs within the TMD framework. Semi-inclusive deepinelastic scattering (SIDIS) and Drell-Yan (DY) measurements have shown empirical evidence for nonzero spin momentum correlations in the initial state [1][2][3][4][5]. Additionally, e + e − annihilation and SIDIS measurements have shown nonzero spin-momentum correlations in the final-state hadronization process [6][7][8]. With the development of robust and theoretically interpretable jet finding algorithms, there have also been several measurements studying nonperturbative hadronization in unpolarized p+p collisions [9,10] as well as spin-momentum correlations in polarized p+p collisions [11,12]. Transverse single spin asymmetries of inclusive hadron production in p+p collisions of up to ∼40% also indicate large spinmomentum correlations within the nucleon and/or the process of hadronization [13]; however, these measurements cannot probe functions within the TMD framework directly because there is not a simultaneous measurement of a small and hard transverse momentum scale. * PHENIX Spokesperson: akiba@rcf.rhic.bnl.gov † Deceased The focus on multidimensional parton structure has brought the importance of soft gluon exchanges in hard interactions to the forefront of QCD research. In particular, the role of color exchanges due to soft gluon interactions with remnants of the hard scattering has brought to light fundamental predictions about QCD as a gauge-invariant quantum field theory. For example, the Sivers TMD PDF [14,15], which correlates the partonic transverse momentum with the nucleon spin, is predicted to exhibit modified universality when measured in the SIDIS and DY processes [16]. The underlying physical cause of this prediction is the interference between color fields in the two processes when an initial-state gluon is exchanged in DY vs. a final-state gluon in SIDIS. The Sivers function has already been measured to be nonzero in SIDIS [1]. The first measurements of the DY-like W boson and DY transverse single spin asymmetry were recently reported [17,18], and are consistent with the predicted sign change of the Sivers function. However, the uncertainties on the measurements are still large enough that a definitive conclusion cannot be drawn; more data will ultimately be necessary as the prediction is for the entire Sivers distribution as a function of the partonic longitudinal and transverse momentum x and k T .
In leading order perturbative QCD processes where a colored parton is exchanged in the hard interaction, and thus color is necessarily present in both the initial and final states, soft gluon exchanges can lead to new effects in a TMD framework similarly to the predicted modified universality of certain TMD PDFs. In hadronic collisions where a final-state hadron is measured and the process is sensitive to a small transverse momentum scale, factorization breaking has been predicted [19][20][21][22]. These processes can offer sensitivity to the non-Abelian nature of QCD. In processes where factorization is broken, the nonperturbative objects can no longer be factorized into a convolution of TMD PDFs and TMD FFs due to the complex color flows that are possible throughout the hard scattering and remnants of the collision and thus connect the initial and final state hadrons. In both cases of factorization breaking and modified universality of certain TMD PDFs, gluon exchanges with the remnants cannot be eliminated via gauge transformations. There have also been recent studies showing that factorization is broken for certain processes even at the collinear multi-loop level [23,24]; however, the focus of this work is to probe TMD factorization breaking effects.
In the past decade, the role of color in hadronic interactions has also been explored in several different types of observables in p+p collisions. The measurement of collective behavior in high multiplicity p+p collisions has prompted new studies of global observables in addition to interference effects between necessarily color connected multiple partonic interactions [25][26][27]. Measurements of dijet+jet and direct photon-jet+jet correlations at the Large Hadron Collider studying color-coherence effects have found effects from color-radiation patterns specific to hadronic collisions [28,29]. New dijet observables which rely on jet substructure have been shown to be sensitive to color flow in tt events [30]. With theoretical advances in jet substructure techniques that may be more sensitive to factorization breaking effects, such as within soft-collinear effective theory [31], future measurements involving jet substructure may be an effective way to probe soft radiation patterns that allow quantification of color flow effects similarly to Ref. [30]. The advancement of both experimental and theoretical techniques has allowed the magnitude of effects from color flow and color connections to be probed in a variety of ways across various subfields of QCD.
A previous measurement [32] in p+p collisions at √ s = 510 GeV found that nonperturbative momentum widths in dihadron and direct photon-hadron angular correlations binned in a fixed range of the away-side associated hadron p T decreased with the hard scale of the interaction, which is qualitatively opposite to what is expected from Collins-Soper-Sterman (CSS) evolution [32] which comes directly out of the derivation of TMD factorization [33]. This study is intended to extend the analysis of Ref. [32] for processes which are predicted to violate factorization; however, this measurement is made at √ s = 200 GeV. It supplants the previous PHENIX analysis at √ s = 200 GeV [34] with an approximately four times increase in integrated luminosity. The two different center-of-mass energies of Ref. [32] and this analysis allow TMD effects to be studied as a function of √ s in processes predicted to violate factorization.

II. EXPERIMENT DETAILS
In 2015 the PHENIX experiment collected data from p+p collisions at √ s = 200 GeV. A total minimum bias integrated luminosity of 60 pb −1 was used for the analysis of dihadron and direct photon-hadron angular correlations, from which data quality assurance and z vertex cuts of |z vtx | < 30 cm were applied. The PHENIX detector measures two particle angular correlations via its electromagnetic calorimeter (EMCal) and drift chamber (DC) plus pad chamber (PC) tracking system located in two central-rapidity arms. The central arms are nearly back-to-back in azimuth and each cover approximately ∆φ ∼ π/2 radians in azimuthal angle and ∆η ∼ 0.7 units in pseudorapidity centered about midrapidity. A schematic diagram of the two central arms is shown in The PHENIX EMCal [36] comprises eight sectors, four in each arm. Six of the sectors are lead-scintillator sampling calorimeters, and the other two are lead-glasš Cerenkov calorimeters. A high-energy-photon trigger in the EMCal is used to identify events with a high-p T photon. Photons are identified with a shower shape cut that removes charged hadrons as well as most high energy clusters that overlap closely with another photon, which helps eliminate π 0 merging effects at high p T . The granularity of the EMCal as well as the shower shape cuts allow for π 0 and η reconstruction up to ∼20 GeV in the diphoton channel, which allows for direct photon measurements as well. In this analysis isolated photons are measured between 5 < p T < 15 GeV/c, and neutral pions are measured between 4 < p T < 15 GeV/c. Previous π 0 , η, and direct photon cross sections as well as π 0hadron or direct photon-hadron correlation can be found in Refs. [32,34,[37][38][39].
The PHENIX tracking system [35] measures nonidentified charged hadrons with a DC and PC tracking system. The PC1, located radially behind the drift chamber, and PC3, located radially in front of the EMCal, allow for tracks identified in the DC to be matched with PC hits. The track matching condition in the PC1 and especially the PC3 reduces secondary tracks from conversions or decays. The ring-imagingČerenkov system, located radially behind the DC in Fig. 1, is also used to remove electrons from the charged hadron sample. The tracking system is also used to suppress hadronic shower contamination when identifying photons by matching tracks in the DC and PC to showers in the EMCal and subsequently removing them. In this analysis charged hadrons are collected from 0.5 < p T < 10 GeV/c. A previous nonidentified charged hadron cross section from PHENIX can be found in e.g. Ref. [40].

III. ANALYSIS
Dihadron and direct photon-hadron correlation functions are constructed following the methods of Refs. [32,34]. The correlation functions are constructed by measuring the raw number of correlated trigger particles and associated hadrons, where the trigger particle refers to either a leading π 0 or isolated photon. To account for the PHENIX acceptance, the yields are corrected by a mixedevent distribution by constructing correlations between trigger particles from one event and associated hadrons from a different event. The mixed event distribution is constructed on a run-by-run basis to quantify any changing efficiencies of the detector with time. The correlation functions are additionally corrected by a charged hadron efficiency to quantify the inefficiency of the DC and PC tracking system. The charged hadron efficiency is determined by simulating single particle hadrons in a GEANT3-based description of the PHENIX detector [41]. After corrections the correlation functions are divided by the total number of trigger particles to construct a per-trigger yield. In this analysis a 9% uncertainty is assigned to the charged hadron efficiency, which amounts to an overall normalization uncertainty on the per-trigger yields. In total the correlation functions correspond to yields within full azimuthal coverage and |η| < 0.35. The correlations are constructed similarly to previous PHENIX analyses in Refs. [32,34,[42][43][44].

A. Statistical subtraction of decay photons
The correlation functions can be constructed for any generic trigger-associated hadron combination. To identify direct photon-hadron correlations, an additional statistical subtraction must be applied to account for the background due to decay photon-hadron correlations. Reference [44] used a method which statistically subtracts the decay photon-hadron per-trigger yield component from a total inclusive photon-hadron per-trigger yield component with the following equation In this equation Y is the per-trigger yield of a particular component denoted by the subscript, and R γ is the relative contribution of direct photons to decay photons defined as R γ ≡ N inclusive /N decay . In Ref. [44], direct photons are defined as any photon not from a decay process, and includes next-to-leading order fragmentation photons.
To reduce the presence of fragmentation photons, Refs. [32,34] implemented an isolation cone criterion in the method, which we also use in this analysis. This has the added benefit of also reducing the decay photon background and thus providing a larger signal to background ratio for the direct photons. Additionally, photons that could be tagged as coming from π 0 or η decays were removed. Because the tagging procedure does not remove all decay photons, Eq. 1 was modified to include these cuts and introduce isolated photon quantities where again the trigger particles are noted in the subscripts for the various per-trigger yields and "iso" refers to both "isolated" and not tagged as coming from a π 0 decay. Here R iso γ is the relative contribution of isolated direct photons to isolated decay photons, defined in a similar way to R γ except with isolated quantities. While the subtraction procedure of Ref. [44] removes all decay photons, the modified isolated subtraction procedure removes background due to isolated decay photon-hadron correlations, which, in the PHENIX acceptance, are due most often to isolated neutral pion decays which decay asymmetrically such that the low-p T photon is not detected.
To implement the isolated statistical subtraction procedure, isolation and tagging cuts are applied at the event-by-event level in the data analysis. Inclusive photon candidates are removed from the analysis if a partner photon with p T > 500 MeV/c is found such that the invariant mass of the pair falls within the π 0 or η invariant mass regions, 0.118 < m γγ < 0.162 and 0.5 < m γγ < 0.6 GeV/c 2 . In addition to the tagging cuts, an isolation cut is implemented. The isolation criterion is that the sum of EMCal energy and p T of charged tracks within a cone radius of 0.3 radians must be less than 10% of the candidate photons energy, similar to Ref. [34]. To reduce the effects of the detector acceptance, candidate isolated photons are also required to be ∼0.1 radians from the edge of the detector in both φ and η. This forces a large portion of the cone to fall within the PHENIX acceptance.
The number of isolated decay photon-hadron pairs are not known a priori because an unknown fraction of the isolated photons still come from decay processes. The isolated decay photon-hadron correlations are determined with a Monte Carlo generated probability density function that quantifies the probability of an isolated π 0 to decay to an isolated photon within the PHENIX acceptance such that the photon was not able to be tagged as coming from a decay process. The functions are used to map the isolated π 0 -hadron correlations to the corresponding daughter isolated photon-hadron correlations for use in the statistical subtraction procedure, where the isolation criterion for the π 0 meson is the same as described above for single photons. In total a 4% systematic uncertainty was assigned for the statistical subtraction process, which includes additional background coming from higher mass state decays. With the probability functions, the decay photon-hadron per-trigger yields can be determined by where P (p π 0 T , p γ T ) is the probability density function described above and N are the number of isolated π 0hadron pairs or isolated π 0 triggers measured as indicated in the subscripts.
To determine R iso γ , the quantity R γ must be corrected for the isolation and tagging cuts. The R γ was previously measured in Ref. [34], so these values are used and corrected with the measured isolation and tagging efficiencies of this analysis. The quantity is the ratio of isolated inclusive photons to isolated decay photons, is dependent only on the photon p T , and can be written as where "niso" refers to "not isolated," tag decay is the tagging efficiency, niso decay is the isolation efficiency, and the various N values are the number of photons that correspond to the subscript and superscript with which they are associated. The bottom-most fraction in Eq. 4 is simply the number of isolated inclusive photons divided by the total number of inclusive photons and can be determined by just counting the photons that pass the necessary cuts. The decay photon tagging efficiency can be written as tag decay = R γ N tag decay /N inclusive and can also be determined by counting the number of photons that pass the cuts. To determine the isolation efficiency, the isolation cut is applied at the level of the parent π 0 and the decay probability functions are used to map the effect to the daughter photons. The isolation efficiency can be written as Therefore, all of the necessary components can be determined to calculate R iso γ . Figure 2 shows the measured values of R iso γ as a function of p γ T in this analysis compared to the previous PHENIX publication [34] at the same √ s. The values show consistency with the previous analysis and are all larger than unity, indicating the signal-to-background of the isolated direct photons to decay photons. These can also be compared to the values of R γ in Refs. [34,44] showing that the tagging and isolation cuts increase the signal-to-background.  The dihadron correlations have two peaks at ∆φ ∼ 0 and π corresponding to intrajet and back-to-back jet correlations, respectively. The isolated direct photon-hadron correlations show away-side yields that are consistently smaller than the corresponding π 0 -hadron yields, which is to be expected because the π 0 -hadron correlations probe larger hard scales due to the π 0 being a fragment of the jet. Note that a 9% charged hadron normalization uncertainty is not explicitly shown on any of figures displaying the per-trigger yields; this is denoted in each of the captions where the uncertainty applies.

B. pout distributions
The hard scattering quantity p out = p assoc T sin ∆φ, the transverse momentum component of the associated hadron with respect to the trigger particle axis, is defined kinematically in Fig. 4. Rather than constructing the per-trigger yields as a function of ∆φ, they are instead constructed as a function of p out in a similar way to the ∆φ correlations. These distributions are the quantity of interest because the nonperturbative TMD structure can be observed from the correlation functions; additionally they have the advantage that the nonperturbative com- ponent of the away-side jet width can be separated from the perturbative component in momentum space with a transition region between the two components.
Reference [32] found that the nonperturbative momentum widths of the p out distributions binned in a fixed range of p assoc T decreased with the hard scale p trig T , contrary to what is qualitatively expected from CSS evolution. The p out distributions were binned in a fixed p assoc T range only, which means that the longitudinal momentum fraction z of the away-side hadron with respect to the away-side jet was generally decreasing as the hard scale increased. Here we bin the p out correlation distributions as a function of the quantity x E , which is defined in Ref. [34] as cos ∆φ (6) and is geometrically shown in Fig. 4. Because full jet reconstruction within PHENIX severely limits available statistics due to the limited acceptance, this quantity can be used as a proxy for z, the longitudinal momentum fraction of the associated away-side hadron with respect to the away-side jet. Although x E is an approximation for z, this alternative binning allows for a clearer comparison of the p out distributions between hard scales as the associated hadrons are then binned in a similar kinematic way that is normalized by the near-side p trig T . Because x E is only a proxy for z, there is an embedded dependence on the correlation between these two quantities for both the π 0 -hadron and direct photon-hadron correlations. Previous correlations measurements from PHENIX have shown that the quantityx h =p assoc T /p trig T , where quantities with a hat indicate partonic quantities, is on average less than unity [34], and this has also been shown in direct photon-jet correlations [45]. Thus for direct-photon hadron correlations, z > x E . For π 0hadron correlations, there is an additional dependence on z T = p trig T /p trig T , which is roughly 0.6 at Relativistic-Heavy-Ion-Collider (RHIC) energies [43]. Thus, for π 0hadron correlations, on average z < x E ; therefore the dihadron and direct photon-hadron correlations are on average probing different values of the away-side hadron momentum fraction z.    The p out distributions for p+p collisions at √ s = 200 GeV binned in 0.1 < x E < 0.5 are shown in Fig. 5. The open points are the π 0 -hadron correlations, while the filled points are the direct photon-hadron correlations. In constructing the correlations, the underlying event was statistically subtracted following a method similar to Ref. [32]. The functions used to statistically subtract the underlying event are shown as fits to the away-side distributions in Fig. 3. Although the correlation functions are binned in x E instead of p assoc T , there is still a clear transition from nonperturbative to perturbative sensitivity in the distributions. This is highlighted by the Gaussian fits to the small p out region [-1.1,1.1] GeV/c in the figure, where the fits clearly fail at describing the correlations at large p out . We also note that Ref. [32] found that the large p out region was described reasonably well with a Kaplan fit; here the distributions are not described by a Kaplan fit due to the smaller center-of-mass energy. When √ s is smaller, it is less likely that a high p T gluon radiation will occur such that p out is large. This causes the p out distributions to fall more quickly towards zero at large p out .
The widths of the Gaussian fits are extracted to quantify the evolution of the nonperturbative away-side jet widths as a function of p trig T . Systematic uncertainties are evaluated by adjusting the fit region by ± 0.2 GeV/c and taking the absolute difference of the resulting Gaussian width. The values are shown in Fig. 6 and Table I and clearly demonstrate that the widths increase with p trig T . This is in contrast to Ref. [32], where the widths decreased as a function of p trig T when the p out distributions were binned in a fixed p assoc T range and thus not in a way to account for the longitudinal momentum fraction of the away-side hadron with respect to the near-side trigger particle. To study the dependence of the nonperturbative momentum widths on the fragmentation quantity x E , the p out distributions were constructed as a function of x E . To have sufficient statistical precision, and to compare to the √ s = 510 GeV data which is described later, the distributions were integrated over a larger range of 7 < p trig T < 12 GeV/c. The p out distributions are shown in Fig. 7 and the nonperturbative structure is fit with a Gaussian function, shown as dashed lines and solid lines for π 0 and direct photon triggered correlations, respectively. The nonperturbative to perturbative transition is still visible; however, the region of both large p out and x E lacks statistical precision. The Gaussian fits are performed in varying regions of p out , depending on the x E bin because the nonperturbative structure strongly depends on the x E bin probed. The Gaussian widths are extracted and shown as a function of x E in Fig. 8, where the systematic uncertainties on the widths are estimated in a similar way to the previous nonperturbative momentum widths. To compare the results measured here to the previous PHENIX data at √ s = 510 GeV, the data from Ref. [32] were rebinned in x E similarly to the results shown here. The nonperturbative momentum widths are extracted from Gaussian fits to the small p out region and the widths from both center-of-mass energies are shown as a function of p trig T in Fig. 9. Note that only the π 0hadron correlations with p trig T >7 GeV/c were reanalyzed; this is because, due to detector capabilities, there was a minimum p assoc T limit of 0.7 GeV/c in the √ s = 510 GeV analysis [32]. For the minimum x E cut of 0.1 to be unbiased, we require that the p trig T >7 GeV/c in the √ s = 510 GeV data. The values at the two different center-of-mass energies are consistent with one another at a similar p trig T within uncertainties. Figure 10 shows the nonperturbative Gaussian widths as a function of x E at the two different center-of-mass energies. The nonperturbative Gaussian widths also show little dependence as a function of x E on the center-of-mass energy. The data shown here will provide additional constraints on processes predicted to break TMD factorization and appear to show that, at the hard scales and energies probed by RHIC, the nonperturbative away-side widths are consistent within uncertainties at two different center-of-mass energies. Similar conclusions have been drawn by the STAR collaboration, where polarized TMD observables are consistent within uncertainties between √ s = 200 and 500 GeV [11,12].
The nonperturbative away-side jet widths are also shown as a function of x T = 2p trig T / √ s in Fig. 11. The Gaussian widths do not appear to scale with x T ; however, they appear to show qualitatively similar behavior to DY interactions. The nonperturbative momentum widths sensitive to a small transverse momentum scale increase with √ s at a similar x T . This behavior as a function of √ τ = x 1 x 2 , where x 1 and x 2 are the partonic momentum fractions of the quark antiquark pair, and √ s can be observed in TMD momentum widths measured from DY data (see e.g. [46]  tive momentum widths clearly rise with √ s, while in the measurements presented here, as well as other polarized TMD observables from RHIC [11,12], nonperturbative momentum widths are consistent with each other as a function of p trig T and √ s.

V. CONCLUSIONS
Dihadron and direct photon-hadron correlations were measured at the PHENIX experiment in p+p collisions at √ s = 200 GeV. These processes have been predicted to violate factorization in a TMD framework due to soft gluon exchanges with remnants that are possible in both the initial and final states in hadronic collisions [21,22]. To better interpret the possible factorization breaking effects that were studied in Ref. [32], the p out distributions  were binned in the hard scattering variable x E to provide more control over the fragmentation dependence as a function of the hard scale of the interaction. When accounting for the away-side hadron p assoc T with respect to the near-side p trig T , the nonperturbative momentum widths are found to increase with the hard scale of the interaction p trig T . This is qualitatively similar to DY interactions in that nonperturbative momentum widths sensitive to a small transverse momentum scale increase with the hard scale of the interaction. This behavior has been verified phenomenologically in both DY and SIDIS processes [47][48][49][50][51], and this measurement shows that this behavior is similar in processes which are predicted to break factorization.
For the time being, the only obvious way to quantify factorization breaking effects is to compare data in processes that are predicted to violate factorization and calculations assuming factorization holds. Such calculations are not available at this time, largely due to the fact that even for unpolarized observables TMD PDFs and TMD FFs are not yet well constrained due to a lack of data. The first global fit of TMD data was very recently reported [52], and in particular the study calls for more data covering a broader range of kinematic variables to constrain future global fits. Future comparisons of magnitudes, shapes, and evolution of these observables to calculations should be made to quantify the magnitude of factorization breaking effects; for example, the rate at which TMD observables evolve in processes predicted and not predicted to break factorization may be different.
Nonetheless, additional observables can be identified to provide more data to quantify effects from factorization breaking. For example, spin asymmetries in the direct photon-jet channel have been predicted to arise due to factorization breaking effects [53]; polarized p+A collisions have also been proposed as an avenue to quantify factorization breaking effects in the same hard scattering channel [54]. Unpolarized processes still have the potential to quantify factorization breaking effects; for example, nearly back-to-back direct photon-quarkonium production could be used to quantify factorization breaking effects as well, depending on whether or not the quarkonium state is produced in a color singlet or color octet state. Measuring processes that may be sensitive to factorization breaking effects, including those presented here, will be important for better understanding the characteristics of QCD as a non-Abelian gauge-invariant field theory. The exploration of the role of color in hadronic interactions is deeply connected to the unique properties of QCD, and future measurements will continue to constrain the magnitudes of these effects.