{'fraction_non_alphanumeric': 0.029269657206794832, 'fraction_numerical': 0.009714169463940596, 'mean_word_length': 4.559875583203732, 'pattern_counts': {'":': 0, '<': 7, '<?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': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A topological realization of spin polarization through vortex formation in collisions of Bose-Einstein condensates

The global spin polarization of hadrons in heavy ion collisions has been measured in STAR (the Solenoidal Tracker At Relativistic heavy ion collider) experiments, which opens up a new window in the study of the hottest, least viscous and most vortical fluid that has ever been produced in the laboratory. We present a different approach to spin polarization from conventional ones: a topological realization of spin polarization through quantum vortex formation in collisions of Bose-Einstein condensates (BEC). This approach is based on the observation that the vortex is a topological excitation in a superfluid in presence of local orbital angular momentum and is an analogue of spin degrees of freedom. The formation processes of vortices and vortex-antivortex pairs are investigated by solving the Gross-Pitaevskii Equation with a large-scale parallel algorithm on Graphics Processing Unit (GPU) to very high precision. In a rotating environment, the primary vortex with winding number one is stable against perturbation, which has a minimal energy and fixed orbital angular momentum (OAM), but the vortices with larger winding numbers are unstable and will decay into primary vortices through a redistribution of the energy and vorticity. The injection of OAM can also be realized in non-central collisions of self-interacting condensates, part of the OAM of the initial state will induce the formation of vortices through concentration of energy and vorticity density around topological defects. Different from a hydrodynamical description, the interference of the wave function plays an important role in the transport of energy and vorticity, reflecting the quantum nature of the vortex formation process. The study of the vortex formation may shed light on the nature of particle spin and spin-orbit couplings in strong interaction matter produced in heavy-ion collisions.

The global spin polarization of hadrons in heavy ion collisions has been measured in STAR (the Solenoidal Tracker At Relativistic heavy ion collider) experiments, which opens up a new window in the study of the hottest, least viscous and most vortical fluid that has ever been produced in the laboratory. We present a different approach to spin polarization from conventional ones: a topological realization of spin polarization through quantum vortex formation in collisions of Bose-Einstein condensates (BEC). This approach is based on the observation that the vortex is a topological excitation in a superfluid in presence of local orbital angular momentum and is an analogue of spin degrees of freedom. The formation processes of vortices and vortex-antivortex pairs are investigated by solving the Gross-Pitaevskii Equation with a large-scale parallel algorithm on Graphics Processing Unit (GPU) to very high precision. In a rotating environment, the primary vortex with winding number one is stable against perturbation, which has a minimal energy and fixed orbital angular momentum (OAM), but the vortices with larger winding numbers are unstable and will decay into primary vortices through a redistribution of the energy and vorticity.
The injection of OAM can also be realized in non-central collisions of self-interacting condensates, part of the OAM of the initial state will induce the formation of vortices through concentration of energy and vorticity density around topological defects. Different from a hydrodynamical description, the interference of the wave function plays an important role in the transport of energy and vorticity, reflecting the quantum nature of the vortex formation process. The study of the vortex formation may shed light on the nature of particle spin and spin-orbit couplings in strong interaction matter produced in heavy-ion collisions.

I. INTRODUCTION
The rotation of a free body induced by a change in its magnetic moment is called the Einstein-de Haas effect which was discovered in 1915 1 . It is a consequence of the angular momentum conservation. A twin effect that the rotation of an uncharged material leads to a spontaneous magnetization was discovered by Samuel Barnett in the same year and named the Barnett effect 2 . These two effects reveal an intrinsic connection between the rotation and spin polarization in materials.
In ultra-relativistic heavy-ion collisions, similar phenomena also exist: huge orbital angular momenta (OAM) can be generated in peripheral collisions with respect to the reaction plane. A part of OAM is transferred to the hot and dense medium in the form of the global polarization of hadrons 3-8 (see, e.g. [9][10][11][12][13][14][15] , for recent reviews). The microscopic mechanism of converting OAM to hadrons' spin polarization is the spin-orbit coupling in particle scatterings 3,8,[16][17][18][19] , resulting in the spin-vorticity coupling of the fluid 20-27 when an ensemble average over randomized collisions is taken. The global polarization hyperons such as Λ (including Λ) can be measured through their weak decays into hadrons that break the parity 3 .
A non-vanishing global polarization of Λ hyperons in Au+Au collisions at √ s N N = 7.  GeV was measured in STAR experiments 28,29 .
It is hard to measure the polarization of mesons from their strong decays due to parity conservation. However, the 00-component of the spin density matrix ρ 00 for the vector meson can be measured through the angular distribution of its decay daughter in its rest frame. In the unpolarized case, ρ 00 should be 1/3, which implies that three spin states in the spin quantization direction (say, the z-direction) with S z = −1, 0, 1 are equally populated.
When ρ 00 is not 1/3, three spin states are not equally populated and the spin is called to be aligned in the spin quantization direction. Recently a large positive deviation of ρ 00 from 1/3 for φ mesons has been measured by the STAR collaboration 30 , which may be explained by fluctuating strong force fields 31,32 .
In this paper, we will present another approach to spin polarization: a topological real-ization of spin polarization through quantum vortex formation in collisions of Bose-Einstein condensates (BEC). Quantum vortex in BEC originates from topological defects and is an analogue of spin degrees of freedom. Therefore the formation process of quantum vorticies in collisions of BEC may shed light on the nature of the particle spin as well as spin-orbit and spin-vorticity couplings in the strong-interaction matter.
Although the collisions of BEC in this work will be described by the Gross-Pitaevskii equation for non-relativistic superfluids, they can be generalized to relativistic superfluids in the future. There are a few examples for relativistic superfluids. One example is the relativistic BEC far from equilibrium which can be formed dynamically in an over-populated initial state of gluons at the early stage of heavy ion collisions 39,40 . Another example is the pion condensation in the cold and dense nuclear matter 41 . At zero temperature, the isospin asymmetric quark matter will undergo a phase transition into the pion condensate when the isospin chemical potential is larger than the pion mass [42][43][44][45][46] . These relativistic superfluids in heavy ion collisions or neutron stars can be described by relativistic hydrodynamics 47 . The infrared dynamics of relativistic scalar fields can be mapped to non-relativistic complex scalar fields, so that the universal scaling of infrared modes and vortex dynamics in superfluids can be studied in both relativistic and non-relativistic regimes 48 . In the present work, we will focus on the vortex dynamics of non-relativistic systems. The relativistic version of the same study will be addressed in the future.
The outline of the paper is as follows. In Sec. II, we will give an introduction to vortex formation in BEC. In Sec. III, we will derive the Gross-Pitaevskii equation and the conservation laws for the particle number, energy-momentum, OAM and vorticity. In Sec. IV, we will discuss the stationary properties of a single vortex in the uniform BEC. Comparing to the vortex with winding number one (called primary vortex), the vortex with a larger winding number is unstable and subject to decay into primary vortices. In Set. V, with a high precision numerical simulation performed on Graphics Processing Unit (GPU), we will demonstrate the stability of the primary vortex as well as the decay process of vortices with larger winding numbers. In Set. VI, we will perform a simulation of the non-central collision process of two BECs. The finite OAM of the initial state can generate an array of vortices, where the energy and vorticity are concentrated around topological defects. With larger collision velocities, vortex-antivortex pairs can also be generated without initial seeds of topological defects. The conclusions and discussions will be made in Set. VII.
Convention. The coordinate three-vector r is denoted in boldface and its three components are denoted as r i with i = 1, 2, 3 or x, y, z corresponding to three basis vectors (e 1 , e 2 , e 3 ) or (e x , e y , e z ). The space derivative operator as a three-vector is denoted as As topological defects, bound vortex-antivortex pairs or unbound vortices can be gener-ated in a two-dimensional system in thermal equilibrium through the Berezinskii-Kosterlitz-Thouless (BKT) topological phase transition 58,59 . The construction pattern of vortices is employed to define a topological order. This long range order is based on the overall property of the system rather than the behavior of a two-point correlation function, the latter is used to define a traditional phase transition while the former is associated with a topological phase transition.

II. A BRIEF
In a system crossing the second-order phase transition, vortices can be spontaneously created by the growing and merging of domains with independent order parameters, known as Kibble-Zurek mechanism (KZM) [60][61][62][63] . The KZM is originally introduced in cosmology, but it can be applied to a wide variety of physical systems. Related to such a mechanism, the formation of Kibble-Zurek solitons in a Bose-Einstein condensate has also been extensively studied 56,57,64-67 .
Vortex-antivortex pairs can also be created in an over-populated boson system through strong turbulence 48 . The dynamics of vortex plays an important role in understanding how the off-equilibrium system evolves toward a quasi-stationary nonthermal fixed point [68][69][70] .
The creation, scattering and decay of vortex-antivortex pairs are strongly correlated with the scaling properties of the self-similar evolution [71][72][73] . The strongly anomalous scaling can be understood through three-body collisions of vortex-antivortex pairs with unbound vortices as well as collisions between vortex-antivortex pairs 48,73 .
Vortices can be formed in cold dark matter halos if they have sufficient number density, angular momentum and self interaction 74 . When de Broglie wavelengths of bosonic dark matter particles are larger than average inter-particle distances, the interference of matter wave inside a halo can also leads to vortex string or vortex rings with non-trivial windings 75 .
The formation, deformation and reconnection of vortex rings in BECs or in viscous fluids are intensively studied [76][77][78] . It would be interesting to study the signature of these vortex structures in the dark matter superfluid, for instance, through gravitational lensing 79,80 . If low-density vortices in galaxies can be confirmed, the superfluidity can potentially distinguish the scenarios of some dark matter models.
During the vortex nucleation in a rotating BEC, the symmetry breaking and quantum nature are studied by comparing the symmetry properties of the true many-body state and its mean-field approximation 81 . At the critical frequency, the exact many-body ground state exhibits strong correlation and entanglement, while in the mean-field approximation there is a jump in the order parameter reflecting the change of symmetry. Far from the critical frequency, the mean-field theory can well describe the formation process of vortices.

III. GROSS-PITAEVSKII EQUATION AND CONSERVATION LAWS
For the quantum state of a boson gas, let us define Ψ † (r) and Ψ(r) as the field operators to create and annihilate a particle at point r respectively. These operators satisfy the bosonic commutation relation, The ensemble average ⟨Ψ † (r)Ψ(r)⟩ gives the expectation value of the number density n(r) with the total number given by N = ∫ drn(r). If the interacting range r 0 and the s-wave scattering length a s are much smaller than the average distance between particles d = (N V ) −1 3 , the binary interaction will dominate the dilute and weakly interacting gas. The dynamics of the system is controlled by the Hamiltonian operator, where ̵ h is the reduced Planck constant, m is the particle mass and V (r −r ′ ) is the two-body potential. In the zero temperature limit, almost all particles are in the ground state to form the Bose-Einstein condensation (BEC). Because the total number is huge, the noncommutability of the field operators can be ignored. In this situation, we can replace the field operator Ψ(r, t) with a classical field ψ(r, t) so the Hamiltonian operatorĤ becomes the classical Hamiltonian H. The classical field ψ(r, t) is also called the wave function of the BEC but there is a difference from the normal wave function in the Schrödinger equation or quantum mechanics for single particles in the normalization: in the BEC case, ψ(r, t) is normalized to the particle number N, while in the latter case, the wave function is normalized to 1. In the rest of the paper, ψ(r, t) will be called either the classical field or the wave function: both are equivalent.
For the behavior of the classical field over a longer distance than the scattering length, the details of the potential are not relevant any more. We can approximate the potential by a contact form V (r − r ′ ) ≈ λδ 3 (r − r ′ ), as long as the Born approximation is applicable and the coupling constant λ can reproduce the s-wave scattering length, λ = 4πa s ̵ h 2 m. So the Hamiltonian for a BEC at zero temperature can be written as, The classical field ψ(r, t) varies slowly in time and space and satisfies a non-linear Schrödinger equation, the so called Gross-Pitaevskii Equation (GPE) From this equation, it is easy to obtain conservation laws about the particle number, energy and momentum where ρ, j, H, ⃗ P, mj and Π are the particle number density, the particle number current density, the energy density, the energy current density, the momentum density and the momentum current density, respectively. The classical field can be parameterized in terms of the number density ρ = ψ 2 and the complex phase angle θ(r, t), which is called the Madelung representation. In the Madelung representation the quantities in conservation laws can be put into a compact form with the help of two vectors, In the particle number conservation law (5), the particle current is defined as j = ρv. In the energy conservation law (6), the energy density has the form, which is composed of a "classical" kinetic part, a "quantum-pressure" part and an interacting part 82 . Correspondingly, the energy current density is given by, In the momentum conservation law (7), the momentum density is given by mj, while the momentum tensor density is given by, which is symmetric for interchange of two indices. With these ingredients defined, Eq.(5,6,7) make up a quantum hydrodynamic description for BECs with spin-0 Bosons. For BECs with spin-1 bosons, a quantum hydrodynamic model is derived in 83 .
For quantum hydrodynamics, the variation in density is not negligible if the "quantum pressure" is comparable to the interaction energy density. The balance between them determines the healing length ξ = ̵ h √ 2mλρ, over which a perturbation in density can relax to its bulk value. When ξ ≪ D with D being the characteristic length scale of density variation, the Euler's equations for the irrotational compressible perfect fluid can be reproduced by just setting K = 0 in Eqs. (11,12,13). From the Euler's equations, we know that the interaction energy density also works as the classical pressure p = (λ 2)ρ 2 . The speed of sound in BEC is determined by the compressibility of the fluid, .
The fluid can be considered incompressible if the velocity is small in comparison with c s .
In the matter wave picture, the flow velocity is proportional to the gradient of the wave function's phase θ which changes 2π over de Broglie wave length λ D , so we have v ∼ ( ̵ h m)2π λ D . If the flow velocity is comparable with the speed of sound, the particle number density may change dramatically, so that the destructive interference of matter waves are observable at the typical length scale λ D ∼ 2 √ 2πξ.
It is interesting to make comparison of these length scales, Therefore the condition λ D ≫ ξ ≫ d ≫ a s holds for the dilute Boson gas in the weakly interacting limit.

A. Conservation of OAM and Vorticity
The Lagrangian of the system can be derived by a Legendre transformation from the Hamiltonian (3), The GPE (4)  The conservation equation for the OAM reads, where M 0i and M ij (i, j = 1, 2, 3) are the OAM density and the i-th component of the OAM current density in the j direction, respectively. They are given by, where j = ρv is the current density of the particle number and Π il is the momentum tensor density given by (13). The OAM density vector is then l = m(r × j).
We can also rewrite the OAM conservation equation (18) into the conservation equation where the vorticity density M 0j and vorticity current density M ij in the j direction can be obtained by replacing r → (1 m)∇ and respectively. They are given by, Eq. (19) can also be derived by taking a curl of Eq. (7) or directly from the GPE. The curl of the momentum current stands for the local density of angular momentum or the vorticity density, which leads to conservation of the total vorticity ∫ d 3 r∇ × j. We can rewrite, then Eq. (19) can be put into a different form, or in the Madelung representation, If the interaction is not in a contact form, there will be an additional term from the in- be the divergence of the number density.

IV. STATIONARY VORTEX SOLUTION TO GPE
For a rotating BEC, there may exist vortex solutions to the GPE. The classical field has a rotation symmetry alone the vortex string, and can be parameterized in the following form, where (r, φ) are polar coordinates in the transverse plane perpendicular to the vortex string, with r being the distance of the field point from the vortex string and φ being the azimuthal angle, s is the winding number which must be an integer to guarantee the single-value wave function, it is required that ρ = 0 at r = 0, and the time dependence of the wave function is fixed by the chemical potential µ = δH δN.
From the field form Eq. (25), we can read out the velocity v = ( ̵ hs) (mr)φ, and then we obtain the OAM density l = ρ(r)s ̵ hẑ, which indicates that each particle contributes s ̵ h to the total OAM, no matter how far it is from the vortex string. The number density varies strongly in the range of the healing length ξ near the vortex string. For convenience we use the profile function f (η) = ρ(η) ρ 0 for the field amplitude where ρ 0 is the density in the bulk and η ≡ r ξ is the rescaled radius. Inserting Eq. (25) into Eq. (4), we obtain the with boundary conditions f (η → 0) = 0 and f (η → ∞) = 1. To satisfy boundary conditions, it is easy to see µ = λρ 0 , then the profile function can be numerically solved. A shooting method is employed to solve this non-linear differential equation. By tuning the first order derivative f ′ (η) at η = 0, we can find the solution which satisfies the condition f (η → ∞) = 1.
We have checked the result with a variational method and a semi-analytical method, all of them can precisely reproduce the profile function for s = 1, 2, 3.
In Fig. 1, we show the profiles of the field, vorticity density and energy densities as  Fig. 1d, For the primary vortex, these quantities show a core structure in the profiles, in other words, they have maxima at the vortex center. While for vortices with s > 1, the profiles of these quantities show a ring structure, which vanish at the vortex center and reach maxima at finite radii. These differences indicate that the primary vortex is topologically different from vortices with larger winding numbers. The primary vortex looks like a particle in the ground state in a rotational BEC, while the vortices with larger winding numbers look like particles in exciting states and will decay into the primary one when they are perturbed by fluctuations 84 . This picture can be better illuminated by an analysis of the energy and OAM carried by vortices 84 .
In grand canonical ensemble, the relative energy of a vortex in the BEC in a cylinder with the volume V = πLR 2 is given by, where H(r) is given in Eq. (11). There is a logarithmic divergence from the kinetic energy, so we introduce a radius cutoff R ≫ ξ to regulate the integral. With the profile solution f (η), the result takes a compact form, where the factor γ s = 1.464, 0.572, 0.358 for s = 1, 2, 3 respectively. For each individual vortex, the relative energy is proportional to s 2 , while the OAM is proportional to s, so vortices with larger winding numbers are not stable and will decay to primary vortices. Even the primary vortex is an excited state of BEC because of its positive energy. Only in a rotating frame with a considerable larger angular velocity ω or in a system with a sizeable OAM can the vortex solution correspond to the global or local minimum of the free energy defined as To know how vortices evolve with time, we need to solve the GPE (4) with proper initial conditions and period boundary condition 85 . We choose the following initial value for the classical field, which represents an array of vortices or rotating fields as shown in Fig. 2 for N  does not have such a periodic feature for finite N 0 , but the difference between gradients at two corresponding points on the boundary can be suppressed significantly by increasing N 0 .
With N 0 ≥ 3, the difference in gradients is negligible.
In the numerical simulation, we set the parameters as ̵ h = 1, m = 4, λ = 0.002, and ρ 0 = 1. With these values of parameters the healing length becomes ξ = ̵ h √ 2mλρ 0 ∼ 7.9 and the chemical potential becomes µ = λρ 0 = 0.002 which set the scales in space and time respectively. We employ a symplectic algorithm at sixth order 86,87 to perform the time evolution on GPU, and find that the simulation with grid size a = 1 and time step ∆t = 0.5 can maintain the energy conservation and time reversibility in a high precision. A brief introduction to symplectic integration algorithm is given in Appendix A.
To look at the time evolution of a single vertex, we have to minimize the interference from other vortices. So we set the block size to be much larger than the healing length A ≫ ξ. In short time, each block can be considered as effectively isolated from others, so its evolution is almost independent. In Fig. 3 2 2 2 2 The decay process can also be demonstrated by the kinetic energy density and vorticity.

VI. VORTEX PRODUCTION IN NON-CENTRAL COLLISIONS
In this section, we discuss the production of vortices in non-central collisions of two BEC fields as wave packets. Experimentally, this setup can be realized by colliding two clusters of cold atoms at extremely low temperature. In our numerical simulation, we split a space area of 518ξ × 518ξ into 4096 × 4096 grids, each grid has a size 0.1265ξ × 0.1265ξ. Similar to nuclei, we choose the initial density of the field to be in the Woods-Saxon distribution where R 0 is the radius of the disc and a is the skin thickness. In our numerical calculation, we choose R 0 ≈ 60ξ and a ≈ 1.26ξ to guarantee the reliability of the results. Due to the repulsive quantum pressure and self-interaction, an isolated disc is not stable on its own. This can be seen in the time evolution of the density profile in Fig. 7: the field expands isotropically and the central density drops to two-thirds of its initial value at t = 45 ̵ h µ.
We consider head-on collisions of two field configurations of disc shape with velocities v 0 e x and −v 0 e x at the initial time. We choose the impact parameter (5 6)R 0 , i.e. the distance between the centers of two discs in the y direction. The initial field configuration then has the form,  With our choices of R 0 and v 0 , the two fields collide before significant expansion. After the collision, the system expands due to the repulsive quantum pressure and the self-interaction.
There is a net radial flow of fields from the collision region to infinity, together with the produced primordial vortices. In the simulation, the periodic boundary condition is automatically guaranteed before the time when the outgoing wave front arrives at the boundary.
We take a very large lattice, so that the formation of vortices and vortex-antivortex pairs in the central region is not affected by the boundary condition.

A. Density distribution
Because of the quantum nature, two colliding wave packets interfere when they approach to each other. The interference fringes in the number density distribution are shown in

B. Conservation of the winding number
During the formation of vortices, the velocity field and the density distribution of the particle number, energy, OAM and vorticity change dramatically, but the total winding number as a measure of the topological defect is conserved. The winding number at each point is determined by the winding integral of the local velocity as, It takes a nonzero integer when a defect is present inside the area of the integral. The sign of the integer is related to the local vorticity. The position is controlled by the surrounding environment.
As shown in purple circles in Fig. 9, there is a defect in the symmetry center of the collision. The position of the defect is fixed and does not move with time because of the symmetric initial condition. Its winding number is always −1 even when the surrounding densities change dramatically. We will focus on this defect and study the dynamical formation of vortices in later subsections.

C. Concentration of energy and vorticity
During the formation of a vortex, the energy and vorticity are being concentrated towards the position of the topological defect or the vortex. Fig. 10 shows the time evolution of the energy density and vorticity density distributions in the central region as marked by the red square in Fig. 8. At t = 0, the main contributions to the energy density are the kinetic At a later time, there is a dramatic redistribution of conserved quantities, such as the particle number, energy, OAM and vorticity in the overlapping region. At t = 12, Fig. 10 shows richer structures than the number density in Fig. 8. There is a distinct belt of the energy density in the center of the destructive fringes where the number density is small but θ and ρ change fast on the belt, so that the magnitudes of v and K are large. The separation of the positive and negative vorticity is visible: the negative vorticity region is connected along the belt covering the symmetry center while the positive regions are disconnected.
The difference becomes more obvious at t = 24. As time goes on, the energy density and the negative vorticity density are being concentrated towards the symmetry center marking the formation of a vortex. At t = 48, the formation of the vortex is almost completed with its profiles in the energy and vorticity density satisfying those of the vortex solution in Fig. 1, while the bulk value of the number density in the vicinity is about ρ ρ 0 ∼ 1.0.

D. Evolution of average densities of particle number, OAM and vorticity
To see how conserved quantities evolve to values of the vortex, we can look at the time evolution of the average densities of these conserved quantities. The evolution is governed by conservation laws Eqs. (5,6,7,19). We will focus on the areas inside the purple circles with the radius R = 6ξ in Fig. 10 and look at the changes of conserved quantities.
The winding number weighted by the number density is defined as which is also the ratio of the total vorticity inside the circle to the average number density on the boundary times 2π ̵ h. We see that ⟨s⟩ D matches s in Eq. (33) for a density distribution with the rotation symmetry.
We see that the number density (black-dashed line) increases, saturates, and then de- in the interval 10 < t < 32. It is interesting to notice that there is a peak around t = 40, which comes from the breaking of the destructive interference fringe at the boundary by filling up particle density in the lower density region except the center of the vertex. After that, the average number density decreases continuously due to global expansion.
The time evolution of the average OAM inside the circle marked in Fig. 10 is shown by the blue-dotted line. We see that its magnitude increases in t < 5 due to the asymmetric inflow of particles from the positive and negative x direction. This part of the OAM comes from the initial geometry. During 5 < t < 10 the inflow of particles from both directions becomes almost symmetric, the average OAM is almost a constant. When the dynamical balance between the pressure and the flow is constructed, there is an exponential increase in the magnitude of the OAM during 10 < t < 40. The exponential increase stops when the destructive interference fringe breaks at about t = 40 when the magnitude reaches a maximum. At the same time, the OAM per particle ⟨l z ⟩ ⟨ρ⟩ (solid-red line) saturates at the value 1 which indicates that each particle carries one unit ̵ h on average. After t = 40 the OAM per particle remains a constant even though the particle number and OAM densities decrease continuously in global expansion.
The green-long-dashed line shows the average winding number in the circle as a function of time which has a similar feature to the average OAM per particle. Both indicate the collective motion inside the circle that leads to the formation of a vortex. We have tested the saturation to 1 for both the average winding number and the average OAM per particle with different radii smaller than a half distance between two topological defects. The saturation in a smaller circle is fulfilled at an earlier time, which means the formation of a vortex starts from the topological defect. As time goes on, more and more particles are involved in a vortical rotation, the effective size of the vortex increases with a rate determined by the particle number density and OAM current nearby.

E. Vortex-antivortex pair production
During the dynamical formation of vortices in non-central collisions, the excess energy is released by emission of sound waves. The front of an outgoing sound wave will generate a radial flow and interfere with the remnant initial flow just outside the formation region.
As a consequence, the destructive fringe start to break. As shown in Fig

F. Remarks about simulation
To check the precision of the numerical algorithm, for each final state of the collision a time-reversal evolution has been performed, and the initial field configuration can be reproduced correctly. We have also inverted the direction of the particle's velocity by taking the complex conjugate of the final state field, from which a proper time evolution can also reproduce the initial field configuration. This indicates that the numerical viscosity is under control and is negligible in the process of vortex production.
The production of vortices and pairs of vortices are very sensitive to the initial velocity.
The larger the initial velocity is the more vortices and pairs of vortices are created in a shorter time. For smaller initial velocities, the formation of vortices depends on the competition between rotation and expansion. If the initial disc is big enough to suppress the expansion for a long time, the central vortex is more likely to be observable, because the central topological defect need more time to absorb enough OAM from the outer environment. With smaller initial velocities and a smaller disc, it is harder to form a notable vortex due to the lack of enough time and OAM.
We have tested the situation with opposite initial OAM by inverting the impact parameter of the collision. In this case, the winding numbers of single vortices (not vortex-antivortex pairs) are flipped. In central collisions with zero OAM in the initial state, there are no single vortices, but the vortex-antivortex pairs can be created if the colliding velocity is large enough.

VII. CONCLUSIONS AND DISCUSSIONS
We present a different approach to spin polarization from conventional ones through quantum vortex formation in collisions of BEC. This approach is based on the observation that the vortex is a topological excitation in a superfluid in presence of local orbital angular momentum and is an analogue of spin degrees of freedom. Therefore the formation process of quantum vorticies in collisions of BEC may shed light on the nature of the particle's spin as well as spin-orbit and spin-vorticity couplings in the strong-interaction matter.
We have investigated the stationary property and the dynamical formation of vortices (including vortex-antivortex pairs) by numerically solving the GPE with a large-scale parallel algorithm on GPU to extremely high precision. We have demonstrated that the primary vortices can be considered as basic degrees of freedom in a system with a sizable global OAM in the initial state. The OAM is present in a rotational BEC or in non-central collisions of two BEC fields. The energy and vorticity density will nucleate around topological defects, leading to the formation of vortices in short time.
In a rotating environment, each particle carries the same local OAM, no matter how far it is from the rotation axis. An isolated vortex with a given winding number prefers to keep its own stationary properties. The density distribution of the vortex around its center responds quickly against perturbation, trying to restore its stationary profile. Even for a vortex with winding number s > 1, the stationary density profile can also be restored in short time. But in long time, a s > 1 vortex is unstable and subject to decay into primary vortices ( s = 1) which are the lowest energy states. The decay can be seen in the splitting of a topological defect. After that, all primary vortices are being decoupled while rotating around each other, and finally form an array of vortices. The decay precess takes much longer time than the relaxation time in the profile-restoration process.
We have observed the formation of vortices in non-central collisions of two BEC discs. A global OAM is present in the initial field configuration. If the initial global OAM is large enough, vortices are generated. The primary vortices are first generated in the central region of the collision as the result of a part of OAM is transferred to vortices. The winding number of each topological defect is conserved in the field evolution. Around a topological defect, the densities of the particle number, energy and vorticity are tuned to match the profile of the primary vortex: the particle number density reaches its local minimum while densities of the energy and vorticity (absolute value) reach their local maxima. With the emission of sound waves, an array of primary vortices are finally formed. This process can also be considered as the localization of a global rotation. If observed in a large area, vortices can be considered as fundamental objects like spin degrees of freedom. Its formation can be viewed as a consequence of the spin-orbit coupling.
In the follow-up evolution, the interaction between sound waves and interference fringes may generate vortex-antivortex pairs. The creation of vortex-antivortex pairs is not sensitive to the impact parameter, but strongly depends on the collision velocity. A vortex-antivortex pair does not carry a local OAM but carries energy and momentum in a finite region. The pair creation process in collisions is a result of redistribution of local energy and momentum densities.
The decay precess breaks the time-reversibility, i.e., when the primary vortices are decoupled, the time-reversal simulation from the final state can not recover the initial state. The reason is unclear by now, because we can not distinguish quantum decoherence from the numerical viscosity. But in short time simulation, the coherence and time-reversibility of the field are maintained to a high precision because the initial state can be well recovered by the time-reversal simulation from the final state. We also checked the influence of the numerical viscosity by inverting the direction of motion in the final state followed by a proper time evolution, and found that the initial state configuration can be recovered. This means that the numerical viscosity is under control in short time evolution.

Appendix A: Symplectic integration algorithm
In this appendix we give a brief introduction to the symplectic integration algorithm to solve the GPE (4). The formal solution to the GPE can be put into the form, whereÛ is the evolution operator depending on the coordinate and momentum operatorsx andp as follows,Û The formal solution is not ready for numerical implementation because there is a time ordering operator T and operators are not commutable. The main difficulty is due to eÂ +B ≠ eÂeB, so we cannot perform each step of the evolution independently. There is a approximation called "splitting operators" to replace the coupled operators with a bunch of separated operators, which is easy to perform and can achieve a very high precision. For example, to O(∆t 7 ), the evolution operatorÛ can be rewritten as, where d α and c α are coefficients precisely designed to achieve the desired precision. Then the evolution from t → t + ∆t is made by a combination of operators acting on the wave function at a given time t. The symplectic integrator has a symmetric form, so that the time reversal symmetry can be preserved with an extraordinary precision,Û (t → t + ∆t)Û (t + ∆t → t) = O(∆t 6 ).
Given the initial wave function in coordinate space with a periodic boundary condition, we first transform it into momentum space, So the effect of the operator exp(−ic α ∆t p 2 2m) is a modification to the phase of the wave function in momentum space, Then we transform the wave function back to coordinator space, The effect of the operator exp(−id α ∆tλ ψ(x) 2 ) is to modify the phase of the wave function in coordinate space, In this way, at each time step ∆t, we have to make Fourier transformation back and forth eight times. In a successive evolution over t, the particle number conservation will be preserved, while the energy conservation will be slightly violated at O(∆t 6 ). A smaller ∆t or a higher order algorithm will help to achieve a better conservation.