{'fraction_non_alphanumeric': 0.044531452697124795, 'fraction_numerical': 0.013585866924546549, 'mean_word_length': 3.8896321070234112, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A microscopic description for polarization in particle scatterings

We propose a microscopic description for the polarization from the first principle through the spin-orbit coupling in particle collisions. The model is different from previous ones based on local equilibrium assumptions for the spin degree of freedom. It is based on scatterings of particles as wave packets, an effective method to deal with particle scatterings at specified impact parameters. The polarization is then the consequence of particle collisions in a non-equilibrium state of spins. The spin-vorticity coupling naturally emerges from the spin-orbit one encoded in polarized scattering amplitudes of collisional integrals when one assumes local equilibrium in momentum but not in spin.


I. INTRODUCTION
A very large orbital angular momentum (OAM) can be created in peripheral heavy ion collisions [9,10,24,45,46,58,60]. Such a huge OAM can be transferred to the hot and dense matter produced in collisions and make particles with spins polarized along the direction of OAM [24,37,45,60]. Recently the STAR collaboration has measured the global polarization Several theoretical models have been developed to study the global polarization. If the spin degree of freedom is thermalized, one can construct the statistic-hydro model by including the spin-vorticity coupling S µν ω µν into the thermal distribution function [5][6][7][20][21][22].
Here S µν is the spin tensor, ω µν = −(1/2)(∂ µ β ν − ∂ ν β µ ) is the thermal vorticity, the macroscopic analog of the local OAM, and β µ ≡ βu µ is the thermal velocity with β = 1/T being the inverse of the temperature and u µ being the fluid velocity. It turns out that the average spin or polarization is proportional to the thermal vorticity if the spin-vorticity coupling is weak.
In order to describe the STAR data on the global Λ/Λ polarization, the hydrodynamic or transport models have been used to calculate the vorticity fields in heavy ion collisions [4, 13-15, 38, 39, 54]. Then the polarization of Λ/Λ can be obtained from vorticity fields at the freezeout when the Λ/Λ hyperons are decoupled from the rest of the hot and dense matter [40,44,52,64].
All these models are based on the assumption that the spin degree of freedom has reached local equilbrium. But this assumption is not justified. The recent disagreement between some theoretical models and data on the longitudinal polarization indicates that the spins might not be in local equilibrium [8,48,63]. Although one model of the chiral kinetic theory can explain the sign of the data [53], it cannot reproduce the magnitude of the data. If the spins are not in local equilbrium, how is the polarization generated in particle collisions? This is also related to the role of the spin-orbit coupling which is regarded as the microscopic mechanism for the global polarization. In one particle scattering such as a 2-to-2 scattering at fixed impact parameter the effect of spin-orbit coupling in the polarized cross section is obvious [24,45], but how does the spin-vorticity coupling naturally emerge from the spin-orbit one? It is far from easy and obvious as it involves the treatment of particle scatterings at different space-ime points in a system of particles in randomly distributed momentum. To the best of our knowledge, this problem has not been seriously investigated due to such a difficulty. In this paper we will construct a microscopic model for the global polarization based on the spin-orbit coupling. We will show that the spin-vorticity coupling naturally emerges from scatterings of particles at different space-time points incorporating polarized scattering amplitudes with the spin-orbit coupling. This provides a microscopic mechanism for the global polarization from the first principle through particle collisions in non-equilibrium.
The paper is organized as follows. In Section II we will introduce scatterings of two wave packets for spin-0 particles. The wave packet method is necessary to describe particle scatterings at different space-time points. In Section III we will study collisions of spin-0 particles as wave packets which take place at different space-time in a multi-particle system. In Section IV we will derive the polarization rate for spin-1/2 particles from particle collisions. As an example, we will apply in Section V the formalism to derive the quark polarization rate in a quark-gluon plasma in local equilibrium in momentum. In Section VI we will discuss the numerical method to calculate the quark polarization rate, a challenging task to deal with collision integrals in very high dimensions. We will present the numerical results in Section VII. Finally we will give a summary of the work and an outlook for future studies.
Throughout the paper we use natural units = c = k B = 1. The convention for the metric tensor is g µν = diag(+1, −1, −1, −1). We also use the notation a µ b µ ≡ a · b for the scalar product of two four-vectors a µ , b µ and a · b for the corresponding scalar product of two spatial vectors a, b. The direction of a three-vector a is denoted asâ. Sometimes we denote the components of a three-vector by indices (1,2,3) or (x, y, z).

II. SCATTERINGS OF WAVE PACKETS FOR SPIN-0 PARTICLES
In this section we will consider the scattering process A + B → 1 + 2 · · · + n, where the incident particles A and B in the remote past are localized in some region and can be described by wave packets. The details of this section can be found in the textbook by Peskin and Schroeder [49]. The purpose of this section is to give an idea of how the wave packets displaced by an impact parameter are treated in the scattering process, and to provide the basis for the discussion in the next section. We work in the frame in which the central momenta of two wave packets are collinear or in the same direction which we denote as the longitudinal direction. We assume that the wave packet B is displaced by an impact parameter vector b in the transverse direction, so the in state can be written as Here we see that the incident particles are treated as two wave packets |φ A and |φ B defined in Appendix A. The definition of the single particle states |k A and |k B can also be found in Appendix A. As we have mentioned that the amplitudes φ i (k i ) center at p i = (0, 0, p iz ) for i = A, B. We assume that the out state is a pure momentum state |p 1 p 2 · · · p n out in the far future. This is physically reasonable as long as the detectors of final-state particles mainly measure momentum or they do not resolve positions at the level of de Broglie wavelengths.
Taking into account the normalization factors for the in-state and out-state, the scattering probability is given by where the normalization of single particle states and wave packets is given in Appendix A.
Since P(AB → 12 · · · n) depends on the impact parameter b, we can rewrite it as P(b).
This probability gives the differential cross section at the impact parameter b, The total cross section is then an integral over the impact parameter . We can also Here we have used the normalization condition for the wave amplitude (A9). We note that the above formula is derived in the frame in which incident particles are collinear in momemtum. We can boost the frame to the center-of-mass frame of the incident particles and the cross section is invariant.
If the number densities of A and B in coordinate space are n A and n B respectively, the collision rate, i.e. the number of scatterings per unit time and unit volume is given by where we have rewritten the rate in a Lorentz invariant way by making use of the fact that In this section we will derive the collision rate in a system of spin-0 particles of multispecies. We will generalize the result of the previous section by treating the incident particles as wave packets. The emphasis is put on the collision of two particles at two different spacetime points.
We will frequently use two frames in this and the next section: the lab frame and the center-of-mass system (CMS) of one specific collision. In the lab frame, the movement of one species of particles follows their phase space distribution f (x, p). There are many collisions taking place in the system. Figure 1 shows one collision of two incident particles as shown in the right panel of Fig. 1, see Appendix C for more details of such a Lorentz transformation. Hereafter we denote the quantities in the CMS by the index 'c'. There is an inherent problem in the collision of incident particles located at different space-time points: the collision time is not well defined. If we assume that the collision takes place at the same time in the lab frame, i.e. t A = t B , after being boosted to the CMS, the time will be mis- The reverse statement is also true: Such an ambiguity in the collision time cannot be avoided but can be constrained by the requirement that the difference ∆t c = t c,A − t c,B should not be large, otherwise the incident particles are irrelevant or the collision is un-causal in the CMS. In the calculation of this paper, we will put a simple constraint ∆t c = 0. In the right panel of Fig. 1, we also see that the impact parameter b is given by the distance of x c,A and x c,B in the transverse direction which is perpendicular to p c,A or p c,B . In the longitudinal direction or the direction of p c,A or p c,B , two space points are also different in general, i.e.p c,A · x c,A =p c,A · x c,B . In the calculation we also require that the distance between two space points in the longitudinal direction, ∆x c,L =p c,A · (x c,A − x c,B ), should not be large, otherwise the incident particles as wave packets lose coherence and cannot interact in the CMS. In the calculation, we will also put a simple constraint ∆x c,L = 0. The CMS constraint ∆t c = 0 and ∆x c,L = 0 is equivalent to the condition ∆t = v bst · ∆x and (v A − v B ) · ∆x = 0 in the lab frame, see Appendix C for the derivation. Since we will work in the CMS of incident particles in each collision, for notational simplicity, we will suppress the index 'c' (standing for the CMS) of all quantities in the rest part of this section. So all quantities are implied in the CMS if not explicitly stated here.
We know that the momentum integral of the distribution function gives the number density in the coordinate space. Similar to Eq. (8), the collision rate in corresponding momentum and space-time intervals can be written as where in the CMS, f A and f B are the phase space distributions for the incident particle A and B respectively, and ∆σ denotes the infinitesimal element of the cross section given by Here we have assumed that the scattering takes place at the same time and the same longitudinal position in the CMS, so we put two delta functions to implement these constraints.
The constant C AB is to make ∆σ have the right dimension of the cross section and will be defined later. In Eq. (10) K is given by In Eq. (11) G i (i = 1, 2) denote distribution factors depending on particle types in the final state, we have G i = 1 for the Boltzmann particles and G i = 1 ± f i (p i ) for bosons (upper sign) and fermions (lower sign). Note that f i (p i ) can be in any other form in non-equilibrium cases. In (11) we have taken the following form for Here we take the Gaussian form for the wave packet amplitude φ i (k i − p i ) as in (A10), where α i denote the width parameters of the wave packet A or B. For simplicity we will set equal width for two incident particles (even for different species), α A = α B = α.
We can also make the approximation of narrow wave packets, so we have |k i | ≈ |k i | ≈ |p i | for i = A, B and then E A,k E A,k ≈ E A and E B,k E B,k ≈ E B , and the energy factors in (11) drop out. By taking the integral over x A and x B and then the integral over on-shell momenta p A , p B , p 1 and p 2 , we obtain the scattering or collision rate per unit volume, We can rewrite the integral over x A and x B in Eq. (14) as two delta functions in Eq. (11). In Eq. (16) we have integrated over y 0 = ∆t = t A − t B and to remove two detla functions, then we are left with the integral over the transverse part y µ T = (0, b) with b being in the transverse direction. Because we work in the CMS in which all kinematic variables depend on the incident momenta in the lab frame, the impact parameter b in the CMS depends on (x A , x B ) as well as (p A , p B ) in the lab frame through a boost velocity. Now we define the constant C AB in (9,14) as C AB ≡´d 4 X = t X Ω int so that the final results have the right dimension. Here t X and Ω int are the local time and space volume for the interaction respectively. Note that C −1 AB´d 4 X (· · · ) plays the role of the average over X or (· · · ) X . If we take the limit t X Ω int → 0, we obtain the local rate per unit volume from Eq. (14), where N AB→12 is the number of scatterings. We emphasize again that all quantities in Eq.
(17) are defined in the CMS of two incident particles (we have suppressed the index 'c').
In this section we will generalize the previous section for spin-0 particles to spin-1/2 ones. Our purpose is to derive the polarization rate from collisions in a system of particles of multi-species. We assume that particle distributions in phase space are independent of spin states, so the spin dependence comes only from scatterings of particles carrying the spin degree of freedom.
As a simple example to illustrate the idea of the polarization arising from collisions, we consider a fluid with the three-vector fluid velocity in the z direction v z that depends on x, which we denote as v z (x). We assume dv z (x)/dx > 0. In the comoving frame of any fluid cell in the range [x − ∆x/2, x + ∆x/2] where ∆x is a small distance, the fluid velocity In a thermal system with collective motion, there are many scatterings whose reaction planes point to almost random directions, but in average the direction of the reaction plane points to that of the local rotation or vorticity. To calculate the polarization in a thermal system with collective motion, we have to take a convolution of distribution functions and polarized scattering amplitudes similar to (17).
In this section we will distinguish quantities in the CMS and lab frame, i.e. we will resume the subscript 'c' for all CMS quantities, while quantities in the lab frame do not have the subscript 'c'. respectively. The quantization direction of the spin state is chosen to be along the direction of the reaction plane in the CMS of the scattering. The polarization rate per unit volume for particle 2 in the final state is given by where P AB→12 denotes the polarization vector and n c =b c ×p c,A is the direction of the reaction plane in the CMS of the scattering which is also the quantization direction of the spin. In the second to the last line of Eq. (18), the summation of 2s 2 M(· · · , s 2 )M * (· · · , s 2 ) over s 2 = ±1/2 gives the polarized amplitude squared for particle 2 in the final state, and the factor 2 arises from the normalization convention for the polarization that makes it in (18) is one of our main results.

PLASMA OF LOCAL EQUILIBRIUM IN MOMENTUM
In this section we will calculate the quark/antiquark polarization rate from all 2-to-2 parton (quark or gluon) collisions in a quark-gluon plasma (QGP) of local equilibrium in momentum but not in spin. We assume that the QGP is a multi-component fluid with the same fluid velocity u(x) as a function of space-time for all partons. The partons in a fluid cell follow a thermal distribution in momentum in its comoving frame with the local temperature four-momentum of the parton and β(x) ≡ 1/T (x).
We consider the scattering, A + B → 1 + 2, where A and B denote two incident partons in the wave packet form localized at x A and x B respectively, and '1' and '2' denote two outgoing partons in momentum states. In order to calculate the polarization rate from the collision of two wave packets displaced by an impact parameter by Eq. (18), we must work in the CMS of the incident partons for each collision. Note that many collisions take place in the system at different space-time, the CMS of each collision depends on the momenta of incident partons which vary from collision to collision. In one collision, the phase space distributions for incident partons (denoted as i = A, B) can be written in the form where x, p are the space-time and momentum in the lab frame respectively, while x c , p c are their corresponding values in the CMS of A and B in this collision which depend on p A and p B in the heat bath (lab frame) through the boost velocity, and u µ c (x c ) denotes the fluid velocity in the CMS as a function of the space-time in the CMS.

A. Polarization rate
We now apply Eq. (18) to 2-to-2 parton scatterings. For simplicity we assume that the phase space distributions of incident partons follow the Boltzmann distribution, i.e. (18). Also we assume that y c,T is small compared with X c so that we can make an expansion in y c,T for the distributions, the details are given in Appendix B. The relevant contribution in the linear or first order in y c,T involves the term y µ c,T [∂(βu c,ρ )/∂X µ c ]p ρ c,A which can be rewritten as where µρ shows up in the y c,T expansion, which can be converted to the spin-vorticity coupling through polarized parton scattering amplitudes encoding the spin-orbit coupling effect, as we will show shortly.
The second term in Eq. (20) invloves the symmetric part of the thermal velocity derivatives in space-time, which is assumed to vanish in thermal equilibrium for the spin, known as the Killing condition [5][6][7]22]. In this paper, however, we do not assume the thermal equilibrium for the spin degree of freedom, so we keep this symmetric term in the calculation.
Keeping the first order term in the y c,T expansion and neglecting the zeroth order term which is irrelevant, Eq. (18) can be simplified as where we have used where the factor 2 arises from the normalization convention for the polarization. Note that in the above formula there is a sum over color degrees of freedom of all incident and outgoing partons. We may write ∆I AB→12 M n c as Here we have used and Ja(k L c,A ) are Jacobians for the longitudinal momenta k L c,A and k L c,A and are given by and k L c,A (j 1 ) and k L c,A (j 2 ) with j 1 , j 2 = 1, 2 are two roots of the energy conservation equation (24) and (25) Latin indices label spatial components in the the CMS. The derivation of (24) is given in Appendix D.
In a system of gluons and quarks with multi-flavors, there are many 2-to-2 parton scatterings with at least one quark in the final state. The quark polarization rate for a specific flavor reads where d 4 P AB→1q (X)/dX 4 is given by Eq. (24), and 2-to-2 parton scatterings are listed in Table I. The antiquark polarization rate can be similarly obtained.

B. Polarized amplitudes for quarks/antiquarks in 2-to-2 parton scatterings
In this subsection we will derive the polarized amplitudes for quarks in 2-to-2 parton scatterings. The Feynman diagrams of all 2-to-2 parton scatterings at the tree level with at least one quark in the final state are shown in Table I. For anti-quark polarization, we can make particle-antiparticle transformation in all processes listed in Table I, for example, gg → q aqa , etc.. In this subsection, we discuss polarized amplitudes for quarks, those for antiquarks can be easily obtained.
In order to obtain the quark polarization, we sum over the spin states of all partons in the scattering except one quark in the final state. For simplicity of the calculation, we assume that the quark masses are equal for all flavors and the external gluon is massless.
We introduce a small mass in the gluon propagator in the t-channel to regulate the possible divergence.
In this subsection, all variables are defined in the CMS, for notational simplicity we will suppress the subscript 'c', for example, p A actually means p cA . Here a and b denote the quark flavor, s i = ±1/2 (i = A, B, 1, 2) denote the spin states, k i (i = A, B, 1, 2) denote the momenta, q, q 1 , q 2 , q 3 denote the momenta in propagators. The processes for antiquark polarization can be obtained by making a particle-antiparticle transformation.
We take the quark-quark scattering q a q b → q a q b with a = b (different flavor) as an example to demonstrate how to derive the polarized scattering amplitude which depends on the spin state of the quark in the final state. The Feynman diagram of this process is shown in Table   I. The spin-momentum configurations are shown in the diagram. We can then write down the corresponding amplitudes following the Feynman rule where g s is the strong coupling constant, i, j, k, l = 1, 2, 3 denote the fundamental colors of quarks, c, d = 1, · · · , 8 denote the adjoint colors of gluons, t c and t d are generators of SU (N c ) in fundamental representation satisfying [t a , t b ] = if abc t c , q = k A − p 1 , and q = k A − p 1 . We obtain the product I 1 I * 2 as In Eq. (29) we have used the notation p · γ ≡ p ρ γ ρ , a sum over all spins except s 2 and over all colors of quarks and gluons have been taken, and C qaq b →qaq b is the color factor for this process given in Table II The can be expressed as a function of (k A ,k A ,k B ,k B ) and (η kA , η kA , η kB , η kB ).
color factors Color factors in scattering processes qaqa→qaqa , C gg→qaqa , C (4) gqa→gqa The polarized amplitudes for quarks in all 2-to-2 parton scatterings listed in Table I  behind such huge number of terms unless we make an appropriate approximation.

C. Evaluation of polarized amplitudes for quarks/antiquarks
The evaluation of contracted traces of quark polarized amplitudes are very complicated.
This has been done with the help of FeynCalc [47,50]. There are about 10 4 terms in the expansion of contracted traces for 2-to-2 parton scatterings.
In this subsection, all variables are defined in the CMS, for notational simplicity we will suppress the subscript 'c' if not explicitly specified, for example, p A actually means p cA .
In order to show the physics in the midst of the huge number of terms, we have to make an appropriate approximation. As we know that the incident particles are treated as wave packets in order to describe scatterings displaced by impact parameters. A realistic approximation is that the wave packets are assumed to be narrow, i.e. the width is much smaller than the center momenta of the wave packet in Eq. (13). In the extreme case that the width of the wave packet is zero, we recover the normal scattering of plane waves. Since the positions of incident particles can be anywhere in plane waves, in average the relative OAM of two incident particles is zero, leading to the vanishing polarization of final state particles. This fact can be verified by settinĝ in the trace part in Eq. (30), then we have The above result is of the zeroth order, now we turn to the first order in the deviation from momenta in (31). We expand and (η kA , η kA , η kB , η kB ) about their central values (η pA , η pA , η pA , η pA ) to the first order in the differences,k where the first order quantities are denoted with ∆ (for example, ∆ A , ∆η kA ). We also The delta functions in Eq. (21) lead to So ∆ 1 in (33) can be determined by and p 0 determined by Note that once p 0 and ∆ 1 are given, E 0 , ∆ 1 , ∆ 2 satisfy So we have a freedom to choose the value of E 0 . Then we use (32) and (33) in the contracted trace part in Eq. (30) and expand it to the first order in ∆-quantities. Still, the final result has many terms but all terms of ∆ 1 , ∆ 2 and ∆ 1 cancel out.
In order to further simplify the contracted trace part in Eq. (30), we use the property that the first order contributions do not have terms of ∆ 1 , ∆ 2 , ∆ 1 by setting which leads to k A + k B = k A + k B = 0 and then Using (38) and (39) in the contracted trace part in Eq. (30) for q a q b → q a q b , we obtain a where we denote the contracted trace part for q a q b → q a q b as Tr 2 qaq b →qaq b , c A ≡ cosh(η kA /2), c A ≡ cosh(η kA /2), s A ≡ sinh(η kA /2), and s A ≡ sinh(η kA /2). We see in (40) that there are four typical terms proportional to (n × p 1 ) ·k A , (n × p 1 ) ·k A , (n ×k A ) ·k A and (p 1 ×k A ) ·k A , in which the first three terms are from the spin-orbit coupling and the last one corresponds to the non-coplanar part of p 1 ,k A andk A . We will show in the next section that (40) is a good approximation for the contracted trace part to the exact result.
It can be proved that ∆I AB→12 M for all 2-to-2 parton scatterings in Table I have the same structure as in (40) for q a q b → q a q b under the approximation in (38,39).
Note that ∆I AB→12 M depends linearly on the direction of the scattering plane n =b ×p A , we can write the contracted trace part in the form ofb · I, as is done in Eq. (23). We take the term (n × p 1 ) ·k A in (40) as an example, which can be rewritten as Therefore I contains the term inside the square brackets on the right-hand side of Eq. (41).
Another example is the term proportional to (p 1 ×k A )·k A , we see that all terms have factors of the form n · V (V =k A ,k A , p 1 ) inside the square brackets, these terms can be rewritten as n · V =b · (p A × V), so I contains the termp A × V.

VI. NUMERICAL METHOD TO CALCULATE QUARK/ANTIQUARK POLARIZATION RATE
In this section we will calculate the polarization rate for quarks in a QGP from Eq. (24).
Here we assume a local equilibrium in particle momentum but not in spin. We will consider two cases: the approximation as in (38,39)  Although the 2-to-2 processes for anti-qaurk polarization are different from those for quarks, it can be shown that the polarization rate for anti-quarks is the same as that for quarks, because all 2-to-2 scatterings for anti-quark polarization can be obtained from those in Table I by making a particle-antiparticle transformation. In the following we discuss only the quark polarization. The same discussion can also be applied to the antiquark polarization.
The local polarization rate in Eq. (24) for quarks involves a 16-dimensional integration, which is a major challenge in the numerical calculation. In the Monte Carlo integration, the number of sample points grows exponentially with the dimension, so even a very rough calculation in high dimensions would need huge number of sample points.
To overcome this difficulty, we split the integration into two parts: a 10-dimension (10D) integration over (p c,1 , p c,2 , k T c,A , k T c,A ) and a 6-dimension (6D) integration over (p A , p B ). We carry out the 10D integration and store the result as a function of p c,A (and p c,B = −p c,A ).
Then we carry out the 6D integration using the pre-calculated 10D integral.
The 10D integral, the last five lines of Eq. (24), depends on p c,A and p c,B = −p c,A which appear in the wave packet function φ A and φ B respectively. So we denote the 10D integral as Θ jk (p c,A ), from Eq. (27) the polarization rate per unit volume for one quark flavor can be rewritten as where the second equality defines W ρν and the sum of A, B, 1 is over all 2-to-2 processes in Table I.

A. The 10D integration
The 10D integral Θ jk (p Then we sample p c,1 and p c,2 . In order to increase the efficiency of the sampling, we should determine the range of p c,1 and p c,2 . We can first determine the ranges of lengths |p c,1 | and |p c,2 | by a numerical search. Then we determine the ranges of directionsp c,1 and p c,2 . For a givenp c,1 , which can be randomly chosen, we find that the largest value of θ ≡ arccos(−p c,1 ·p c,2 ) betweenp c,2 and −p c,1 occurs when Hence we obtain the range of θ as c,A ) for later use. It takes a couple of days to finish the calculation. We find that when |p c,A | > 2.5 GeV, the 10D integral is almost zero. This is due to the fact that if α |p c,A |, the incident wave packets can be almost regarded as plane waves which give vanishing polarization.

B. The 6D integration
Now we carry out the remaining 6D integration over p A and p B in (42). As we have mentioned in Section V that we assume partons with p µ A = (E A , p A ) and In the preceding subsection, we calculated the 10D integral Θ jk (p c,A is in the z direction. We have to transform the tensor Θ jk (p c,A is rotated to the real direction of p c,A determined by Eq. (C1). The rotation matrix R ij is defined by p c,A,i = R ij p (z) c,A,j , with which we define the transformation for the tensor Our numerical results show that the tensor W ρν has the form where we see that ρ and ν should be spatial indices or W 0ν = W ρ0 = 0. The form of (45) will be verified in the numerical results in Section VII. Then from (42) we obtain the polarization rate per unit volume for one quark flavor where , and for spatial indices we have the 3D form with u being the spatial part of the four-velocity u ρ .

VII. NUMERICAL RESULTS
In this section we will present our numerical results. The approximation in (38,39) is inspired by the first order contribution in the narrow wave packet approximation. In order to see how effective the approximation is, we compare in Fig. 2 the results of the 10D c,A ) for the scattering processes q(q) + q → q(q) + q and g + q → g + q in two cases: with and without the approximation. Here the process q(q) + q → q(q) + q stands for a sum over 5 different processes in Table I. Note that we do not show the results for g +g → q +q for which all elements of Θ jk (p  c,A ) for q(q) + q → q(q) + q and g + q → g + q in two cases: (1) with the approximation in (38,39) and (2) (38,39). The unit of Θ jk (p We then work out the rest 6D integral and obtain W ρν in Eq. (45). In the 6D integration we have to determine the maximum value of |p A | and |p B | or the integration range of |p A | and |p B |. In Fig. 3, as an example, we show the dependence of W 31 y on |p A | max = |p B | max for q(q) + q → q(q) + q, where we choose b 0 = 2.2 fm, z = 0 fm and T = 0.3 GeV. We see in the figure that the value of W 31 y is very stable when |p A | max = |p B | max > 8T . The numerical results for W ρν show the structure of (45). We can write W ρν in an explicit matrix form As an example, we show in Fig. 4 the results for all components of W 31 as functions of the cutoff b 0 for the quark polarization. We see in the figure that W 31 x and W 31 z are two or three orders of magnitude smaller than the positive values of W 31 y , which gives the polarization in the y direction. As we can see in the figure that W 31 y increases with the cutoff b 0 . The reason for such a rising behavior is due to the Taylor expansion of f A (x c,A , p c,A )f B (x c,B , p c,B ) to the linear order in y c,T = (0, b c ) as in App. B. There should exist an upper limit for b 0 above which the coherence of the incident wave packets is broken and the results are not physical. Such an upper limit can be set to be the order of the hydrodynamical length scale ∼ 1/∂ µ X u ν and should be larger than the interaction length scale 1/m D . It can be proved that W 31 for the anti-quark polarization is the same as that for the quark one. The numerical results show that the magnitude of all element W ρν are equal so we denote it as W . We have constructed a microscopic model for the global polarization from particle scatterings in a many body system. The core of the idea is the scattering of particles as wave packets so that the orbital angular momentum is present in scatterings and can be converted to spin polarization. As an illustrative example, we have calculated the quark/antiquark polarization in a QGP. The quarks and gluons are assumed to obey the Boltzmann distribution which simplifies the heavy numerical calculation. There is no essential difficulty to treat quarks and gluons as fermions and bosons respectively.
To simplify the calculation, we also assume that the quark distributions are the same for all flavors and spin states. As a consequence, the inverse processes that one polarized quark is scattered by a parton to two final state partons as wave packets are absent. So the relaxation of polarization cannot be described without inverse processes and polarized distributions. We will extend our model by including the inverse processes in the future.

IX. SUMMARY AND CONCLUSIONS
The global polarization in heavy ion collisions arises from scattering processes of partons or hadrons with spin-orbit couplings. However it is hard to implement this microscopic picture consistently to describe particle scatterings at specified impact parameters in a thermal medium with a shear flow. On the other hand the statistic-hydro model or Wigner function method are widely used to calculate the global polarization in heavy ion collisions. These models are based on the assumption that the spin degrees of freedom have reached a local equilibrium. So there should be a spin-vorticity coupling term in the distribution function to give the global polarization proportional to the vorticity when it is small. However it is unknown if particle spins are really in a local equilibrium. In this paper we aim to construct a microscopic model for the global polarization from particle collisions without the assumption of local equilibrium for spins. The polarization effect is incorporated into particle scatterings at specified impact parameters with spin-orbit couplings encoded. The spin-vorticity coupling naturally emerges from particle collisions if we assume a local equilibrium in particle momenta instead of particle spins. This provides a microscopic mechanism for the global polarization from the first principle through particle collisions in non-equilibrium.
As an illustrative example, we have calculated the quark polarization rate per unit volume from all 2-to-2 parton (quark or gluon) scatterings in a locally thermalized quark-gluon plasma in momentum. Although the processes for anti-quark polarization are different from those for quarks, it can be shown that the polarization rate for anti-quarks is the same as that for quarks because they are connected by the charge conjugate transformation. This is consistent with the fact that the rotation does not distinguish particles and antiparticles.
The spin-orbit coupling is hidden in the polarized scattering amplitude at specified impact parameters. The polarization rate involves an integral of 16 dimensions, which is far beyond the capability of the current numerical algorithm. We have developed a new Monte-Carlo integration algorithm ZMCintegral on multi-GPUs to make such a heavy task feasible. We have shown that the polarization rate per unit volume is proportional to the vorticity as the result of particle scatterings, a non-equilibrium senario for the global polarization. So we can see in this example how the spin-vorticity coupling emerges naturally from particle scatterings. In this appendix, we will give definitions and conventions for the single particle state in coordinate and momentum space and those for the wave packet.

Single particle state in coordinate and momentum space
For simplicity we first consider the single particle state of spin-0 particles, then we generalize it to spin-1/2 particles.
A position eigenstate is denoted as |x and satisfies following orthogonality and completeness conditions The normalization of the state |x is then where Ω is the space volume.
A momentum eigenstate is denoted as |p and satisfies following orthogonality and completeness conditions where E p = |p| 2 + m 2 is the energy of the particle. Note that p |p is Lorentz invariant.
The normalization of |p is then From Eq. (A1) and (A3) we can define the inner product x|p as With the above relation we can check where we have inserted the completeness relation in (A3). We can express |x in terms of |p and vice versa, (A7)

Single particle state as a wavepacket
In the real world a particle is always localized in some finite region, so its state can be represented by a wavepacket |φ which is a superposition of plane wave states, and φ(k) is the amplitude and can be normalized to unity, The energy dimension of |φ is 0. A typical form for φ(p) satisfying Eq. (A9) is the Gaussian which is centered at p 0 . The wavepacket function in coordinate space is where we have used Eq. (A5).
If we displace the particle state by b in coordinate space, the new wavepacket function is given by where the new wavepacket state is For spin-1/2 particles, the single particle state |k, λ has a spin index λ which is the spin along a quantization direction. The orthogonality and completeness conditions in (A3) now become k , λ |k, The wavepacket has the form and satisfies the normalization condition φ, λ|φ, λ = 1 similar to Eq. (A9).

Appendix C: Lorentz transformation
In the lab frame two colliding particles have on-shell momenta p A = (E A , p A ) and p B = (E B , p B ). The Lorentz transformation for the energy-momentum from the lab frame to the CMS of two colliding particles is where i = A, B, v bst is the boost velocity or the velocity of CMS in the lab frame and is given by and is the Lorentz contraction facror corresponding to v bst . Equation (C1) defines the Lorentz transformation matrix Λ µ ν . The reverse transformation to (C1) from the CMS to the lab frame can be obtained by flipping the sign ofv bst , The above defines the Lorentz transformation matrix The difference of two space-time points in the CMS are expressed in lab frame variables, where ∆t = t A − t B and ∆x = x A − x B . We then express the impact parameter as Let us look at the CMS constraint δ(∆t c )δ(∆x c,L ) in Eq. (10) (we have recovered the subscript 'c'). The condition ∆t c = 0 leads to while the conditionp c,A · ∆x c = 0 leads to where we have used We carry out the integration over the impact parameter and show how to remove the delta functions by integration in Eq. (21).
Substitute Eq. 23 into Eq. 21, we have the integration of b c in the following form Then we carry out the derivatives on a j , a k and a l , where we have used w 0 = ab 0 with b 0 being the upper limit or cutoff of b c , J i (i = 0, 1, 2) are Bessel functions, and a 3 a 2 a k δ lj + a 2 a l δ jk + a 2 a j δ lk − 3a l a j a k .
Note that the overall minus sign of Eq. (D3) cancels the one in Eq. (21).
We carry out the integration to remove the delta functions. First we integrate over k c,B and k c,B to remove six delta functions in three momenta, the result is to make following replacement in the integrand We are left with two delta functions for energy conservation which can be removed by the integration over k L c,A and k L c,A , where 'L' means the longitudinal direction along p c,A . To this purpose, we express the energies in terms of longitudinal and transverse momenta So two delta functions for energy conservation become where the Jacobians of two delta functions are given by and k L c,A (i = 1, 2) and k L c,A (i = 1, 2) are two roots of the energy conservation equation E c,A + E c,B − E c,1 − E c,2 = 0 and E c,A + E c,B − E c,1 − E c,2 = 0, respectively. The explicit forms of k L c,A (i = 1, 2) and k L c,A (i = 1, 2) are where C 1 and C 2 are given by with H being defined by Appendix E: Some formula for Dirac spinors The Hamiltonian for a Dirac fermion with the mass m is given by where γ µ = (γ 0 , γ) are Dirac gamma-matrices, α ≡ γ 0 γ, and σ = (σ 1 , σ 2 , σ 3 ) are Pauli matrices. The energy eigenstate can be found from the equation where E p = p 2 + m 2 , the sign ± in the right-hand side corresponds to positive/negative energy state, χ and φ are Pauli spinors which form a Dirac spinor (χ, φ). We can express χ in terms of φ and vice versa, where η = ±1 correspond to the positive and negative energy state respectively. So the positive energy solution becomes where s = ±1 is the spin orientation of the Pauli spinor and n = (sin θ cos φ, sin θ sin φ, cos θ) is the spin quantization direction. The spin eigenstates along n are given by which satisfy The negative energy solution can be put into the form The Dirac spinor for anti-particles can be defined by or defined in terms of the positive energy solution, v(s, p) = iγ 2 u * (s, p) The two Dirac spinors in (E8) and (E9) are actually the same up to a sign. Now we rewrite the Dirac spinor of a moving particle in the way of a Lorentz transformation of the one in the particle's rest frame. The Lorentz transformation matrix for the Dirac spinor is given by wherep ≡ p/|p| is the momentum direction, η p is the rapidity satisfying E p = m cosh(η p ), |p| = m sinh(η p ), v p = tanh(η p ), E p + m = 2m cosh 2 1 2 η p , E p − m = 2m sinh 2 1 2 η p . So u(s, p) can be expressed by a Lorentz boost of u(s, 0) for the particle at rest, In the same way we can rewrite v(s, p) as where we have usedū(s, q) =ū(s, 0)Λ −1 1/2 (−q),v(s, q) =v(s, 0)Λ −1 1/2 (−q), The spin projector is defined by where n σ is the Lorentz boost of the polarization vector (0, n) in the particle's rest frame satisfying n · p = 0 and n 2 = −1. In the particle's rest frame, we have Π rest (s, n) = 1 2 (1 + sn · Σ) Π(s, n)u(s, p) = u(s, p), As an example, we can explicitly verify the first one as Π(s, n)u(s, p) = 1 2 Λ 1/2 (−p)u(s, 0) where we have used Λ −1 1/2 (−p)γ σ Λ 1/2 (−p) = Λ ν σ (−p)γ ν and Λ ν σ (−p) = Λ ν σ (p). Using the spin projector, we have the following relation Π(s 0 , n) s u(s, p)ū(s, p) = Π(s 0 , n) (p · γ + m)| p µ =(Ep,p) = u(s 0 , p)ū(s 0 , p), where p · γ ≡ p µ γ µ .
Appendix F: Polarized amplitudes for quarks in 2-to-2 parton scatterings In this appendix, we give polarized amplitudes for quarks in all 2-to-2 parton scatterings listed in Table I. We assume the same quark mass m for all flavors and that the external gluon is massless. We introduce a mass into internal gluons or gluon propagators in the t and u channel to regulate the possible divergence.
All kinematic variables are defined in the CMS in this appendix, for notational simplicity we will suppress the subscript 'c' for all variables, for example, p A actually means p cA . The values of color factors, denoted as C AB→CD for the process A + B → C + D, are given in Table II. 1. q a q b → q a q b with a = b Following the Feynman diagram in Table I, we obtain the difference in the squared amplitude between the spin state s 2 = 1/2 and s 2 = −1/2 for q b in the final state, where q = k A − p 1 and q = k A − p 1 are momenta in the propagators.