Phenomenological constraints on $A_N$ in $p^\uparrow p\to \pi\, X$ from Lorentz invariance relations

We present a new analysis of $A_N$ in $p^\uparrow p\to \pi\, X$ within the collinear twist-3 factorization formalism. We incorporate recently derived Lorentz invariance relations into our calculation and focus on input from the kinematical twist-3 functions, which are weighted integrals of transverse momentum dependent (TMD) functions. In particular, we use the latest extractions of the Sivers and Collins functions with TMD evolution to compute certain terms in $A_N$. Consequently, we are able to constrain the remaining contributions from the lesser known dynamical twist-3 correlators.

Schematically, one writes the (polarized) differential cross section for A ↑ + B → C + X as where S T is the transverse spin vector of hadron A, f a/A(t) is the twist-t parton distribution function (PDF) associated with parton a in hadron A (similarly for f b/B(t) ), while D C/c(t) is the twist-t fragmentation function (FF) associated with hadron C in parton c. The twist-3 correlators can either be of the 2-parton or 3-parton type and are categorized into intrinsic, kinematical, and dynamical functions [1,26]. The intrinsic functions are twist-3 Dirac projections of collinear 2-parton correlators, while the kinematical functions are first transverse momentum moments of transverse momentum dependent (TMD) 2-parton functions. The dynamical functions are 3-parton correlators. The factors H, H , and H are the hard parts for each term, and the symbol ⊗ denotes convolutions in the appropriate parton momentum fractions. One then can calculate the TSSA A N as 1 (z) (i.e., the first p ⊥ -moment of Collins function H ⊥ 1 (z, z 2 p 2 ⊥ ) [57] 2 ) the intrinsic function H(z), and the dynamical function 3Ĥ FU (z, z 1 ). Later, a phenomenological analysis of A N incorporated this twist-3 fragmentation piece along with the QS term [59]. The procedure employed by the authors of Ref. [59] was (1) to use a QCD equation of motion relation (EOMR) to write H(z) in terms of H ⊥(1) 1 (z) and an integral ofĤ FU (z, z 1 ); (2) then fix H ⊥(1) 1 (z) in terms of the Collins function and F FT (x, x) in terms of the Sivers function, both of which had been previously extracted from TMD processes using a simple Gaussian ansatz for the respective transverse momentum dependence [60,61]; (3) finally, fitĤ FU (z, z 1 ) to the RHIC data. The result was a successful description of A N data from RHIC for both neutral and charged pion production, where the twist-3 fragmentation piece dominated the asymmetry [59]. This was a crucial step towards clarifying the mechanism behind TSSAs in hadronic collisions.
Since that study, there has been a recent development, due to the work in Ref. [26], that has provided an additional operator constraint, a Lorentz invariance relation (LIR), on the twist-3 FFs H ⊥(1) 1 (z), H(z), and an integral ofĤ FU (z, z 1 ) that enter A N , which was not known at the time of the computation in Ref. [59]. From a purely theoretical standpoint, the work in Ref. [26] has elucidated that multi-parton (dynamical) correlators are the fundamental objects that cause transverse spin observables. For example, one can solve the relevant LIR and EOMR and write H ⊥(1) 1 (z) and H(z) in terms of integrals ofĤ FU (z, z 1 ). That is,Ĥ FU (z, z 1 ) is the "base" function from which H ⊥(1) 1 (z) and H(z) are derived. However, from a practical phenomenological standpoint, one can use LIRs to eliminate "unknown" functions in terms of "known" functions. To be more specific, the only twist-3 functions that enter A N in p ↑ p → π X that we a priori have any information on are the kinematical correlators because they are connected to TMD functions. 4 Therefore, we can use these relations to re-express A N in terms of the maximum number of kinematical correlators. 5 To this end, we present a new analysis of the TSSA in p ↑ p → π X, where we compute the QS term and the twist-3 fragmentation piece as was done in [59]. However, we do not aim at a fit of the 3-parton FFĤ FU (z, z 1 ). Rather, we first employ the LIR and EOMR to rewrite the cross section, eliminating unknown functions where possible. We then determine the contribution only from the kinematical correlators, where we use the latest TMD evolved extractions of those functions. Consequently, we are able to provide a constraint on the parts that remain from the lesser known dynamical functions. The Letter is organized as follows: in Sec. 2 we review the relevant analytical formulae needed in our analysis, and, in particular, rewrite the fragmentation term using the aforementioned EOMR and LIR. Next, in Sec. 3 we conduct our numerical study of A N and compare our results to experimental data. Finally, in Sec. 4 we summarize our work.
2. The Qiu-Sterman and fragmentation contributions to p ↑ p → π X We consider TSSAs in the single-inclusive production of pions from proton-proton collisions, where we have indicated the momenta and polarizations of the particles. The spin-averaged differential cross section dσ in Eq.
(2) at leading-order can be written as where i is a sum over all partonic interaction channels, a can be a quark, anti-quark, or gluon and likewise for b, c, α s is the strong coupling constant, and f 1 (x) (D 1 (z)) is the standard twist-2 unpolarized PDF (FF). We have made explicit that parton c fragments into a pion. The well-known hard factors for the unpolarized cross section are denoted by S i U [62,63] and can be found in, e.g., Appendix A of Ref. [9]. They are functions of the partonic Mandelstam variablesŝ = xx S ,t = xT/z, andû = x U/z, where S = (P + P ) 2 , T = (P − P h ) 2 , and U = (P − P h ) 2 .
Let us now turn to the spin-dependent differential cross section d∆σ(S T ) in Eq. (2). All three terms in Eq. (1) enter into the analysis. However, as stated in Section 1, we will focus on the qgq SGP (QS) piece of the first term and the third (fragmentation) term: The definitions of the relevant functions can be found in Refs. [1,26]. First, we give the expression for the QS term, which reads [7,9] E h d∆σ QS (S T ) where the Levi-Civita tensor is defined with 0123 = +1, and the hard factors are denoted by S i F FT and can be found in Appendix A of Ref. [9]. There is an operator identity that relates the QS function F q FT (x, x) to the first k T -moment of the Sivers function [45], where In Eq. (7) we have indicated that the Sivers function is either the one extracted from semi-inclusive deep-inelastic scattering (SIDIS) or the Drell-Yan (DY) process [64,65]. Next, we look at the fragmentation term, which was first fully calculated in Ref. [22] and reads where M h is the pion mass, and h 1 (x) is the standard twist-2 transversity PDF. The functions H ⊥(1) 1 (z), H(z), and H FU (z, z 1 ) are, respectively, the kinematical, intrinsic, and dynamical unpolarized twist-3 FFs discussed in Section 1. The hard factors associated with them are represented by S i with the corresponding subscript, and they can be found in Appendix A of Ref. [22]. 6 Note that H ⊥(1) 1 (z) is the first p ⊥ -moment of the Collins function, These collinear twist-3 FFs are related to each other through an EOMR [17,22,23,[66][67][68], It is important to mention that there is a similar EOMR involving TMD FFs [68,69], where Therefore, comparing Eqs. (11) and (13) one can make the identification 8 (using the variable substitution y = z/z 1 ) That is, the r.h.s. of Eq. (14) is actually a TMD FF integrated over p ⊥ . We will revisit this observation below.
Another set of formulae relating the intrinsic, kinematical, and dynamical functions are LIRs [10,26,[72][73][74]. For twist-3 FFs, these expressions were derived for the first time in Ref. [26] using identities among non-local operators and taking into account constraints from Lorentz invariance. For our purposes, we focus on the one relating H ⊥(1) 1 (z), H(z), andĤ FU (z, z 1 ) [26]: We now employ both the EOMR (11) and LIR (15) to re-express the fragmentation contribution to A N in Eq. (9) in terms of the kinematical function H ⊥(1) 1 (z) and the the dynamical correlatorĤ FU (z, z 1 ) (viaH(z)), where we have used the identity (14) for the last term in brackets. The hard factors are given bỹ 6 Note that in Ref. [22],Ĥ(z) ≡ H ⊥(1) 1 (z). 7 Note that D(z) = z 2 d 2 p ⊥ D(z, z 2 p 2 ⊥ ) for a generic FF. 8 This connection was made previously in Ref. [70], where H(z, z 2 p 2 ⊥ ) andH(z, z 2 p 2 ⊥ ) were calculated in a spectator model, and in Ref. [71] in the context of the asymmetry A sin φ S UT in SIDIS. which explicitly read S qq→qq S qq→qq Sq q→qq S qq →qq Note that the factor 1/(−x t − xû) in Eq. (17) is cancelled in each of the above expressions (18)-(23). Such a simplification also occurred for A N in p ↑ → π X after using the LIR (15) [26] and appears to be a non-trivial cross-check of the result. We also mention that one can obtain the channels involving anti-quark fragmentation by charge conjugating the above partonic processes ab → cd. The hard parts are the same as the ones above, i.e.,S¯a¯b →cd =S ab→cd . At this point, we emphasize again thatH(z) in Eq. (16) is the twist-3 qgq TMD FFH(z, z 2 p 2 ⊥ ) integrated over p ⊥ . Thus, we stress that one could obtain some information onĤ FU (z, z 1 ) (viaH(z) orH(z, z 2 p 2 ⊥ )) through transverse spin observables within the collinear twist-3 and/or TMD factorization formalisms. For example, the processes p ↑ p → πX and p ↑ → π X depend onH(z). On the other hand, at small-q T , where q is the momentum of the virtual photon, the unintegrated correlatorH(z, z 2 p 2 ⊥ ) enters (along with several other terms) the asymmetries in SIDIS [68] as well as the analogous ones in e + e − → h a h b X [69]. At the same time, a measurement of A sin φ S UT in SIDIS integrated over P hT is sensitive to a single term 9 that couplesH(z) to the transversity PDF h 1 (x) [68]. 10 We want to reiterate once more thatĤ FU (z, z 1 ) (viaH(z) orH(z, z 2 p 2 ⊥ )) does show up in several different observables that can provide more information on this function in the future. In the next section, we perform an analysis of the TSSA for p ↑ p → πX and comment on what constraints onH(z) it can provide.

Phenomenology
We now give our numerical estimate for A N . Using Eqs. (4), (6), (16), we are able to express, respectively, the numerator and denominator of A N in Eq. (2) as where z min = −(T + U)/S , x min = −(U/z)/(T/z + S ), and x = −(xT/z)/(xS + U/z). In Eq. (24) the quantities H π/c,i (x, x , z) and F a,i (x, x , z) are We focus on the contributions in Eqs. (26), (27) from the functions F FT (x, x) and H ⊥(1) 1 (z), for which we have information on from the Sivers function and Collins function, respectively, that have been extracted from TMD processes [60,61,[76][77][78][79][80][81][82]. In our analysis we will ignoreH(z) in Eq. (26), which is equivalent to approximating Eq. (13) as This is known as a Wandzura-Wilczek-type approximation [83], which for TMDs was explored in Refs. [84,85]. We emphasize that Eq. (28) is a statement thatH(z)/z is parametrically smaller than H ⊥(1) 1 (z), not thatĤ (z, z 1 ) is zero. In fact,Ĥ FU (z, z 1 ) must be nonzero because it was shown in Ref. [26] that H ⊥(1) 1 (z) is an integral ofĤ FU (z, z 1 ) . That is,Ĥ FU (z, z 1 ) = 0 implies H ⊥(1) 1 (z) = 0, and consequently, Eq. (26) would vanish identically. Moreover, we know from current extractions of the Collins function that H ⊥(1) Here, the purpose of our computation is not to offer a complete analysis of A N ; it is to use recent TMD evolved extractions of known (i.e., kinematical) inputs to the observable, along with a new constraint from the LIR (15), to assess how well we are currently able to describe the data and ascertain what contributions remain from the dynamical functions. This will help guide a future fit of these correlators, in particularĤ FU (z, z 1 ) (orH(z)), where one would be able to confirm or refute the approximation (28). The functionĤ FU (z, z 1 ) was originally extracted in Ref. [59] before the LIR (15) was derived, and, therefore, that work must be updated to include this constraint.
As mentioned above, we will ignoreH(z) and compute the terms in Eqs. (26), (27) that involve F FT (x, x) and H ⊥(1) 1 (z), using the latest fits of the Sivers and Collins functions that incorporate TMD evolution. We elaborate a bit here. TMD evolution is conventionally carried out in coordinate space, and thus one usually defines the Fourier transform of the TMD in a two-dimensional b-space [86]. Take the quark Sivers function as an example [80,87], which is probed at the momentum scale Q. TMD evolution evolves the above Sivers function from an initial low scale Q 0 to a high Q: f ⊥,q(α) 1T (x, b; Q). At the same time, one can further perform an operator product expansion to relate the TMD at the initial scale Q 0 to a corresponding collinear function [87]. For the Sivers function at small b, one has a connection to the QS function: where the coefficient C was derived in [88,89]. Thus, in the global analysis of the Sivers asymmetry that incorporates TMD evolution, one can directly extract the QS function F FT (x, x). Similarly we can directly extract H ⊥(1) 1 (z) through an analysis of the Collins asymmetry with TMD evolution. For details, see Ref. [82].
In the following numerical calculations, we use H ⊥(1) 1 (z) extracted in Ref. [82], evolved through the "diagonal" piece of its evolution equations [16,90]. We found that the QS function extracted in Ref. [80] leads to a negligible contribution from Eq. (27) to A N , so we will ignore it in the plots below. Instead, we focus on the fragmentation term in Eq. (24) driven by (26), which gives the dominant contribution to the asymmetry. This confirms the original findings in Ref. [59]. A few additional comments are in order about the size of the QS term. The smallness of this term is because of the large-x behavior of the Sivers function [60,61,91], which can be explored and constrained in future experiments, like those at the 12 GeV upgrade of Jefferson Lab (JLab 12) [92,93]. More specifically, one can only generate a non-negligible QS term if the β parameters (that control the large-x behavior of the Sivers/QS function) are different for up and down quarks [59,91], whereas in Ref. [80] one has β u = β d . While there does exist an extraction of the Sivers function where β u β d [61], this extraction is not at the NLL accuracy that we use for other functions in our analysis. Moreover, if such a Sivers function is used, as explored in Ref. [59], it most likely will not overwhelm the large effect from the fragmentation term. Nevertheless, it will be interesting to include both the QS functions and the fragmentation terms in the global analysis in order to understand A N fully. Such a study will be explored in the future.  (28)) compared to (a) STAR Collaboration 2004 data [33] on A N for π 0 , (b) STAR Collaboration 2012 data [36] on A N for π 0 at y = 3.68, (c) STAR Collaboration 2008 data [39] on A N for π 0 at y = 3.3, (d) STAR Collaboration 2008 data [39] on A N for π 0 at y = 3.7, (e) BRAHMS Collaboration 2007 data [35] on A N for π ± (black closed data π + , blue open data π − ) at θ = 2.3 o , (f) BRAHMS Collaboration 2007 data [35] on A N for π ± (black closed data π + , blue open data π − ) at θ = 4 o . The solid lines correspond to the results using the central parameters from Ref. [82] for h 1 (x) and H ⊥(1) 1 (z). The shaded regions correspond to an estimate of 90% C.L. error band from Ref. [82] due to uncertainties in h 1 (x) and H ⊥(1) Predictions for A N in π ± production as function of x F at y = 3.7 for positive x F using Eq. (16) (along with the approximation (28)). The error bands correspond to the theoretical uncertainty from both h 1 (x) and H ⊥(1) 1 (z).
In Fig. 1, we show the result of our calculation compared with the BRAHMS charged pion and STAR neutral pion data for A N as a function of x F [33,35,36,39]. The error bands for A N in those plots are based on the uncertainties in h 1 (x) and H ⊥(1) 1 (z) from the fit in Ref. [82]. In the large-x F (i.e., large-x and large-z) region these functions are not well-constrained, and consequently, one obtains large errors due to these inputs. We see that, although the theoretical calculations undershoot or overshoot A N in some places, the central curves do a reasonable job in describing the data. We are especially encouraged by these plots given that the contribution fromH(z) still needs to be included. This clearly demonstrates that this function must be nonzero. Through this computation, we now have a constraint onH(z) and leave a fit of this function to A N data for future work. We emphasize again that the unintegrated version of this correlator also enters multiple asymmetries in SIDIS and e + e − → h a h b X, whileH(z) itself can be directly measured in A sin φ S UT in SIDIS integrated over P hT . Moreover, since F FT (x, x) and H ⊥(1) 1 (z), h 1 (x) enter the TMD evolution equations for the Sivers and Collins asymmetries, respectively, in SIDIS and e + e − → h a h b X, one can eventually perform a global analysis that includes all these observables along with A N in proton-proton and lepton-nucleon collisions (where F FT (x, x),H(z), H ⊥(1) 1 (z), h 1 (x) all enter). This would better constrain the large-x F behavior of these functions and greatly reduce the error bands in our plots since we have data from RHIC in this region. We found that the uncertainty in h 1 (x) in this regime is what dominates the error over that from H ⊥(1) 1 (z), see Fig. 2(a). Thus, it is evident that the A N data would allow us to drastically improve the extraction of transversity. Also, future measurements at JLab12 can improve the situation in the large-x region [93]. In order to demonstrate the powerful capability of RHIC future measurements [2], we present our predictions for A N in π ± production at 500 GeV in Fig. 2(b). One can clearly see that large-x F measurements of A N will reduce the uncertainty of the large-x behavior of transversity and, together with other data sets, allow us to explore the missing contribution fromH(z). In addition, we also give our result for A N as a function of P hT in Fig. 3 compared with the STAR data from Ref. [94]. One can see that our calculations exhibit a flat behavior, similar to that shown in Ref. [59]. The reason is that in the forward region, wheret becomes very small, the qg → qg channel dominates, and the hard function S qg→qg H ⊥ 1 ∝ 1/t 3 compensates the twist-3 (P hT ) −1 fall off of the asymmetry. Again, one has to keep in mind that there is still a term missing,H(z) from our analysis, which needs to be fit to data as has been emphasized above. It is also important to emphasize that the experimental data has a very large uncertainty which prevents an unambiguous identification of the P hT -dependence. These open issues can only be addressed by future experimental measurements and theoretical work.
We end this section with a brief comment about the fragmentation contribution to A N in p ↑ A → π 0 X. Recently, a calculation of this term was carried out in Ref. [95] that included gluon saturation effects in the unpolarized nucleus. The authors found that the first two terms in braces in Eq. (9) are proportional to A −1/3 (see also [96]), while the third term is proportional to A 0 . Since in Ref. [59] one finds that this third term is negligible (see Fig. 3 of Ref. [59]), the  (28)) compared to STAR Collaboration 2013 data [94] on A N for π 0 at √ s = 500 GeV as a function of P hT (a) 0.16 < x F < 0.24, (b) 0.24 < x F < 0.32, and (c) 0.32 < x F < 0.4. The solid lines correspond to results using the central parameters from Ref. [82] for h 1 (x) and H ⊥(1) 1 (z), we can obtain a new estimate for the contribution of the third term in Eq. (9) to A N . To be specific, we replace H π/c,i (x, x , z) in Eq. (24) with where, as before, S iĤ FU can be found in Appendix A of Ref. [22], and we include only terms involving H ⊥(1),π/c 1 (z). 11 Our estimate is shown in Fig. 4. We see that the contribution to A N from the third term in Eq. (9) is actually moderate in size and certainly not negligible. Since this part of A N in p ↑ A collisions is proportional to A 0 [95], the fragmentation term for A N is not inconsistent with the STAR data [97] on the asymmetry in p ↑ A → π 0 X. Recently, the author of Ref. [98] has found that, in contrast to Ref. [99], the QS piece contribution to A N in pA collisions is small (and may very well vanish), leaving the fragmentation term as the only source of A N in pA collisions 12 .

Summary
We have given an estimate of A N in p ↑ p → π X due to the kinematical twist-3 functions, which are first (k T or p ⊥ ) moments of TMD functions. Using the newly derived LIR (15) from Ref. [26], we have written A N in terms of a maximum number of these correlators, which helps optimize the phenomenology. In particular, we computed the terms in A N that involve H ⊥(1) 1 (z), using the most recent TMD evolved extraction of that function [82] and confirmed the conclusion of Ref. [59] that the fragmentation piece can be the dominant source of the asymmetry. Note that the dynamical twist-3 FFĤ FU (z, z 1 ), which played a key role in the analysis of Ref. [59], is also crucial to our work here since H ⊥(1) 1 (z) (andH(z)) is an integral of this function. We have simply "re-shuffled" the pieces by incorporating LIRs (and EOMRs) in order to write the cross section in terms of a maximal number of known functions. We also found that the existing extractions of the QS function lead to a negligible contribution to A N , although one must be mindful that this function is not well-constrained at large x. Our results give a reasonable description of the data, although the central curves cannot fully account for all of the asymmetry. However, we know the pieces that remain from the dynamical twist-3 functions, in particularĤ FU (z, z 1 ) (viaH(z)), must be nonzero. Therefore, we must fit this function to A N in order to obtain a complete result for this observable. In addition, we can reduce the rather large theoretical uncertainties by including both the transversity and Sivers functions simultaneously in such a fit. We leave both of these matters for future work, which are needed to definitively resolve the puzzle of what causes A N .