Recursive model for the fragmentation of polarized quarks

We present a model for Monte Carlo simulation of the fragmentation of a polarized quark. The model is based on string dynamics and the ${}^3P_0$ mechanism of quark pair creation at string breaking. The fragmentation is treated as a recursive process, where the splitting function of the subprocess $q \to h + q'$ depends on the spin density matrix of the quark $q$. The ${}^3P_0$ mechanism is parametrized by a complex mass parameter $\mu$, the imaginary part of which is responsible for single spin asymmetries. The model has been implemented in a Monte Carlo program to simulate jets made of pseudoscalar mesons. Results for single hadron and hadron pair transverse-spin asymmetries are found to be in agreement with experimental data from SIDIS and $e^+e^-$ annihilation. The model predictions on the jet-handedness are also discussed.


I. INTRODUCTION
The fragmentation of quarks into hadron jets is a nonperturbative QCD process and as such is still poorly understood theoretically. For this reason many simulation models based on the recursive splitting process q → h þ q 0 [1,2] have been developed to account for the main jet properties, prepare high-energy experiments and analyze their results. The most elaborate one is the Lund symmetric model (LSM) [3], based on the semiclassical dynamics of a string or color flux tube. This model is the basis of the fragmentation part in the event generator PYTHIA [4], which successfully describes the experimental data for different scattering processes. Up to now, however, the available codes of quark fragmentation do not include the quark spin degree of freedom, at least in a systematic way, like they do for flavor. In particular they cannot simulate the Collins effect [5] and the dihadron transverse-spin asymmetries observed in semiinclusive deep inelastic scattering (SIDIS) and in e þ e − annihilation [6]. Progresses toward the inclusion of the quark spin has been done in Refs. [7][8][9], where it was proposed to extend the LSM by assuming that the quarkantiquark pairs created during string breaking are in the 3 P 0 state, which means that their spins are parallel and they have one unit of orbital angular momentum, antiparallel to the total spin. The model of Refs. [7][8][9] is a quantum version of the classical string þ 3 P 0 mechanism proposed in order to explain the inclusive hyperon polarization [3] and used to explain the single-spin asymmetry in p↑ þ p → π þ X [10].
The road map proposed in Refs. [7][8][9] has now been pursued, and in this paper we give full account of the model which has been developed and of the program we have written to simulate the fragmentation of a polarized quark. The main results obtained will also be given, as well as their comparison with the data. Preliminary results were presented in [11,12].
The paper is organised as follows. The theoretical framework of the recursive model on which the code is based is summarized in Sec. II. The details of the string model, first in the spinless case, and then introducing the spin and the 3 P 0 model, are spelled out in Sec. III. The implementation of the model in a Monte Carlo code is the subject of Sec. IV, while the results of the simulation are compared with corresponding data from the COMPASS and BELLE Collaborations in Sec. V. Section VI is dedicated to the simulation of jet-handness, and our conclusions are drawn in Sec. VII.
For completeness, we remind that a different model of polarized quark fragmentation, based on the spin density matrix formalism within the extended quark-jet framework of Field and Feynman, has been recently proposed in Ref. [13].

II. THE GENERAL RECURSIVE MODEL
We consider the hadronization occurring in e þ e − annihilation or W AE decay in two jets. In SIDIS the antiquarkq B is replaced by the target remnant represented by a diquark. This process is supposed to occur via a quark chain diagram shown in Fig. 1 and modeled as the set of splittings i.e., as the iteration of the elementary splitting where the flavor content of the hadron h is qq 0 . The index r ¼ 1; 2…N in Eq. (2) is the rank of the hadron or of the splitting quark. The production of baryons is not included in the present code. In the following k denotes the 4-momentum of a quark, p that of a hadron. In Eq. (3) momentum conservation gives p ¼ k − k 0 . In the recursive model one assumes that the initial 4-momenta k 1 ≡ k A and kB ≡ −k B are on mass shell and generated beforehand. 1 In the q AqB centre-of-mass frame we orient theẑ axis (named jet axis) along k A . In SIDIS this axis usually differs from that defined in the laboratory frame by the virtual photon momentum, due to the primordial transverse momentum of the struck quark inside the nucleon. The "light cone" components of p are p AE ¼ p 0 AE p z and the transverse ones p T ¼ ðp x ; p y Þ (similarly for k and k 0 ). The mass-shell constraint writes p þ p − ¼ m 2 h þ p T 2 ≡ ϵ 2 h , where ϵ h is the hadron transverse energy.
The energy-momentum sharing between h and q 0 in Eq. (3) is drawn at random following the splitting distribution where the longitudinal splitting variable Z ¼ p þ =k þ is the fraction of forward light cone momentum of q taken by the hadron h. ðdZ=ZÞd 2 p T ¼ d 3 p=p 0 is relativistically invariant.
III. THE STRING + 3 P 0 MODEL A. Review of the spinless string fragmentation model Hadronization of a quark pair q AqB is considered as successive breakings of a massive string stretching between q A andq B , which we call here a sting. Each breaking creates a new quark-antiquark pair. A semiclassical treatment of this process leads to a recursive model with a very specific form of the splitting function. We will start with the simple classical ð1 þ 1ÞD yoyo model [14] where the created quarks have no mass, no spin and no transverse momentum. Then the complexity will be increased step by step by introducing masses, transverse momenta and spin.

The ð1 + 1ÞD yoyo model
In this model everything occurs in the ðt; zÞ hyperplane. One assumes that the sting has a uniform probability Pdzdt to break in the space-time area dzdt. From the quantum point of view, the "string fragility" P is taken into account by adding an imaginary part −iℏP=2 to the string tension κ ≃ 1 GeV=fermi [15][16][17]. The complex string tension κ C ¼ κ − iℏP=2 is analogous to the complex mass m − iℏγ=2 of an unstable particle. The decay products are small strings which oscillate like yoyos. Figure 2 shows the corresponding space-time history. The string world sheet (hatched domain) is bordered by quark world lines. The breaking points Q 2 ; Q 3 ; Á Á Á Q N , completed by the return points Q 1 and Q Nþ1 , form an a-causal chain, i.e., the 2-D vector Q r Q s is spacelike. The one-point breaking density in the ðt; zÞ plane is The exponential is the probability that no string breaking occurred in the past light cone of Q. Breaking at Q r creates a quark pair fq rqr g. q r andq rþ1 meet at point H r to form the yoyo h r , which represents a stable hadron or a resonance, depending on its mass. Its 2-momentum p r is given by where the "check" symbol means interchanging the time and z components, i.e.,v ¼ ðv z ; v 0 Þ for any vector v ¼ ðv 0 ; v z Þ. In principle, a yoyo h r can further break, simulating the decay of a resonance. However one limits the model to the production of these "primary" hadrons. It already reproduces salient properties of the hadronic final states: two back-to-back jets, Feynman scaling, charge retention effect and rapidity plateau. A weak point of the model is its fully continuous mass spectrum starting at m ¼ 0.
Recursive treatment.-In the "P z ¼ −∞" frame, or using X − ¼ t − z as time variable, the hadron emission points H 1 ; Á Á Á H N are ordered in time according to their ranks in Eq. (2). This allows to treat the sting decay as a recursive quark fragmentation, identifying q r of the q rqr pair with q r of Eq. (2). In the string picture the 4-momentum conservation k r ¼ p r þ k rþ1 applies to the canonical momentum of the quark, given by At any point q of the rth quark trajectory we have where k mec r , given byǩ mec r ¼ κqQ r , is the mechanical momentum of the quark (see Fig. 3). It is on-mass-shell, but varying along the quark line. k string r is the momentum flow, from right to left, of the string through Oq and plays the role of a linear 2-potential. It is given byǩ string r ¼ κOq. Any recursive model can be uniquely defined by the double density of consecutive quarks in momentum space. In the ð1 þ 1ÞD yoyo model it is with b L ¼ P=ð2κ 2 Þ and k ≡ k can . Note that k þ > 0 and k − < 0 for all quarks, except for k − The exponential factor is the probability that no string breaking occurs in the past light cone of H, as shown in the right picture of Fig. 4. From Eq. (9) we obtain the single quark density equivalent to Eq. (5), and the splitting function Multiperipheral feature.-The ð1 þ 1ÞD yoyo model can be cast in the form of a multiperipheral model with quark exchanges, 2 as pictured in Fig. 1. In Fig. 4 a vertex is represented both in the string and in the multiperipheral picture. The squared vertex function is and the squared propagator A connection between a QCD multiperipheral model and the string fragmentation model is also discussed in Ref. [18].
2. Introducing transverse momenta and actual hadron masses: The Lund symmetric splitting function (LSSF) Classically, string breaking can only create massless quarks with zero transverse momenta. To overcome this limitation one assumes that quarks pairs are created by a tunneling mechanism analogue to the Schwinger mechanism of e þ e − pair creation in a strong electric field. The quark and the antiquark of the rth pair have then opposite transverse momenta, k rT and −k rT , which are absorbed by hadrons h r and h r−1 . Figure 2 is then an approximate classical picture, a point Q r representing only the middle of a tunneling path as shown more in detail in Fig. 5. Giving to the hadrons their actual masses and using the principle of "left-right" symmetry or "quark line reversal" (hereafter referred to as "LR symmetry"), the authors of [3] came to the Lund symmetric splitting function (LSSF), The inputs of the model are the parameters a q and the function w q 0 ;h;q ðk 02 T ; p 2 T ; k 2 T Þ, which depends on the quark flavors q and q 0 , the hadron species h and the transverse momenta. In the most general model a q ≡ a q ðk 2 T Þ. The function u q ðk 2 T Þ is a normalizing factor. For LR symmetry w must be symmetrical under fq; k 2 T g⇌fq 0 ; k 02 T g together with h⇌h.
Like the ð1 þ 1ÞD yoyo model, the Lund symmetric model can be cast in a multiperipheral form represented in Fig. 1. Equations (9)-(13) for the quark densities, the vertex, the propagators and the splitting function become where F is given by Eq. (14). w q 0 ;h;q is normalized such that any timelike curve passing by O in Fig. 2 is crossed by one and only one QQ 0 segment. The power-law factors lead to a multi-Regge behavior for large rapidity gaps [15]. In a semiclassical approach, the quantum actions of the quarks produce such factors, with a q ðk 2 8]. Note that the vertex and the propagator are not invariant under the full Lorentz group, but under the subgroup generated by (a) the rotations aboutẑ, (b) the Lorentz boosts alongẑ, (c) the reflection about any plane containingẑ. Indeed, the string axis defines a privileged direction of space.
We take w of the form C q 0 ;h;q is proportional to the ðq 0 qÞ wave function in flavor space. It acts upon the hadron species distribution. f T ðk 2 T Þ is a fast decreasing function of k 2 T (e.g., a Gaussian).ǧðϵ 2 h Þ acts upon the correlation 3 between k T and k 0 T , since (14). In PYTHIA, such a correlation is absent, due to the particular where a q ðk 2 T Þ ¼ a is taken flavor-and k 2 T -independent.
FIG. 5. Tunneling process of a q rqr pair. 3 Such correlations are present in the standard multiperipheral model, where hk 0 B. The classical string + 3 P 0 mechanism One assumes that the string breaking, in which a q rqr pair is created (r is the rank of the splitting), occurs via a tunnel effect and that, at the end of this process, q r andq r are on the string axis (the z axis), with transverse momenta k rT and −k rT respectively, zero longitudinal momenta, and separated by the vector The string between q r andq r has been "eaten" by the pair. The modulus of d r is fixed by energy conservation and its orientation is that of the initial color flux, i.e., from q A tō q B . The quark pair has a relative orbital momentum One furthermore assumes that the q rqr pair is in the 3 P 0 state (which possesses the quantum numbers of the vacuum). In such a state the spins are parallel and opposite to L r : It follows from (23), (24), and (25) that the polarizations of q r andq r are correlated to their transverse momenta: Besides (25), which correlates s q r and sq r , there is a correlations between s q r and sq rþ1 coming from the internal wave function of the meson h r . In particular, if h r is a pseudoscalar meson (π, K, η or η 0 ), as required by the 1 S 0 internal wave function. Figure 6 depicts the spin and k T correlations in the recursive decay of the sting when only pseudoscalar mesons are emitted and assuming that q A is polarized along þŷ (as represented by an counterclockwise arrow). According to (25) and (27), q 2 andq 2 are both polarized along −ŷ (clockwise arrow) and their relative orbital momentum L 2 is along þŷ (counterclockwise arrow). Thenq 2 and q 2 move respectively in the þx and −x directions, in accordance with (26). The transverse momentum −k 2T ofq 2 , which is toward þx, is absorbed by h 1 , resulting in a Collins effect with hp 1;x i > 0, more generally hp 1T × S A;T i ·ẑ > 0.

C. Quantum treatment of the quark spin
We encode the quark spin degree of freedom with Pauli spinors and, using the multiperipheral approach, transform the vertex function V and the propagator D of Eqs. (16)-(17) into 2 × 2 matrices acting on quark spin. w, u and the quark density U of Eqs. (16)- (19) become density matrices (Hermitian and semipositive definite) in spin space. Full Lorentz invariance would require the use of Dirac spinors, but Pauli spinors are sufficient to satisfy the invariance under the above mentioned subgroup. Note that it does not take into account the whole spin information (2 q-bits) carried by an off-mass-shell Dirac particle.

General formalism
We first consider a general mutiperipheral model, not necessarily combined with the string model. The amplitude for reaction (1) is To save place, the gothic letters gather several variables: for a quark q ¼ fq; kg, where q is the flavor; for a hadron h ¼ fh; p; s h g, where h is the hadron species and js h i belongs to an adopted spin basis (e.g., helicity basis). Thus, DðqÞ ≡ Dðq; kÞ and Vðq 0 ; h; qÞ ≡ V q 0 ;h;s h ;q ðk 0 ; kÞ. jSi is the Pauli spinor of polarization S ¼ ðS T ; S L Þ, with T and L referring to the transverse and longitudinal polarizations of the quark respectively. jS B i is related to the polarization SB of the antiquarkq B by which is the analog of the Dirac spinor vðk; SÞ ¼ −γ 5 uðk; −SÞ of an antiparticle. The functions DðqÞ and Vðq 0 ; h; qÞ may be chosen as input of the model. However they can be "renormalized" by the transformation DðqÞ → Λ L ðqÞDðqÞΛ R ðqÞ; without changing the physical result, so different inputs lead to the same model. Here we use the "renormalized input" method of [9], where the 4 × 4 matrix is the density operator of pairs of consecutive quarks. i, j label the spin states for q and i 0 , j 0 label the spin states for q 0 . In Eq. (31) we have introduced the partially transposed matrix Equation (31) generalizes Eq. (15). In fact N is a density matrix in spin space but a classical density in momentum space, since we still treat momenta classically (see footnote 1). The single-quark density operator, generalizing Eq. (17), is Invariance under reflection about the ðx; zÞ or ðy; zÞ plane requires whereñðkÞ ≡ẑ × k T =jk T j. U is Hermitian and semipositive definite: U 0 ðqÞ and U 1 ðqÞ are real with U 0 ðqÞ ≥ jU 1 ðqÞj. We assume the strict inequality so that a solution of Eq. (33) for DðqÞ exists. DðqÞ is not uniquely determined by Eqs. (33) and (34) but, using the "renormalization" (30), we can take the positive definite solution DðqÞ ¼ U −1=2 without loss of generality. Introducing the splitting matrix the polarized splitting function to be used in Eq. (4) becomes where ρðqÞ ¼ ð1 þ S q · σÞ=2 is the spin density matrix of quark q, normalized to unit trace. F obeys the normalization condition From the practical point of view Eq. (36) is used to draw the species h, the spin state s h and the momentum p of hadron h.
The spin density matrix of the left-over quark q 0 is Thus Eqs. (36) and (38) are the basis for the recursive generation of a polarized quark jet.

Combination with the string model
For the vertex V we take which generalizes Eq. (16). gðq 0 ; h; qÞ ¼ g q 0 ;h;s h ;q ðk 0 T ; k T Þ is a 2 × 2 matrix acting on quark spin. It also contains the flavor and k T dependence of V. 4 Factorizing the quark propagator as then Eq. (33) becomes where we have defined The term w of Eq. (19) has been replaced by the matrix product g † g.
Particular choices of a q and gðq 0 ; h; qÞ.-We take a q ¼ a q 0 ¼ a ¼ constant and gðq 0 ; h; qÞ of the form This procedure allows to generalize Eq. (21) when the quark spin is taken into account. The term μ q þ σ z σ · k T is the 2 × 2 analogue of the numerator m q þ γ · k of the Feynman propagator. The function f T provides the cutoff in k T . Inspired by the Schwinger mechanism, we take the μ q is a complex parameter having the dimension of a mass. With Imðμ q Þ > 0, the factor μ q þ σ z σ · k T , introduced in [7], reproduces the classical string þ 3 P 0 mechanism. Γ is a 2 × 2 matrix which depends on k T and k 0 T at most as a polynomial. In Eq. (54) and in the next sections we will restrict ourselves to pseudo scalar mesons and, to zero order in k T and k 0 T , we will take which is the analogue of γ 5 . In Ref. [12] the slightly T =2Þ was made. It gives practically the same result. Vector and axial mesons can also be introduced as shown in Ref. [7].
Particular choices ofǧðϵ 2 h Þ.-As in the spinless case, this function acts upon the spin-independent ðk T ; k 0 T Þ correlation which adds to the one mediated by the quark spin. Let us give four examples: Choice (a) favors k T · k 0 T > 0. In choice (b) this correlation is reinforced for c < 0 and weakened for c > 0. Choice (c) suppresses it, like in PYTHIA or in the simplified model of [7]. In choice (d) the factor e cb L ϵ 2 h =2 restores it.
D. The polarized q → h + q 0 splitting function according to the 3 P 0 mechanism and a unique complex mass parameter for all quark flavors, i.e., μ q ≡ μ ¼ ðReμ; ImμÞ. With our choice ofǧ two successive quark transverse momenta k T and k 0 T are correlated. Gathering Eqs. (36), (44), and (50) we obtain the polarized splitting function to be used in (4) The free parameters a and b L play the same role as a and b in the PYTHIA event generator. They govern the suppressions of F at large and small Z respectively. The parameter b T is linked to the spread of the quark transverse momenta produced at the string cutting points.ρ int ðqÞ is an intermediate density matrix which we have not normalized. The corresponding polarization vector is Working out the trace operations in Eq. (51) with Γ h ¼ σ z , the splitting function is explicitly given by The tilde denotes the "dual" of a transverse vector, for instancẽ p T ¼ẑ × p T . The vector S int is the polarization vector of the intermediate spin matrixρ int ðqÞ given in Eq. (53).
Finally, using Eqs. (38) and (50), the spin density matrix ρðq 0 Þ of the quark q 0 is calculated as

IV. MONTE CARLO IMPLEMENTATION
In this section we describe the simulation code handling the fragmentation of a polarized quark into pseudoscalar mesons (π, K, η 0 and η 0 ). It is a stand alone program not yet interfaced with existing event generators. Presently, the flavor and the spin density matrix ρðq A Þ of the fragmenting quark q A are chosen at the beginning of the simulation. The initial quark energy can either be fixed or chosen event by event reading the values from an external file. The output consists in a file with the relevant information on all the hadrons generated in the fragmentation, later on analyzed to obtain the azimuthal angle distributions and the analyzing powers.
A. The Monte Carlo program structure A preliminary task, before starting the generation of the events, is to calculateû q ðk T Þ from Eqs. (47) and (44), then u −1=2 q ðk T Þ and tabulate its values. The initial kinematics for lepton-proton DIS is defined event by event according to the hard subprocess l þ q 0 → l 0 þ q A . We consider the center of mass frame of the system composed by the virtual photon and the target proton. We orient the z axis along the virtual photon momentum and consider first the case where q A has no primordial transverse momentum k Tprim . This reference frame coincides also with the center of mass frame of the final hadronic system whose light cone momenta P AE are defined by where For the e þ e − annihilations, W coincides with the center of mass energy and it is the same for all the events of a simulation.
In our reference frame the quark q A travels along the forward light cone and one can identify k þ A ≡ P þ (implying k þ B ¼ 0). We only consider the fragmentation of this initial quark and neglect the jet initiated byq B , which in the DIS case is a diquark and travels along the backward light cone with momentum k − B ≡ P − (implying k − A ¼ 0). We simulate the splitting process q → h þ q 0 recursively, starting with q ¼ q A , following the steps: (1) generate a new q 0q0 pair (2) form h ¼ qq 0 and identify the type (π; K; η 0 or η 0 ) of the pseudoscalar meson (3) generate Z according to the p T -integrated splitting function and calculate p þ ¼ Zk þ (4) generate p T according to the splitting function at the generated Z and calculate k 0 test the exit condition: if it is not satisfied continue to step 7, otherwise the current hadron is removed and the decay chain ends (7) calculate the hadron four-momentum and store it (8) calculate the spin density matrix of quark q 0 using Eq. (55) with Γ h ¼ σ z and come back to step 1. We iterate steps 1-8 until the exit condition, described below, is satisfied. More details on the different steps are given in the following. (a) Quark flavor and hadron type generation (steps 1 and 2).
The meson identification at step 2 uses the isospin wave function and also suppresses the η 0 meson production with respect to π 0 to account for their mass difference. We have chosen Nðη 0 Þ=Nðπ 0 Þ ≃ 0.57 as suggested in Ref. [1]. (b) Exit condition (step 6). After the r th splitting, the 4-momentum of the remaining string is P rem;rþ1 ¼ kB þ k rþ1 . Then The remaining squared energy to be used in the generation of the next hadrons is W 2 rþ1 ¼ P þ remðrþ1Þ P − remðrþ1Þ − P 2 T;rþ1 . If the last hadron is generated with a very small value of Z, P − remðrþ1Þ can become negative. In this case the last hadron is rejected and a new one is tried. If W 2 rþ1 falls below a given mass M 2 R the chain terminates and the last hadron generated is erased (exit condition at step 6). We take M R ¼ 1.5 GeV=c 2 in order to leave enough energy for the production of one baryon, which is not simulated. The observables investigated here are not sensitive to this value. (c) Recursive splitting (steps 3 and 4). The energymomentum sharing in the splitting qðk T ; k þ Þ → hðp T ; Zk þ Þ þ q 0 ðk 0 T ; ð1 − ZÞk þ Þ is performed using the splitting function given in Eq. (54). The differential probability to produce a hadron with light cone momentum fraction Z is given by the integral of the splitting function over p T and writes The choice of the values of the parameters a, b L , b T and μ entering Eq. (60) will be discussed in the next section. The terms which affect mostly the distribution of Z are (i) the parameter a which suppresses large values of Z by a power law (ii) the exponential exp ð−b L m 2 h =Z − b T ξðZÞk 2 T Þ which depends on m h and on k 2 T . This shifts Z toward larger values when k 2 T is large. To be more precise, the first rank hadron h 1 is generated in the splitting A Þ, hence Z 1 is drawn according to Eq. (60) with vanishing k T and only the mass m h enters the exponential. At the next step the hadron of rank two h 2 is generated in the splitting , therefore Z 2 is shifted towards larger values with respect to Z 1 because now in Eq. (60) enters a not vanishing k 2 2T . At this point the splitting of q 3 is similar to that of q 2 and no differences are expected for hadrons of rank two or higher.
After the generation of Z we draw p T according to the differential probability where k T is fixed from the previous hadron generation. For ImðμÞ > 0, the S int ·k 0 T term of Eq. (61) pushes k 0 T in the direction ofẑ × S int and contributes to the Collins effect through p T ¼ k T − k 0 T . This reproduces the classical 3 P 0 mechanism. An other consequence of the string þ 3 P 0 mechanism is the spin-mediated correlation between k T and k 0 T , which is negative for pseudoscalar meson emission, as can be seen in Fig. 6 between k 2T and k 3T . On the other hand, the exponential factor in Eq. (61) naively forces the relation p T ∼ ð1 − ξðZÞÞk T and since 0 < ξðZÞ < 1 for every value of Z, the effect is that the transverse momenta of two successive quarks tend to be aligned, as already mentioned. Hence in our model there are two effects at work, which are opposite for the case of pseudoscalar mesons: the 3 P 0 mechanism and the dynamical correlations between the transverse momenta of quarks in the decay chain. Since the p T distribution strongly depends on k T due to the exponential factor in Eq. (61), we expect differences between the p 2 T distributions of the first and second rank hadrons and no further change for higher rank hadrons.

B. Values of the free parameters
The values of the four free parameters a, b L , b T , and jμj 2 have been tuned comparing the simulation results for unpolarized quark fragmentations with the measured p 2 T distributions of charged hadrons produced in SIDIS off unpolarized deuteron [20] and with a set of unpolarized p T -integrated fragmentation functions from global fits [21], in order to find a qualitative agreement. The slope of the p 2 T spectrum is sensitive to b L and b T while its detailed shape for p 2 T → 0 is sensitive to jμj 2 . The slopes are not affected by a, which changes the fragmentation functions at large hadron fractional energy z h and has been fixed comparing with the p T -integrated FFs. In this work we have used a ¼ 0.
The parameter ImðμÞ has been fixed comparing the simulated and the measured Collins asymmetries extracted from e þ e − annihilation data, as explained in the next section. Knowing already jμj 2 , ReðμÞ is given up to a sign, which cannot be fixed by transverse spin asymmetries alone. In all the simulations we use ReðμÞ ¼ 0.42 GeV=c 2 , ImðμÞ ¼ 0.76 GeV=c 2 .
Note that all the results but those in Sec. V D have been obtained with a vanishing primordial transverse momentum.

C. Kinematical distributions
Let us first look at the kinematical (spin-independent) distributions with the chosen values of the parameters. Figure 7(a) shows the distributions of the longitudinal splitting variable Z for the first four rank hadrons generated in the fragmentation chain of a u quark. As can be clearly seen, the distribution of the first rank hadron is shifted towards smaller values of Z with respect to the distributions of higher rank hadrons. This is due to the k T dependence of the splitting function discussed above. Also, the distributions of higher rank hadrons are similar, as expected.
It is interesting to compare Fig. 7(a) with Fig. 7(b) showing the z h distributions for the first four rank hadrons. The shapes of the z h distributions are similar for rank 1 and 2 and change sensibly with the hadrons rank r because of the relation z h r ≃ Z r ð1 − Z r−1 Þ…ð1 − Z 2 Þð1 − Z 1 Þ. By definition the Z and z h distributions for the rank one hadron coincide.
The k 02 T distributions for the different splittings are very much the same and, since the initial quark has vanishing k T , the k 02 T distribution of the leftover quark in the first splitting coincides with the p 2 T distribution of the first rank hadron. As a consequence the p 2 T distribution of the first rank hadron is softer than the distributions of higher rank hadrons, as shown in Fig. 8(a). The slope of the h 1 distribution is almost twice the slope of h 2 distribution.
1T i. This difference between the p 2 T of rank 1 and rank 2 hadrons is a common feature of all recursive fragmentation models. A trace of it is the inequality hp 2 T ðh þ Þi < hp 2 T ðh − Þi for medium values of z h in u jets, shown in Fig. 8(b). However the data suggest the opposite. This discrepancy could be reduced with the introduction of resonances, like vector mesons, and their decay. The decrease of hp 2 T i at small z h is due to the factor expð−b L p 2 T =ZÞ in the splitting function [Eq. (51)].

V. RESULTS ON THE TRANSVERSE SPIN ASYMMETRIES
In order to study the transverse spin effects measured in SIDIS off transversely polarized protons and in e þ e − annihilation, fragmentation events have been generated for initial quarks fully polarized along a fixedŷ axis orthogonal to the string axis. Only the results of the dominant u quark fragmentation are shown in the following. We have checked that the results for pion production in d quark fragmentation are related to those of the u quark fragmentation by isospin symmetry.
For the SIDIS case we have used a sample of real COMPASS events. The x B and Q 2 of these events serve to fix the initial kinematics of our simulation event-by-event.
For the study of the asymmetries in the azimuthal distributions of the hadrons produced in e þ e − annihilation a second sample of events has been generated with a fixed c.m. energy W ¼ 10 GeV corresponding to the BELLE energy.
In this section we present the results on the single hadron and the dihadron transverse spin asymmetries and discuss the kinematical dependences of the corresponding analyzing power. The effect of the primordial transverse momentum is also described. The Monte Carlo (MC) results are compared only with the COMPASS and BELLE data, which are in quite good agreement with the corresponding results from HERMES [22] and Jefferson Lab experiments [23] and from BABAR [24] and BESIII [25] experiments respectively.

A. Single hadron transverse spin asymmetries
The well known Collins effect [5] is the left-right asymmetry in the distribution of the hadrons produced in the fragmentation of a transversely polarized quark with respect to the plane defined by the spin and the momentum of the quark. The azimuthal distribution of the hadrons of the jet is given by where S AT is the quark transverse polarization. The angle ϕ C ¼ ϕ h − ϕ S A is the Collins angle, where ϕ h and ϕ S A are the azimuthal angles of the hadron momentum and of the quark spin. The analyzing power a q A ↑→hþX is the ratio between the spin dependent part of the FF (the Collins FF) and the unpolarized quark FF. Experimentally, the Collins effect has been observed in SIDIS, where the Collins FF couples with the transversity PDF, and in e þ e − , where the measured azimuthal asymmetry can be written in terms of products of two Collins FFs. The sin ϕ C dependence of Eq. (62) is effectively obtained in our simulation. The fact that our model is formulated at the amplitude level guaranties that no other azimuthal modulation is present and that positivity constraints are satisfied. Using simulated events, the analyzing power a q A ↑→hþX is calculated as 2hsin ϕ C i and in general it is function of both z h and p T . The Collins asymmetry measured from e þ e − annihilation data has been used to fix the value of the free parameter ImðμÞ ¼ 0.76 GeV=c 2 . More specifically we have compared the mean value of the Collins analyzing power for positive pions in transversely polarized u jets from simulations with the mean value 0.258 AE 0.006 obtained in Ref. [26] from BELLE data. Figure 9 shows the Collins analyzing power a u↑→hþX as function of z h for charged pions and kaons (left panel) and as function of p T for charged pions (right panel). The main feature is that the analyzing power has opposite sign and almost equal magnitude for oppositely charged mesons, as qualitatively expected from the 3 P 0 model. The mean values for hadrons with p T > 0.1 GeV=c and z h > 0.2 are given in Table I.
Also, the analyzing power vanishes for small z h and is almost linear in the range 0.2 < z h < 0.8. A linear dependence on z h is also suggested by the BELLE data [26] when the analyzing power for the favored fragmentation is assumed to be opposite to that for unfavored fragmentation.
The sign and the monotonic dependence of the analyzing power on z h can be understood by writing a u↑→hþX as the sum of different rank hadron contributions weighted by the number of hadrons of that rank. The analyzing power can be written as where the variable "t" can be either z h or p T . N h r is the number of hadrons of type h and of rank r and a u↑→h r þX is the analyzing power associated with rank r, both calculated at the same value t. The analyzing power for the different rank hadrons is shown in Fig. 10. It has opposite sign for even and odd ranks, as suggested by Fig. 6, and decreases with the rank. Such decrease is due to the depolarization of the recurrent quark which, with the current choice of parameters, turns out to be a weak effect. Indeed in each splitting roughly 10% of the recurrent quark transverse polarization is lost. The main cause of decay of the analyzing power at small z h is the mixture of contributions from even and odd ranks. The fact that the z h dependence is roughly linear, not another power law, is a priori accidental. Concerning the sign of the analyzing power, for an initial u quark, a fast positive pion can be produced at first rank or at rank r > 1 following r − 1 π 0 's or η's. On the contrary a negative pion can never be produced at first rank because of its charge. Furthermore the contribution of larger ranks is  smaller because N h r ðrÞ decreases with rank due to the finite W. Thus the signs of the π þ and π − analyzing powers are fixed by the contributions of the first and second ranks, respectively. The same considerations can be made for charged kaons.
From the left panel of Fig. 9, we notice also that the slope for negative mesons, which are unfavored in u chains, is slightly larger than the slope for positive ones. This effect is easily explained by the fact that the absolute value of the analyzing power for a rank two hadron is somewhat larger than the analyzing power for a rank one, as can be seen from Fig. 10. Finally we can see that the slope for π − and K − are similar, as expected because both start to be produced from rank two.
Concerning the analyzing power as function of p T , shown in the right panel of Fig. 9, there are clearly different behaviors for positive and negative mesons. An interesting feature is the change of sign of the analyzing power for positive pions at p T ≃ 0.9 GeV=c. The rank analysis at this value of p T shows that the number of π þ of rank 1 and 3 is roughly the same as the number of π þ of rank 2 and 4. Moreover positive pions with large p T are more likely produced as rank two, following a rank one π 0 or η, than as rank one. This effect combines with the alternate sign of the analyzing power for even and odd rank hadrons to give a u↑→π þ þX ðp T ¼ 0.9 GeV=cÞ ≃ 0. The number of higher rank pions decreases quickly and they give only a small contribution to the asymmetry.
Similar trends are observed in the Collins asymmetry for charged pions produced in SIDIS off transversely polarized protons as measured by COMPASS [27]. The comparison with Monte Carlo results is shown in Fig. 11 as function of z h (left plot) and as function of p T (right plot). The Monte Carlo values in both panels are those of Fig. 9 multiplied by an overall scale factor λ 1 ¼ 0.055 AE 0.010 obtained from χ 2 minimization. In the u-dominance hypothesis for a proton target and neglecting the effect of the primordial transverse momentum, λ 1 is the ratio of the x B -integrated u-quark transversity and the x B -integrated unpolarized u quark density, multiplied by the depolarization factor D NN of lepton-quark scattering.
As apparent from the right panel of Fig. 11, the Monte Carlo describes qualitatively the p T dependence of the experimental points, which do not exclude a change of the π þ asymmetry sign for p T > 0.9 GeV=c.
The agreement between Monte Carlo and COMPASS asymmetries as function of z h is satisfactory for positive pions, whereas for negative pions it is poor for z h > 0.6.

B. Dihadron transverse spin asymmetries
The properties of the analyzing power a u↑→h 1 h 2 þX due to the Collins effect in the h 1 h 2 pair production in a u jet have also been studied. Such analyzing power has been found to be related to a u↑→h AE þX in a recent experimental work in SIDIS [28,29] and its magnitude can be obtained by e þ e − data [30].
In general the distribution of oppositely charged hadron pairs in the same jet, as function of the relevant variables used here, is given by

Comparison with BELLE data
In order to compare with the e þ e − data we have evaluated the quantity ϵðM inv Þ ≡ ha u↑→π þ π − þX i × a u↑→π þ π − þX ðM inv Þ, where ha u↑→π þ π − þX i is the analyzing power averaged over all the kinematical variables, including M inv . For this comparison, in the simulation the analyzing power a u↑→π þ π − þX has been estimated using Φ ¼ ϕ B where ϕ B is the azimuthal angle of the vector p 1T − p 2T . Figure 12 shows the results for ϵðM inv Þ from the simulation when z h1;2 > 0.1 with no cut in p T as circles whereas for those represented by squares we have required p T > 0. 3 GeV=c. The open triangles show the values of ϵ as measured by BELLE [30]. Both in the simulation and in the data the analyzing power shows a saturation for large values of the invariant mass while for small values it tends to zero. It has to be noted that the Monte Carlo data sample has a different invariant mass spectrum with respect to BELLE data. In particular in the BELLE data sample the statistics is larger in the region of the ρ meson while in the program there are no resonances and most of the statistics is at higher values of M inv . Still, both in BELLE and in simulation results, no structure can be seen.
We recall that, in order to cancel, or minimize, the effects due to the primordial transverse momenta, the dihadron asymmetry is normally written in terms of the azimuth of the relative transverse momentum For the BELLE results we are considering here, the vector characterizing the pair is where P T ¼ p 1T þ p 2T is the global transverse momentum of the pair. Defining as "pure" dihadron asymmetry the one defined with respect to the vector R T , the asymmetry extracted from the BELLE data is a combination of the pure dihadron asymmetry and of the global Collins effect of the pair.

Comparison with COMPASS data
In Fig. 13 we show the comparison between the Monte Carlo and the COMPASS dihadron asymmetry for h þ h − pairs measured in SIDIS off transversely polarized protons as function of z (left) and M inv (right). The dihadron asymmetry is extracted using Φ ¼ ϕ R where ϕ R is the azimuthal angle of the vector R T , thus it can be regarded as a pure dihadron asymmetry. Both in COMPASS data and in simulations the cuts z h > 0.1, x F > 0.1, R T > 0.07 GeV=c and jp i j > 3 GeV (i ¼ 1, 2) have been applied.
The left plot of Fig. 13 concerns the dependence on z. The Monte Carlo points are scaled by a factor λ 2 estimated by comparing with the COMPASS asymmetry as function of z. From a χ 2 minimization we obtain λ 2 ¼ 0.055 AE 0.008 in perfect agreement with the value of λ 1 obtained in the single hadron asymmetry case, as expected. The results from the Monte Carlo are in good agreement with the experimental data.
The right plot of Fig. 13 shows the dependence of the analyzing power on M inv . The same cuts as those for the dihadron asymmetry as function of z have been applied. After scaling by the same parameter λ 2 , the Monte Carlo points describe quite well the trend of the data.

C. Comparison between single hadron and dihadron transverse spin asymmetries
Following the work done in Ref. [28] we have studied the relationship between the Collins and the dihadron analyzing powers for hadron pairs in the same u quark jet, as function of the relative azimuthal angle Δϕ ¼ ϕ 1 − ϕ 2 . In that analysis using only the events with at least one h þ and one h − two kinds of asymmetries had been extracted: the "Collins-like" (CL) asymmetries A sin ϕ C CL1ð2Þ for positive and negative hadrons and the dihadron asymmetry for oppositely charged hadron pairs A sin ϕ 2h;S CL;2h . In each bin of Δϕ, the CL asymmetry is the Collins asymmetry of h þ (h − ) of the pair.
As in Ref. [28] we calculate a u↑→h þ h − þX using Φ ¼ ϕ 2h , where ϕ 2h is the azimuthal angle of the vectorp 1T −p 2T andp T ≡ p T =jp T j. Due to the relation the considered asymmetry is a combination of the "pure" dihadron asymmetry and of the global Collins asymmetry of the hadron pair. However, as already discussed in Ref. [28], the azimuthal angle ϕ R is strongly correlated with ϕ 2h , and the dihadron asymmetry measured from 2hsin ϕ 2h;S i with ϕ 2h;S ¼ ϕ 2h − ϕ S A , is essentially the same as the "pure" dihadron asymmetry, which could be verified with the code as well. The blue squares in Fig. 14(a) represent the dihadron analyzing power a u↑→π þ π − þX calculated in the Monte Carlo as function of Δϕ. The blue curve is the result of the fit with the function c ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ð1 − cos ΔϕÞ p as suggested in Ref. [28]. The plot in Fig. 14(b) shows the asymmetry A sin ϕ 2h;S CL;2h as measured in COMPASS. As can be seen, the agreement is good. Note that the A sin ϕ 2h;S CL;2h asymmetry is smaller than a u↑→π þ π − þX by a factor of 0.1 analogous to λ 2 but for the higher cut x B > 0.032 adopted in this experimental analysis.
The same considerations hold also for CL analyzing power A sin ϕ C CL1ð2Þ of h þ and h − shown in the top plot of Fig. 14(a) with red circles and black triangles respectively. The corresponding COMPASS data are shown in top plot of Fig. 14(b): again, the trend is very similar. The MC points are fitted with functions of the type δ 1ð2Þ þ c 1ð2Þ sin Δϕ, as suggested from Ref. [28], and the results are represented by the red and the black dashed lines. The red and the black dashed lines in Fig. 14(b) represent the fits to the experimental CL asymmetries as shown in Ref. [28], which are consistent with vanishing δ 1ð2Þ parameters.

D. Introducing the primordial transverse momentum
In the previous sections we did not consider the primordial transverse momentum of the initial quark. In this section we show the results when the initial quark q A does have a primordial transverse momentum. Figure 15 depicts the string direction in the DIS γ Ã -nucleon center of mass frame when the struck quark has primordial transverse momentum k T prim , inherited from the quark motion in the nucleon. k T prim , also written k T=q γ , is defined with respect to the γ Ã momentum q γ . The target remnant has the opposite −k T prim . The string is stretched between q A and the target remnant. Its axis is therefore rotated from the γ Ã -nucleon axis. The effect of a random k T prim is the broadening of the spectra of hadrons transverse momenta with respect to q γ . This should partly smear the single hadron asymmetry. The primordial momentum k T prim is generated according to the probability where hk 2 T prim i is a free parameter. The fragmentation of the initial transversely polarized quark q A is performed using the rotated string axis asẑ axis. Then we go in the laboratory frame with a boost along the γ Ã -nucleon axis.
In the small angle approximation, the rotation in the string center of mass frame is practically equivalent to make the following shift in p T (which is relative to the string axis) where p T=q γ is the hadron transverse momentum with respect to the γ Ã axis and x F ¼ ð2p z =WÞ c:m: is the Feynman scaling variable. Since x F ¼ z h − ϵ 2 h =ðz h W 2 Þ, Eq. (69) almost coincides for large x F with the often used equation p T=q γ ¼ z h k T prim þ p T . The shift is zero at x F ¼ 0 and opposite to k T prim in the backward hemisphere as can be guessed from Fig. 15.
From Eq. (69) follows at fixed x F hp 2 The effect of k T prim is clearly seen in Fig. 16 showing the hp 2 T=q γ i as function of z h for positive hadrons when the fragmenting quark has hk 2 T prim i ¼ 0.3 ðGeV=cÞ 2 . The large z h region, where z h ≃ x F , is more sensitive to the primordial transverse momentum and the effect decays for smaller values of z h . It turns out that the difference between hp 2 T i for positive and negative hadrons shown in Fig. 8(b) is somewhat reduced due to the x 2 F hk 2 T prim i term but still the negative hadrons are produced with larger transverse momenta.
In Fig. 17 we show the effect of the primordial transverse momentum on the Collins analyzing power as function of z h (left plot) and p T=q γ (right plot) for positive and negative pions. The analyzing power for hk 2 T prim i ¼ 0.3 ðGeV=cÞ 2 (full points) is compared to that for vanishing primordial transverse momentum (open points). The reduction of the analyzing power is visible at large z h (left plot) and at low p T=q γ (right plot). We note also that the change of sign of the analyzing power as function of p T=q γ is no more there. The same effects are also observed for charged kaons. Table II shows the mean values of the single hadron and dihadron analyzing powers for charged pions for different values of hk 2 T prim i. At variance with the Collins asymmetry for single hadrons, the asymmetry for pairs of oppositely charged hadrons is practically not affected by the noise introduced by k T prim .

VI. RESULTS ON THE JET HANDEDNESS
The present model can treat both longitudinal and transverse polarizations at the same time. In particular it can predict jet handedness [31,32] which for a particle pair h 1 h 2 and for a longitudinally polarized quark q A can be parametrized in the form The simplified model of Ref. [7] predicts such an effect with an analyzing power proportional to Imðμ 2 Þ. The same factor appears in the present model. We have made a simulation for π þ π − pairs in the jet of an initial longitudinally polarized u quark and calculated the analyzing power a ⃗ u→π þ π − þX JH as 2hsinðϕ 2 − ϕ 1 Þi. Figure 18 shows the dependences of a ⃗ u→π þ π − þX JH on the invariant mass M inv of the pion pair (left plot) and on the sum of their fractional energies z 1 þ z 2 (right plot). While we do not observe a strong dependence on M inv , the handedness analyzing power increases with z 1 þ z 2 . This is expected since at large z 1 þ z 2 both hadrons have nearly fixed ranks (rank 1 for π þ and rank 2 for π − ). Comparing Fig. 18 with Fig. 13, where we remind that the Monte Carlo analyzing power is scaled by the factor λ 2 , we find for the jet handedness an effect smaller than the dihadron asymmetry by one order of magnitude. Note that we obtain an opposite asymmetry by reversing the sign of ReðμÞ.
Up to now, attempts to observe jet handedness were not conclusive, see, e.g., Ref. [33]. Several reasons can explain this failure: (i) the sign of the asymmetry may vary too much with the charges, the rapidity ordering or the invariant mass of the h 1 − h 2 pair. (ii) the observable cosðϕ 2 − ϕ 1 Þ is very sensitive to a redefinition of the jet axis. It can be easily blurred by a too large experimental uncertainty on the orientation of the jet axis or by gluon radiation.
Like for the Collins effect, the blurring effect can be eliminated by involving one more particle. Indeed, for three particles h 1 , h 2 and h 3 of the jet, the pseudoscalar quantity where P ¼ p 1 þ p 2 þ p 3 , is independent of the jet axis and we may take hJi as helicity-sensitive estimator (the estimator hsignðJÞi was proposed in Ref. [32]). However it requires the clean measurement of three particle momenta and its amplitude depends on a 6 kinematical variables, e.g., z 1 , z 2 , z 3 , jp 1;⊥P j, jp 2;⊥P j, and jp 3;⊥P j.

VII. CONCLUSIONS AND PERSPECTIVES
We have developed a stand alone Monte Carlo code for the simulation of the fragmentation process of a polarized quark (u, d, or s). The theoretical framework is provided by the string fragmentation model where the quark-antiquark pairs in the string cutting points are produced according to the 3 P 0 mechanism. The quark spin is included through spin density matrices and propagated along the decay chain reproducing the string þ 3 P 0 mechanism.
With respect to the Lund symmetric model, this model requires an additional complex mass parameter whose imaginary part directly affects the single hadron Collins asymmetry. The three free parameters present in the string fragmentation framework and the absolute value of the complex mass have been tuned by comparison with unpolarized SIDIS data and with global fits of p T -integrated fragmentation functions.
The analyzing powers have been extracted from the simulated events both for the single hadron and for the hadron pairs. The results of the simulation show a Collins analyzing power of opposite sign for oppositely charged mesons. The dependence on the kinematical variables z h and p T has been investigated, finding a reasonable  Fig. 17 (left) for positive and negative pions with cuts z h > 0.1 and p T=q γ > 0.1 ðGeV=cÞ have been applied. We show also the mean values of the asymmetry for π þ π − pairs with the same cuts. agreement with experimental results. A clearly different from zero analyzing power for hadron pairs of opposite sign in the same jet is also obtained from the same simulated data. The Monte Carlo results are compared to BELLE and COMPASS dihadron asymmetries finding again a satisfactory agreement. Thus a unique mechanism generates both the Collins and the dihadron asymmetry. Furthermore with the same model we predict also a jet handedness effect in the fragmentation of a longitudinally polarized quark.
Such a model can be a guide to optimize the estimators of quark polarimetry. An interface of our Monte Carlo program with the PYTHIA event generator is foreseen. A further improvement of the model is the inclusion of resonances, in particular of vector mesons, and the generation of their hadronic decays.