Emergent Flow Signal and the Colour String Fusion

: In this study, we develop the colour string model of particle production, based on the multi-pomeron exchange scenario, to address the controversial origin of the ﬂow signal measured in proton–proton inelastic interactions. Our approach takes into account the string–string interactions but does not include a hydrodynamic phase. We consider a comprehensive three-dimensional dynamics of strings that leads to the formation of strongly heterogeneous string density in an event. The latter serves as a source of particle creation. The string fusion mechanism, which is a major feature of the model, modiﬁes the particle production and creates azimuthal anisotropy. Model parameters are ﬁxed by comparing the model distributions with the ATLAS experiment proton–proton data at the centre-of-mass energy √ s = 13 TeV. The results obtained for the two-particle angular correlation function, C ( ∆ η , ∆ φ ) , with ∆ η and ∆ φ differences in, respectively, pseudorapidities and azimuthal angles between two particles, reveal the resonance contributions and the near-side ridge. Model calculations of the two-particle cumulants, c 2 { 2 } , and second order ﬂow harmonic, v 2 { 2 } , also performed using the two-subevent method, are in qualitative agreement with the data. The observed absence of the away-side ridge in the model results is interpreted as an imperfection in the deﬁnition of the time for the transverse evolution of the string system.


Introduction
In recent decades, a significant effort has been made in the study of a unique state of matter called [1] quark-gluon plasma (QGP).The first experimental evidence of QGP formation was claimed by CERN in 2000 [2].In 2005, this statement was quantitatively confirmed by RHIC experiments [3][4][5][6], whose studies raised questions about the properties of QGP.
Among numerous illuminating hints of a QGP signal, there was evidence [4,5] of the azimuthal anisotropy of particles produced in semi-peripheral heavy ion collisions.The measurements were done in a model-independent way [7] using transverse flow harmonics, v n , of different orders, n: directed flow (v 1 ), elliptic flow (v 2 ), triangular flow (v 3 ) etc.A considerable magnitude of v 2 was observed [4], whose contribution is assumed to be a dominant one following the almond shape of the intersection region of two nuclei.The values of v 2 were comparable to the calculations [8] obtained when relativistic hydrodynamical fluid was considered in the early times of the collision.Therefore, it is the collective motion of thermally equilibrated partonic degrees of freedom that is believed to convert spatial anisotropies into momentum space under pressure gradients of the compressed medium.This rejected the view that QGP is a gas, but nevertheless, left still to establish whether this represents a perfect or viscous fluid.
Since then, plenty of flow measures have been proposed that vary in sensitivity to nonflow effects.Those are pair correlations [9], multi-particle correlations [10] measured with cumulants [11] or with subevent cumulant method [12], symmetric [13] and asymmetric [12] cumulants considered in subevent method and a peculiar combination of flow-mean particle transverse momentum correlations [14].
The surprising evidence that connects these studies and the present paper is the experimental observation of collective behaviour in small systems.Namely, the unexpected ridge signal [15,16] and associated flow harmonics [17] seen in high-multiplicity protonproton (pp) interactions at the LHC.These observations should not be explained by the fluid nature of the medium produced in such a small droplet of matter as it should not be able to reach thermal equilibrium prior to hadronization [18].Nevertheless, there are successful attempts [19] to apply hydro description even in this case, which raises questions about the use of the same approach in all colliding systems [20].
On the one hand, from the first principles of the theory of strong interaction, the description of multi-particle production is complicated by the feature that the majority of the particles are produced in soft processes.It means that their transverse momenta, p T , do not exceed about 1 GeV, therefore, the perturbative calculations in quantum chromodynamics (QCD) are inapplicable in this regime.This forces us to work in phenomenological approaches that can effectively describe the transition from, for instance, a few colliding hadrons to hundreds and thousands of particles produced.
One of the successful techniques is based on the concept of formation of colour strings between colliding partons which then fragment into observable hadrons [21].This method appeared to be effective both in phenomenological calculations (see e.g., dual parton model [22], string percolation model [23], colour-glass condensate and glasma model [24]) and as the basis of many Monte-Carlo event generators, such as, e.g., EPOS [25], PYTHIA [26], HIJING [27], AMPT [28] models.In this paper, we consider the colour string picture of multi-particle production to model the transverse flow effect measured in pp collisions.It is the strings' interaction in the form of fusion [23] that plays a primary role in the model.
Our interest in this topic is inspired by the approach developed in the series of papers [29][30][31] in which the flow signal appears due to the quenching of particle transverse momentum in a string medium [30].Namely, it is assumed that particles emit gluons while passing through the strings, which is similar to the energy losses of charged particles moving in the electromagnetic field in quantum electrodynamics (QED).Hence, the appearance of forbidden azimuthal angles changes the distribution of particles in the transverse plane.Here, the fusion of strings creates a heterogeneous medium density along the path of a particle.
In addition, we consider that the fusion accelerates the overlapped strings, resulting in a boost of the produced particles.It creates distinct directions of motion in the transverse plane for particles originating from the rearranged colour field.This approach originates from the pioneering paper [32] and has already been tested in more recent models [33,34].The interplay of these two mechanisms, driven by string fusion, allows one to obtain specific joint correlations in pseudorapidity and azimuthal angle, η and φ, respectively, of the same origin in collision systems of different types.
In this paper, we extend our model [35,36] that provided a thorough description of particle correlations and fluctuations in η, for inelastic pp interactions at top RHIC and LHC energies.In this paper, we follow up [37] and test the model application to the problem of describing flow in pp interactions, aiming to obtain characteristic η-φ correlations at the pp system centre-of-mass energy √ s = 13 TeV.
In Ref. [36], we emphasized that the previous model implementation lacks the shortrange correlations introduced.It did not allow us to fully describe the ALICE experiment data, but we were able to catch the trend.In the present investigation, on the contrary, we aim to take advantage of this drawback of the model since getting rid of non-flow effects is the main challenge in flow studies.However, there is still a need to partially modify the model.
In the present study, the basic model features, as in Ref. [36], are the 3D (threedimensional) dynamics of strings and string fusion.The strings consideration is slightly modified, but the main point stays intact: longitudinal contraction and shift of strings in the rapidity space, inspired by [38], and their transverse motion governed by an attractive Yukawa potential via the σ-meson exchange, following [39].
The main set of changes introduced in this study concerns the string fusion mechanism [40,41].Namely, we abandon the consideration of the simplified cellular approach [42] to find the overlapped strings by looking after their centres.Instead, we divide the transverse plane into pixels that are much smaller than the string's size.The contribution of a pixel to particle production is determined by the degree of strings fusion that occurs inside it.We apply these changes with an eye toward future studies of nucleus-nucleus (AA) collisions.
It is worth noting separately that the mechanisms of momentum loss and string boosts applied to AA collisions provide a good description of the observed flow signal [43,44].However, the additional introduction of the transverse motion of strings would naively result in the formation of a single colour flux tube after string fusion is taken into account.Therefore, the problem with the previous implementation of the model [36] is that a modified but uniform colour field could not serve as a source of azimuthal anisotropy.
The approach with the fine granularity of the configuration-momentum space that we propose in the present paper solves this issue, but we first test it in the description of pp collision results.Although the interpretation of pp results on azimuthal anisotropy provokes some tension in the flow community, it is technically more convenient for us to start the model development with the small colliding system, since keeping the fine cellular structure is much more challenging in AA, CPU-time wise.Furthermore, there are other improvements of the model mostly at the stage of proton assembly and event simulation implemented in this work.
To summarize, we develop the Monte-Carlo model based on the colour strings approach that addresses the problem of azimuthal particle correlations via the string-string interactions and particles' momentum loss.Thus, we anticipate that our preliminary estimates will demonstrate whether the nature of collective effects observed in pp collisions can be revealed using such a simplified but authentic model.
The paper is organized as follows.Section 2 presents the model framework, with special attention given to improvements introduced to the model compared to its previous version [36].The description of the flow mechanism implemented in the model is included in a separate Section 3. In Section 4, we demonstrate the inference of free model parameters based on comparison with ATLAS experiment pp collision data at √ s = 13 TeV [45].
Section 5 introduces the flow measures that we calculate in the model.Results obtained are presented in Section 6 and discussed in Section 7.

Colour String Model
The colour string model [21] is a phenomenological approach that is used to describe the soft particle production regime where the perturbative QCD calculations are not applicable.Strings that are extended colour field objects (also called tubes of gluon field) are stretched between partons participating in the collision.We consider them all to be parallel with neither kinks nor twists.It is believed that strings have a finite size in the transverse dimension, which enables them to overlap and interact.String longitudinal dynamics is described by a yo-yo solution, which means that a string periodically stretches and contracts with time.The energy of the colour field inside the tube grows as the string extends, while its massive endpoints slow down.It causes a probabilistic break-up of the colour field via the creation of quark-antiquark, q-q, pairs and the breaking of strings into observed particles.

Multi-Pomeron Exchange in pp Collisions
In this paper, we consider inelastic pp interactions at √ s = 13 TeV in the context of multi-pomeron exchange, with each pomeron represented by a cylindrical diagram.The uncut diagrams contribute to the elastic cross-section, while their unitarity cut creates two colour strings [22] that fragment into observed particles.Therefore, in this approach, the number of strings, n str , in an event is determined by the number of exchanged pomerons, n pom , as n str = 2n pom .
With the natural assumption of the Gaussian distribution of the transverse positions of partons inside the proton, one can calculate a probability of parton-parton interaction that can be interpreted as the probability of string formation [46].The resulting distribution of the number of pomeron exchanges coincides with the Regge-theory parametrization [47] that neglects the three-pomeron vertices: where 5 is the quasi-eikonal parameter related to the small-mass diffraction dissociation of incoming hadrons, ∆ = α(0) − 1 = 0.2 is the residue of the pomeron trajectory, α(0) is the intercept of the pomeron trajectory, γ = 1.035GeV −2 and R 2 = 3.3 GeV −2 characterise the coupling of the pomeron trajectory with the initial hadrons, α = 0.05 GeV −2 is the slope of the pomeron trajectory.
Figure 1 presents the event distribution in the number of strings, n str , obtained with these parameters and used in our simulation.To determine the values of parameters in Equation ( 1), the data on multiplicity distributions and cross sections were used.We substitute into Equation (1) the parameters from Ref. [48], where it is assumed that strings that were initially formed in an event can overlap to various degrees and fuse, forming the string clusters with the modified fragmentation characteristics, which changes the obtained multiplicity.This is in contrast to another approach used in Ref. [49], where string fusion is effectively implemented already at the moment of string formation.Therefore, there are a number of free strings in Ref. [48] and a number of string clusters in Ref. [49].The values of parameters in Equation (1) for these two methods will differ, but both the models should lead to the same charged particle multiplicity distribution, P(N ch ) (see Equation (2) below), corresponding to the data.We follow the first approach that allows us to track the 3D evolution of string density.Thus, the proton no longer consists of a single quark and a diquark and forms not only two strings in pp interaction, as we move away from the conventional picture.Instead, we take into account the multi-parton configurations of sea quarks whose independent and simultaneous interactions result in the formation of the number of pomerons in an event.
In this approach, P(N ch ) is determined as a convolution of P n pom (1) and the multiplicity distribution at fixed number of pomerons, P n pom (N ch ), P n pom P n pom (N ch ). ( Due to the complexity of the analytical formulation of P n pom (N ch ), we treat it using a Monte-Carlo simulation.In Ref.
[35], we demonstrated that for quite a simple case of non-interacting strings, the correlation observables calculated numerically and analytically are in good agreement.

Proton Assembly
Our current study is based on the assumption that all partons of a proton form strings with all partons of another proton.Therefore, we call an event the creation of n str between two protons, each with the number of partons, n part , which is equal to the number of strings in an event, n part = n str .
To fulfil this requirement, we prepare an extensive set of protons with all the possible even numbers of partons, n part = 2k, k =1 , . . ., 50.In this study, we consider the parton composition of a proton as one valence quark, one valence diquark, and (n part − 2) sea quarks.The algorithm provided in Appendix A is more time-and CPU-efficient than the permutations of partons between two selected protons that were implemented in Ref. [36].

Event Simulation
The first step of a simulation of an event is to sample the number of exchanged pomerons, n pom , from Equation (1).Thus, we know the number of strings that should be created in an event, n str = 2n pom .By creating the event, we mean searching for two protons, with the same number of partons, n part = n str , so that all their partons can form strings under certain conditions (see Sections 2.3.1 and 2.3.2).
First, we select two random protons from the prepared set (Section 2.2) with the certain n part .Second, we permute partons in one of the two protons by performing n part replacements and checking whether this leads to the formation of n part good pairs (see the requirements in the next two subsections).If not, we select another pair of protons.

String Formation: Electric Charge
In this study, we ensure the conservation of electric charge.Thus, we accept a string candidate if the electric charge of the string, Q str , defined as the sum of the partons' charges that form a string, equals one of the possible integer numbers: −1, 0, +1, +2.

String Formation: Energy and Mass of Decay Products
We calculate the string energy, E str , by summing the energies of partons at the ends of the string, E part1 and E part2 , as where p part is the initial momentum of a parton and m part is a dynamically defined parton mass (for details of their definitions, see Appendix A.3).We accept the string candidate if E str is sufficient to decay at least in two hadrons at rest, with masses M daughter1 and M daughter2 , i.e., E str ≥ M daughter1 + M daughter2 .
In order to test this condition, it is necessary to identify the species of the pair of daughter particles based on the flavours of the string ends, e.g., the quark-diquark string should decay at least to a baryon and a meson.For completeness, we provide the list of minimum permitted combinations of daughter products, based on the flavours of quarks that we consider in the model: We use m nucleon = 0.9396 GeV, m π = 0.1396 GeV, m K = 0.4977 GeV, m D = 1.8696GeV for nucleons, pions, kaons, and D-mesons, respectively.

String Transverse Dynamics
In general, the colour confinement in QCD, a non-abelian gauge theory, is viewed as the gathering of the colour field between two colour charges in the flux tube of finite transverse size [50].However, lattice QCD demonstrates that the presence of a colour string modifies the QCD vacuum.For example, the correlator between the quark-antiquark chiral condensate, q q , and the Wilson loop, W, is appeared to be not a constant as a function of the transverse distance from a string [51].Indeed, at large distances it reaches the value of unity, meaning that there is no string influence.In the meantime, its values decrease in the vicinity of the string, which indicates a partial restoration of chiral symmetry in this region of the space.
We follow the approach of Refs.[39,52], where the authors interpret these lattice calculations (left-hand side of Equation ( 4)), as a field created by a cloud of σ-mesons (right-hand side of Equation ( 4)): where r ⊥ = r 2 ⊥ + s 2 str is the regularized distance in the transverse plane, r ⊥ is a 2D distance between strings in a pair of strings, s str = 0.176 fm [39] is a genuine string width, unlike the effective string width, which is a result of quantum fluctuations, K 0 is zero-th modified second-kind Bessel function corresponding to a massive scalar propagator in two dimensions and m σ = 0.6 GeV [39] is the mass of the σ-meson that is proposed to be a mediator of the force between strings.
In this paradigm, one can consider string-string interactions in the created Yukawa potential, which results in the non-relativistic attraction between them similar to nuclear forces.Thus, the equations of motion of the string system in an event are defined by the 2D Yukawa interaction [39], ¨ where r ij and rij correspond to r ⊥ and r ⊥ , respectively, from Equation (4).The i and j subscripts indicate that the quantities are constructed for the i-th and j-th strings.Here, g N is the QCD string self-interaction coupling and σ T is the string tension, whose dimensionless product is selected as g N σ T = 0.2 [52], and K 1 is the first modified second-kind Bessel function.Strings are considered to be moving as a whole according to Equations (5), without whatever kinks.To solve this system of differential equations, the set of initial conditions is required.The Gaussian distribution of width 0.5 (model parameter) is used to sample the initial 2D transverse coordinates of string centres in an event.This simplification is to decrease the program's running time: instead of applying the Glauber approach at the partonic level and finding which partons are to form strings, we assume that some configuration of strings has already been created in an event.
Transverse strings' evolution, governed by Equation (5), can be terminated at some proper time, τ stop , which affects the final string density.In Ref. [36], we showed that the largest number of overlaps between strings occurs after the time that we called τ deepest .The latter is the time it takes for a string transverse configuration to attain the global minimum of the potential energy, evolving according to Equation (5).Ref. [36] considers other cases, including no transverse evolution and transverse evolution until the conventional time τ stop = 1.5 fm [52] that is unsuccessful in describing the data, especially the p T -N correlation, where p T is the event-averaged transverse momentum of a particle and N denotes the particle multiplicity.Thus, in the present study, we consider τ stop = τ deepest .It varies event-by-event and depends not only on the initial transverse positions of the strings, but also on the number of strings, and thus on the collision energy.It is worth noting that τ deepest cannot be found for any event.This means that for some configurations of strings, the minimum of the potential energy of the system can only be reached at τ stop > 1.5 fm.In this case, we set τ stop = 1.5 fm as the typical time for string hadronization and restrict system evolution to it.
Figure 2 shows how a few events (one event a row) with different numbers of strings look before and after the transverse evolution of string density.Figure 2 (left half, left to right) demonstrates the event projections to the transverse plane, X-Y, before and after the transverse dynamics of strings.Figure 2 (right half, left to right) present the projections to the X-rapidity plane before and after the transverse evolution.The corresponding time of the transverse motion, τ transv = τ deepest , is indicated as relevant.One can see that the string system is highly compressed and the 2D distribution (Figure 2, second left) approaches a more spherical shape in comparison to the initial positions of the strings (Figure 2, most left).In the rapidity dimension, one obtains quite uniform distribution up to large rapidities (Figure 2, most right).

String Longitudinal Dynamics
In the model considered here, special attention is given to the simulation of partonic degrees of freedom for colliding protons (see Appendix A).It is because the initial rapidities of the string's ends, y part init , are determined by using the momenta of partons that form it as However, the string tension, σ T , slows down the massive partons flying outwards with momentum p part according to dp part /dt = −σ T , where t denotes the time.As a result, the rapidity of the string end decreases [38] by value where τ long is the proper time of the string longitudinal evolution.
The aim is to relate the τ long from Equation ( 7) with the time for the transverse evolution, τ transv , to synchronise the dynamics of the string system in rapidity and X-Y dimensions.However, one has to take into account that the compressions and stretchings of a string are periodic (yo-yo solution [21]) and are followed one by another.Moreover, the movement of massive endpoints of different masses, m part1 and m part2 , in the case here is not symmetrical (we denote by 1 or 2 one of the string ends).Therefore, we define the maximum proper time for each string end, τ part1,2 max , after which the string's endpoint fully stops the deceleration (acceleration) after acceleration (deceleration) as where ỹpart1,2 init are the initial rapidities of the string ends converted to the string rest frame.The simulated events (one event a row) with 2, 16, 38, and 74 strings (top to bottom) for an event projection to the transverse plane, X-Y, before and after the transverse evolution and to the X-rapidity plane before and after the transverse evolution (left to right) for the time of the transverse evolution, t transv , as indicated.The evolution runs till t transv = t deepest (second left and most right).Z axis is not to scale.See text for more details.

String Longitudinal Dynamics
In the model considered here, special attention is given to the simulation of partonic degrees of freedom for colliding protons (see Appendix A).It is because the initial rapidities of the string's ends, y part init , are determined by using the momenta of partons that form it as The simulated events (one event a row) with 2, 16, 38, and 74 strings (top to bottom) for an event projection to the transverse plane, X-Y, before and after the transverse evolution and to the X-rapidity plane before and after the transverse evolution (left to right) for the time of the transverse evolution, τ transv , as indicated.The evolution runs till τ transv = τ deepest (second left and most right).Z axis is not to scale.See text for more details.
To convert the rapidity of a string end from the laboratory frame, y part init , to the string rest frame, ỹpart init , one has to know the rapidity of the string in the laboratory frame, y str .It is calculated as follows.First, let us define strings' momentum, p str = p part1 − p part2 , and strings' energy, E str (3).Thus, the rapidity of a string, y str , is As soon as y str is known, one can recalculate the rapidities of the string ends in the string's rest frame: ỹpart1,2 init = y part1,2 init − y str (10) and substitute Equation (10) into Equation (8).
It is necessary to account for the periodicity of string motion.Namely, after the time 2τ part1,2 max (8), the sign of dp part1,2 /dt flips.Therefore, it is crucial to correctly connect the τ deepest , which is found from the transverse dynamics for the whole event, with τ part1,2 max , which varies for two ends of the string and for each individual string in the event.
To summarize, the time of the movement of a string end only during the last period should be used in Equation (7).The rapidities of the string ends after such loss, y part1,2 fin , are found by subtraction of y part1,2 loss from y part1,2 init with the correct signs.It is worth remarking that the longitudinal evolution changes not only the length of the strings but also their positions with respect to midrapidity.Figure 3 demonstrates how the string density of the model events from Figure 2 (already after the transverse evolution) changes after the longitudinal losses are taken into account.For clarity, we repeat Figure 2 (most right) in Figure 3 (most left).In Figure 3 (left half, left to right), one can see the event projections on the X-rapidity plane before and after the longitudinal dynamics.Figure 3 (right half, left to right) shows the Y-rapidity plane projections also before and after the longitudinal dynamics.for an event projection to the X-rapidity plane before and after the longitudinal evolution and to the Y-rapidity plane before and after the longitudinal evolution (left to right) for the time of the longitudinal evolution, t long , as indicated.The evolution runs till t long = t deepest (second left and most right).Z is axis not to scale.
It can be seen that the string system in the final state (Figure 3, second left and most right) has a significant compression in the longitudinal direction.The cloud of a high density of strings is no longer infinitely extended as it is often assumed in the string models [53]; thus, our model loses the translational invariance in rapidity.The transverse evolution of the strings, as described in Section 2.2, also leads to varying string density, which makes an event strongly heterogeneous.

String Fusion
As a result of the reductions and shifts of strings in the rapidity dimension (Section 2.5), the system of strings is a spaghetti-like structure [39], although the lengths and positions of strings vary with respect to mid-rapidity.Furthermore, due to the attractive motion of strings in the transverse plane (Section 2.4), the clusters of fully or partially overlapped The simulated events (one event a row) with 2, 16, 38, and 74 strings (top to bottom) for an event projection to the X-rapidity plane before and after the longitudinal evolution and to the Y-rapidity plane before and after the longitudinal evolution (left to right) for the time of the longitudinal evolution, τ long , as indicated.The evolution runs till τ long = τ deepest (second left and most right).Z axis is not to scale.
It can be seen that the string system in the final state (Figure 3, second left and most right) has a significant compression in the longitudinal direction.The cloud of a high density of strings is no longer infinitely extended as it is often assumed in the string models [53]; thus, our model loses the translational invariance in rapidity.The transverse evolution of the strings, as described in Section 2.2, also leads to varying string density, which makes an event strongly heterogeneous.

String Fusion
As a result of the reductions and shifts of strings in the rapidity dimension (Section 2.5), the system of strings is a spaghetti-like structure [39], although the lengths and positions of strings vary with respect to mid-rapidity.Furthermore, due to the attractive motion of strings in the transverse plane (Section 2.4), the clusters of fully or partially overlapped strings are formed (taking into account that they have a finite transverse size).Since this overlap changes with rapidity, one obtains a non-trivial string density 3D distribution.
In this study, we follow the consideration from Ref. [23], where the colour field of fused gluonic tubes changes, which impacts particle production.This can also be interpreted as the appearance of string clusters with higher effective tension.
Apparently, there are several options for addressing the overlap and fusion of strings [42].In our previous study [36], we employed a cellular approach that involves the fusion of strings that have the centres in the same cell in the transverse plane.In this case, the cell size is equivalent to the diameter of a string.Thus, a cluster of k overlapped strings in Ref. [36] was replaced by one string with a decreased mean multiplicity per rapidity unit, an increased mean transverse momentum of produced particles and an increased probability of producing heavier particles.The transverse position of the cluster of k fused strings can be assumed to be the mean value of their k transverse positions.
Therefore, if one considers AA collision in the approach [36] discussed just above, one can expect that after the transverse dynamics (Section 2.4) the centres of all strings appear in the same cell.This leads to the replacement of all the variety of the degrees of strings' overlaps by a single string with enhanced tension.In this scenario, it is impossible to create the anisotropy of produced particles since the information about matter density fluctuations is lost.Let us note that in the models without attractive string dynamics in the transverse plane, this is not an issue since an event picture is less dense.Thus, in the current study, we do not consider simplifying string fusion on a coarse lattice but work using fine division of the transverse plane into small enough pixels.
In the approach considered here, the transverse plane is tessellated into square bins (pixels), each of which counts the occupation numbers of strings (number of strings that cover this pixel).In order to have a fine grid, one has to select the area of the bin, S bin , so that the latter is much smaller than the string transverse area, S 0 .In the current implementation, we use a string diameter, d str , of 0.5 fm and a pixel width, d bin , of 0.05 fm.
The resembling procedure is done for the rapidity space with a slice size of 0.1.In the previous study [36], we applied a longitudinal grid separately to each string.In this paper, we use a uniform splitting for the entire longitudinal dimension.
Finally, one can calculate the number of strings that occupy each 3D bin in the mixed configuration-momentum space.Thus, in a sense, we move away from the concept of particle-producing strings towards the concept of particle-producing pixels.
Note that the program's running time is significantly impacted by the number of 3D bins, as it must iterate over all 3D cells.Hence, it is necessary to restrict the volume permitted for simulations.Actually, the collision energy determines the longitudinal (in rapidity) size of the defined space, as the beam rapidity rises with √ s.In turn, the required transverse area depends on the initial distribution of strings and their final positions after the transverse dynamics.Since, in this investigation, we focus on the evolution of the system until the global minimum of its potential energy, after the motion according to Equation (5), the system becomes even more compact.

Fusion and the Kinetic Energy of Strings
Since the total energy of the system of strings in an event to be conserved, we assume that when strings overlap, a redistribution of their potential, U, and kinetic, T, energies occur.It is because the overlap of colour strings modifies the strength of the colour field inside the strings, which affects the string tension.Therefore, the partial overlap of a few strings leads to a decrease in the potential energy of each of them in this region.This modification causes an opposite change in the kinetic energy of strings, pulling them towards each other similar to that considered in Ref. [32].The formation of a string cluster, as a consequence of the strings attraction increases the string tension in comparison to a free string.Interestingly enough, the alternative option, also mentioned in Ref. [32], includes the decrease in the effective string tension and, therefore, should lead to the string repulsion (see also Ref. [33]).
As a result, one needs to parameterize the gain of kinetic energy, ∆T i,j , that the i-th string in an event obtains from the overlap with the j-th string.The functional form of ∆T i,j should reflect the feature that the full overlap of two strings means that their fusion is completed and those strings stop.On the other hand, ∆T i,j should decrease with the distance between the strings' centres.Thus, the following dependence on d i,j -the 2D distance between the i-th and j-th strings-is proposed: Here, d i,j ≤ 2r 0 , where r 0 = 0.25 fm is the string radius, while χ is a dimensional constant measured in GeV/fm that is a free model parameter.Hence, to get a new value of the kinetic energy of a string we iterate over its neighbouring overlapping strings in every rapidity slice.
We assume that the string motion according to Equation ( 5) is non-relativistic; therefore, in the following, we neglect the initial momentum of a string.As a result, the i-th string gains the extra transverse momentum from the j-th string: where m i str = m part1 + m part2 is the sum of the dynamically obtained masses of partons forming the i-th string (see Appendix A.3).
The azimuthal direction of the vector from the centre of the i-th string to the centre of the j-th string, φ i,j , is determined by the distances on X and Y axes between the centres of the strings, ∆X i,j and ∆Y i,j .The projections ∆p i X and ∆p i Y of the total momentum that the i-th string gains from all its overlapping neighbours, ∆p i T , are found by summing the projections of ∆p i,j T as T sin(φ i,j ).
Thus, the 2D vector of the string's transverse momentum, − → ∆p i T , is determined by its overlap with other strings in each rapidity slice.It is this extra momentum that is transmitted to particles produced by the string (see Section 3.1).

Particle Production
In this study, we define string hadronisation in each 3D bin.The strength of the colour field inside it determines the average number of charged particles produced and their average transverse momentum.Particles' rapidities are found from the uniform distribution in each rapidity slice.

Mean Multiplicity with String Fusion
Let us consider a rapidity slice of a free single string.The colour field, E 0 , inside it can be represented as a sum of colour fields inside all transverse (2D) bins that tessellate the area of this string, Thus, the field inside a 2D bin, E bin , is proportional to the ratio of its area, S bin , to the area of a string, S 0 .
Let us now consider a random 3D bin in the space that is populated with k strings.The resulting colour field, E tot , inside the string, is where E i is the field of a string in the bin.
For simplicity, the present study considers all strings to have the same transverse area, S i 0 , which does not vary along the rapidity direction, albeit the study can be more complex [54].The colour field of a free string, E i 0 , is assumed to be constant and is defined by the confinement properties [55].Therefore, in the case considered here, Equation ( 15) can be simplified to The average multiplicity from a free string, ν 0 , is proportional to the field of a string, E 0 .Thus, by transitivity, for the k strings that overlap in a 3D bin, one obtains ν tot ∝ E tot .
One can define the average multiplicity from a 3D bin with the length ε rap and transverse area S bin as with a free parameter µ 0 defining the average multiplicity of a unit of rapidity of a free string and S 0 = πr 2 0 .For comparison, in our previous study [36], we calculated the mean multiplicity from the cluster of k strings (without division in the transverse plane), µ clust , as µ clust = µ 0 ε rap √ k, following Ref.[40].
The actual multiplicity per 3D bin, N bin , is sampled from the Poisson distribution with a given mean, µ bin (17).Then, N ch is a sum of the multiplicities of N bin from all the 3D bins.

Mean Transverse Momentum with String Fusion
The mean transverse momentum of the particles produced by a free string, p 0 , does not depend on the string's length in rapidity or its transverse area.Thus, p 0 remains unchanged despite the division of a strings into pixels.Therefore, we keep the modification of p 0 for the cluster of k fused strings, p T k , from Ref. [36].Namely, where 19 ] is found in Refs.[49,56] from the fits to data and p 0 is a free model parameter.
Particles' transverse momentum is sampled from distribution corresponding to Schwinger mechanism of quark-antiquark pair creation in the colour field of a string that leads to its break-up and the formation of final hadrons [57][58][59], with the mean transverse momentum, p T k (18), from the cluster.

Probabilities of Producing Different Particle Species
The particle species with masses m sp are determined by Schwinger-like probabilities with the modified effective string tension, σ eff = 4p 2 0 k 2β , according to ∝ exp(−πm 2 sp /σ eff ), which is consistent with Equation (19).Typically, in the models that rely on the Schwinger mechanism of particle production [60], σ eff slightly differs from σ T (see Equation ( 7)) to effectively take into account the mechanism of particles re-scatterings.
Pions, kaons, protons and ρ-resonances are included in the model, with the ρ-resonance decaying into two charged pions.With particle masses, the longitudinal component of the particle momentum, p z , is calculated and, thus, the pseudorapidity is where | p| = p 2 T + p 2 z is an absolute value of the particle momentum.

The Origin of Particle Flow in the Model
A purpose of introducing the fine grid in the transverse plane as described in Section 2.6 above is to account for particle azimuthal flow in our model.We follow the consideration in Refs.[29,32] on the way of introducing collective behaviour in the colour string model without a hydro phase.There are three main ingredients.
First, in this type of models, the strings in an event form some clusters [23] distributed anisotropically.Alone, these fluctuations in string density do not produce whatever flow.However, this initial state anisotropy is crucial for the two mechanisms described below as controlling the strength of the effect of those mechanisms.
Second, we consider the change in the strings' kinetic energy that occurs when the strings overlap in 3D space.As a result, the boost from a string is transferred to the particles that it is fragmented to [32].The boost creates correlations in p T and φ between particles produced from strings that interact with each other.What is essential is that string hadronisation is considered at the moment of the minimum of the potential energy of the system of strings (at τ deepest ).This results in the high similarity of the geometry of events with the close number of strings, which leads to the same picture of particle boosts.For instance, Figure 2 (second top, second left) shows a typical event with 16 strings.We argue that all events with this n str resemble each other up to the event rotation and some statistical fluctuations.It means that the events have a non-zero flow of comparable value; therefore, the signal survives after averaging over such group of events.
Third, the particles passing through single strings or string clusters lose some part of the energy due to gluon radiation in the medium [29].When a particle loses its entire momentum, it means that this particle failed to escape.Thus, the azimuthal directions that are forbidden appear in the event.Consequently, the particles can no longer move in a whatever direction and φs of those particles become narrowed and correlated.Moreover, the complicated geometry of the string medium leads to different path lengths in different azimuthal directions and to different p T losses.
The first mechanism of those listed just above is naturally introduced into our model: transverse dynamics of strings result in the formation of string clusters of different degrees of overlap.Longitudinal dynamics make the fluctuations of string density dependent on the rapidity coordinate.
The other two mechanisms are the new features implemented in the model and are described in Sections 3.1 and 3.2 just below.Let us note that different degrees of string overlaps cause variations in the magnitude of energy loss and particle boosts bin-by-bin in the transverse plane-rapidity space.As a result, the particle production in a given event becomes highly asymmetric in azimuthal angle and pseudorapidity.

Transverse Momentum Boost
Particle's transverse momentum sampled from Equation ( 19) receives a Lorenz boost if the particle originates from the string piece that was accelerated due to string fusion [32].
The procedure iterates over all 3D bins and finds the strings that cover each bin.A bin is assigned the momentum projections, ∆p bin,i X and ∆p bin,i Y , that are found as fractions of the momentum projections of the i-th string, ∆p i X and ∆p i Y , defined in Equation (13) as where d bin = 0.05 fm is the bin size in X and Y in the transverse plane, δy i = |y part1 fin − y part2 fin | is the length of i-th string in rapidity which is calculated as the difference between the rapidities of the ends of the string, y part1 fin and y part2 fin , and highly fluctuates.Second, we find the part m bin,i of the mass of the i-th string, m i str , which belongs to this bin in a similar way Using Equations ( 21) and ( 22), we can find the energy of i-th string contained in this 3D bin, ∆E bin,i , as If this 3D bin is covered by k strings in an event, their ∆p bin,i X,Y from Equation ( 21) and ∆E bin,i from Equation ( 23) should be included in its total momentum, p bin,total X,Y , and the total energy inside it, E bin,total , as We perform particle production from each 3D bin in the bin rest frame as described in Section 2.8.In this frame, the particle's azimuthal angle, φ, is sampled from the uniform distribution from −π to π.To obtain the laboratory reference frame, one boosts the four-momenta of the produced particles using the boost vector with p bin,total X , p bin,total Y and E bin,total .
Thus, we obtain the correlated transverse motion of particles produced by a 3D bin.

Transverse Momentum Loss
When a particle traverses a 2D bin with a certain density of strings (occupation number), it loses part of its initial momentum, p init , due to gluon radiation, reaching the value, p fin .In an analogy with the QED process of photon radiation by charged particles moving in the external electromagnetic field, it can be expressed in the following way [29] where κ is a quenching coefficient that is a free model parameter.It is necessary to apply the Equation ( 26) iteratively since σ eff = 4p 2 0 k 2β varies bin-by-bin based on the number of overlapped strings, k.That is why l, a 2D path that a particle accomplishes, can be approximated by d bin √ 2 (the length of the 2D bin's diagonal) for each step.Since the density of strings fluctuates with rapidity (see Sections 2.4 and 2.5), l differs not only at different azimuthal angles, φ, but also within different rapidity slices, ε rap .Let us remark that, in general, the value of p fin can become negative after a number of iterations of Equation (26).We interpret this as the impossibility of the particle to escape in a given azimuthal direction; thus, it is absorbed by the string environment.In this case, the energy of the particle is spent on producing another particle to replace the former one.For a new particle, the model regenerates the transverse momentum and azimuthal angle of this particle; thus, the particle gets a chance to escape in a new direction.As soon as the required number of particles to be produced from this 3D bin (from Poisson distribution with the mean from Equation ( 17)) is known, new particles are to be regenerated until the particle leaves the string medium.
To find the trajectory of a particle in the transverse plane, the Bresenham's line algorithm [61] is applied for each rapidity slice.

Inference of the Model Parameters
We determine the model parameters by comparing the model distributions with the ATLAS experiment data on inelastic proton-proton interactions at √ s = 13 TeV [45] (Figure 4).We describe not only the global observables such as N ch -distribution, p T and η spectra, but also the p T -N ch correlation function.The dN ch /dη| η≈0 region (Figure 4, upper left) is adjusted by finding the appropriate value of the mean particle multiplicity µ 0 (see Equation ( 17)) from a 3D bin of a free string.In turn, the width of the dN ch /dη distribution is controlled by the value of the string tension, σ T (see Equation ( 7)).
By tuning p 0 (see Equation ( 18)), we settle the p T at low N ch (Figure 4, lower right).With the proper selection of κ (see Equation ( 26)), we are able to modify the slope of the p T -N ch correlation function at moderate and high N ch .The p T is influenced by the value of χ (see Equation (11)) at high N ch , when string fusion is most significant.
In the charged particle multiplicity distribution (Figure 4, lower left) the events with N ch < 2 are removed as in the data; thus, we plot P(N ch ) = P full (N ch )/(1 − P full (0) − P full (1)).
Here, P full denotes the charged particle multiplicity distribution that includes events without registered particles (N ch = 0).These events are highly influenced by diffractive processes that are not considered in the model.Moreover, even with the N ch ≥ 2 selection the experimental results at low N ch and low p T are affected.Therefore, our predictions could not be directly compared in this region, but we are mostly interested in the high N ch events as those events are more relevant in the studies of collective behaviour.The resulting p T -spectrum is presented in Figure 4 (upper right).

Flow Measurements and the Quantities of Interest
One can estimate the flow signal by performing the Fourier expansion of the singleparticle distribution in the azimuthal angle, φ, [7].It is necessary to take into account the intrinsic symmetry of particle production in every event, the latter being determined by its reaction plane, Ψ RP .The plane is formed by the direction of the beam and the impact parameter and creates a preferred azimuthal orientation in an event.Thus, the relative particle azimuthal angles, φ-Ψ RP , should be used in the calculations instead of the measured φ.Therefore, in a given part of phase-space, one can expand the distribution into a Fourier series as where n is the Fourier harmonic number and v n is the corresponding Fourier coefficient.The full set of v n s describes the amplitudes of particle distribution asymmetry in the transverse plane averaged over particles in one event.Let us note that both v n and Ψ RP fluctuate event-by-event and, typically, only the moments of the corresponding distributions are measured experimentally.The validity of the Fourier expansion in the real case of finite event multiplicity (especially in pp collisions) is questionable.Moreover, the reaction plane, Ψ RP , cannot be directly measured, so one may substitute it by proxy [62] called the event plane, Ψ EP .However, there is no unique event plane in an event; instead, one determines a set of event planes, Ψ n , depending on n as where the numerator and denominator are calculated from the distributions of the particles in φ per event.
More recent studies have shown [63] that imprecise estimation of the event plane significantly spoils the flow signal in this approach.Therefore, more sophisticated measures to be used such as two-particle correlations that, under certain assumptions, naturally exclude the dependence on To plot the results on cumulants (Figures 6 and 7) in a unified way for different event classifications, we correlate N sel ch with the charged-particle multiplicity N ch calculated on particles with |η| < 2.5 and p T > 0.4 GeV.We find the event-mean N ch for every event class based on N sel ch and plot it along the X-axis.Figure 6 shows the model results for c 2 {2} (Figure 6, upper) and v 2 {2} (Figure 6, lower) calculated for particles with |η| < 2.5 and 0.3 < p T < 3.0 GeV, as a function of N ch .Results are presented for event classes with a 0.5% width of the N sel ch distribution calculated for different event classifications based on p T , as indicated.
The values of c 2 {2} are positive and show an increasing trend towards high-multiplicity collisions.The associated second flow harmonic, v 2 {2}, repeats this behaviour.We interpret it as follows.The string fusion plays a more significant role with the increased string density, therefore, both quenching and boosting of particles become stronger, leading to a greater flow signal.
A similar dependence was presented by ALICE experiment (see [69]) for the twoparticle correlation, Υ(∆φ), integrated over rapidity, that grows with the event multiplicity.The results were also presented for the near-side ridge, but those differ in particle selection and acceptance.
For large N ch in Figures 6 and 7, one can see a slight splitting of the results with different N sel ch definitions.The primary goal of the analysis with the different selections is to test the sensitivity of the flow to particles' transverse momenta.To make it more visible, we show in Figure 8 the p T dependence of two-particle cumulants, c 2 {2} and c 2 {2} (2subs) , for 0-5% event class of N sel ch distribution for particles with 0.3 < p T < 3.0 GeV and |η| < 2.5.To achieve a similar number of particles for all model points, we perform the calculation of cumulants in the p T -intervals of a varying size.Figure 8 shows that c 2 {2} and c 2 {2} (2subs) values gradually increase as the functions of the momentum of particles that are used in their calculations.This is in qualitative agreement with the ATLAS results [66] showing a similar trend.
In the model workflow, the p T behaviour of c 2 {2} and c 2 {2} (2subs) can be explained as follows.First of all, the quenching mechanism that was introduced produces larger anisotropy for particles with higher p T .It is because it becomes increasingly more difficult for particles to escape the string matter while saving sufficiently large transverse momentum.In turn, string-string interactions and their acceleration in the transverse plane increase the p T of correlated particles, leading to larger anisotropy from these boosts.

Discussion
In this study, we developed the workflow of the colour string model to address the challenging question of the origin of the collective flow observed in inelastic proton-proton interactions.It is interesting to draw parallels between our simplified model and the cognate event generators also based on colour string fragmentation such as in the EPOS [25] and PYTHIA [26] models.
The core concept of our model is the description of pp inelastic interaction via the multi-pomeron exchange, Equation (1), find which resembles the parallel scatterings of partons that occur in the EPOS4 event generator.The way to find momenta of string endpoints combines methods of energy-momentum sharing and saturation implemented in the EPOS4 model.All this is strikingly different from the consideration of hard scattering a starting point of multi-parton interactions in the PYTHIA model.The concept that was not implemented in any existing event generator is attractive string dynamics [39] introduced to the workflow of our model.It modifies the shape of the string density in an event resulting in the collapse of particle-producing sources.This leads to the fusion of overlapped strings into an event hot spot, which is analogues to the formation of core in the EPOS4 model.On the other hand, the concept of the formation of the fused particle source with higher tension is realised in the PYTHIA generator by using the rope mechanism.Finally, in the EPOS4 generator, evolution of the core is considered in hydro-regime, which results in a transverse flow signal.In the presented model, particle boosts (similar to the effect of string shoving mechanism in PYTHIA) and their momentum losses in string medium are the sources of azimuthal anisotropy of produced particles.
Moving on to discussing the obtained results, let us note that an event picture prior to hadronisation looks like a highly inhomogeneous string medium.Originating from the multi-pomeron exchange, the system of particle-producing strings is disordered by the longitudinal and transverse dynamics of the system.One observes the peculiar grouping of strings in the mixed configuration-momentum space.Some of them are isolated forming the debris of pp interaction.Other strings participate in the event hot spot, partially overlapping, while the most lucky ones form a dense cluster of strings, overlapping to the highest degree.It is the presence of such a core that determines the crucial collective features of the soft particle production in pp collisions.In particular, in the model, the particle anisotropy appears due to this complex structure of event shape and the effect of particle momentum quenching.
Thus, a recognisable core-corona picture of the string system [70] is observed.This picture represents the interplay of these two characteristic event regions that governs the transverse flow in the model.Namely, in the corona region, it is only the loss of the particle's momentum that plays a role.On the other hand, the impact of the complex core structure is quite complicated to deduce.The relatively low occupancy regions of the core create correlated multi-directional particle boosts, while the hottest region determines a single dominant direction and, therefore, strongly boosts particles with ∆φ ≈ 0.
For example, the model result for the angular correlation function, C(∆η, ∆φ), reveals the familiar near-side ridge that is formed by the collimated emission of particles, boosted by the clusters of strings.
On the other hand, the away-side ridge is missing from the model results.We interpret that as due to an excessive approach of strings at the τ deepest moment, which results in the formation of the over-condensed string region.
In this case, the production of particles with ∆φ ≈ π separation in the model framework with p T boost should come from the peripheral, low-occupancy parts of the core.It seems that, in the current implementation, this area is too scarce to create the away-side ridge.A similar observation was made in Ref. [71], where the absence of the away-side ridge is caused by complications in defining the cutoff between the core and corona.
The model results obtained for the second-order two-particle cumulants are in qualitative agreement with the data by the ALICE and ATLAS experiments on inelastic pp interactions at √ s = 13 TeV.In our future studies, higher statistics would be needed to calculate multi-particle cumulants and to access high-multiplicity events.The flow signal obtained in the model grows with the transverse momentum of charged particles in high-multiplicity events.This reflects the flow origin from the particles formed by interacting strings.
In conclusion, let us note that in this study, it was expected that the τ deepest time-the time of the string system evolution in the transverse plane-that results in the highest density of strings would provide the best description of the collective behaviour of particles in pp interactions.However, it appeared that, in terms of this model, τ deepest controls a degree of core-corona separation.The obtained ∆η-∆φ correlation result indicates that one would have to tune τ deepest more precisely, which is a subject for further studies.and the initial parton energy as The modifications of partons' momenta (A3) and energies (A4) change the partons' masses from "bare" (current) values, m di or m u/d/s/c , to the "dressed" ones, m dressed q/di , according to This approach can be naturally extrapolated to any arbitrary number of sea quarks.Namely, we assign x 0 /3 and e 0 /3 to each of the valence quark and the diquark and split the remaining 1/3 between all the sea quarks in a proton.This procedure increases the masses of quarks and diquark (for valence ones to a greater extent), which makes the distributions of the string ends' rapidities more realistic.However, it is also possible that after the procedure described, parton's energy decreases compared to parton's momentum.Therefore, one cannot calculate parton's dressed mass, m dressed q/di .In this case, the remaining two different groupings of valence quarks into diquark should be considered (see Appendix A.1).If none of the combinations allows finding m dressed q/di , the proton is to be regenerated.

Figure 1 .
Figure 1.Model event distribution, P(n str ), of the number of strings for inelastic pp interactions at the centre-of mass energy √ s = 13 TeV.

Figure 2 .
Figure 2.The simulated events (one event a row) with 2, 16, 38, and 74 strings (top to bottom) for an event projection to the transverse plane, X-Y, before and after the transverse evolution and to the X-rapidity plane before and after the transverse evolution (left to right) for the time of the transverse evolution, t transv , as indicated.The evolution runs till t transv = t deepest (second left and most right).Z axis is not to scale.See text for more details.

Figure 2 .
Figure 2.The simulated events (one event a row) with 2, 16, 38, and 74 strings (top to bottom) for an event projection to the transverse plane, X-Y, before and after the transverse evolution and to the X-rapidity plane before and after the transverse evolution (left to right) for the time of the transverse evolution, τ transv , as indicated.The evolution runs till τ transv = τ deepest (second left and most right).Z axis is not to scale.See text for more details.

Figure 3 .
Figure 3.The simulated events (one event a row) with 2, 16, 38, and 74 strings (top to bottom) for an event projection to the X-rapidity plane before and after the longitudinal evolution and to the Y-rapidity plane before and after the longitudinal evolution (left to right) for the time of the longitudinal evolution, t long , as indicated.The evolution runs till t long = t deepest (second left and most right).Z is axis not to scale.

Figure 8 .
Figure 8. Model results for two-particle cumulants, c 2 {2} and c 2 {2} (2subs) , as a function of p T of charged particles for 0-5% event class based on N sel ch calculated for 0.3 < p T < 3.0 GeV, |η| < 2.5 in inelastic pp interactions at √ s = 13 TeV.See text for details.