Setting the string shoving picture in a new frame

Based on the recent success of the \angantyr model in describing multiplicity distributions of the hadronic final state in high energy heavy ion collisions, we investigate how far one can go with a such a string-based scenario to describe also flow effects measured in such collisions. For this purpose we improve our previous so-called \textit{shoving} model, where strings that are close in space--time tend to repel each other in a way that could generate anisotropic flow, and we find that this model can indeed generate such flows in \AA\ collisions. The flow generated is not quite enough to reproduce measurements, but we identify some short-comings in the presented implementation of the model that, when fixed, could plausibly give a more realistic amount of flow.


Introduction
High energy proton-proton and heavy ion (HI) collisions are frequently analysed assuming very different dynamical mechanisms. Models based on string formation or cluster chains and subsequent hadronization, implemented in event generators like PYTHIA [1,2], HER-WIG [3,4], or SHERPA [5,6], have been very successful in describing particle production in e + e − annihilation and pp collisions. For pp collisions this picture includes multi-parton subcollisions, and at high energies the scattered partons are mainly gluons, which form colour connected gluon chains. This picture corresponds to multiple pomeron exchange (including pomeron vertices), with a gluon chain pictured as a cut BFKL pomeron ladder. It is also consistent with models based on reggeon theory, e.g. the models by the groups in Tel Aviv [7], in Durham [8], or by Ostapchenko [9]. The gluons here also hadronize forming strings or cluster chains.
On the other hand, many features in HI collisions can be described assuming the formation of a thermalized quark-gluon plasma (QGP), in particular collective flow [10][11][12][13][14][15] and enhancement of strange particles [16]. The formation of a hot thermalized plasma is often described as follows: At high energy the gluons in a nucleus form a "color glass condensate" (CGC) [17,18]. This state is described by a classical colour field, where the strength is saturated due to unitarity, with a "saturation scale" Q s ∼ (1/α s )(x/x 0 ) −λ . When two nuclei collide the overlapping non-Abelian fields form a "glasma"; for reviews of the CGC and the glasma see e.g. [19][20][21]. The glasma state contains parallel longitudinal colour-electric and colour-magnetic fields, associated with induced magnetic charges in the projectile and target remnants. These fields also build up a topological charge giving rise to CP-violating effects [22,23]. The glasma is instable, and in this scenario it turns rapidly into a QGP, which soon thermalizes [24]. This transition is often motivated referring to "Nielsen-Olesen instabilities" in the QCD Fock vacuum [25,26], or "Weibel instabilities" in an electro-magnetic plasma [27].
Lately the difference between pp and nuclear collisions has become much less clear. Both collective effects and increased strangeness production have been observed in high multiplicity pp events (see e.g. [28] and [29]). This has raised the question whether a plasma can be formed in pp collisions, or alternatively, if these effects in HI collisions also can be described in a scenario based on strings.
In high multiplicity pp events, the density of strings is quite high. The width of a stringlike flux tube is estimated to be of the order of 1/Λ QCD , which implies that the flux tubes will overlap in space. Recently we have demonstrated, that in pp collisions both collectivity and strangeness enhancement can be explained as consequences of a higher energy density in systems of overlapping strings or "ropes". This gives both a transverse pressure, which can cause long range collective flow [30], and enhanced strangeness following from the higher energy release in the breakup of a rope [31]. At the same time string-based models have successfully described the general particle distributions in nuclear collisions [32,33]. In this paper we will show that some of the collective effects seen in AA collisions can be described by the transverse pressure from overlapping strings, in the same way as in pp collisions, without assuming the formation of a thermalized state with high temperature.
The strangeness enhancement in AA expected when nearby strings form "ropes", will be discussed in a later publication. We expect this effect to be similar, but also further enhanced in AA compared to pp, which was described in ref. [31].
The two scenarios discussed here are qualitatively different. It is important to note that the difference lies not in the initial gluonic states in the projectile and the target. In the pp MC models mentioned above, the parton distributions are fitted to HERA data, and thus include the observed suppression for Q 2 < Q 2 s (x), as a result of saturation at small x. They also agree with BFKL evolution including saturation effects, as formulated in the dipole model DIPSY [34,35], which has the same evolution with rapidity as the CGC, including the suppression for q ⊥ < Q s (x), as an effect of Debye screening at high gluon density. In the limit N c → ∞ this follows because the gluon emissions split the dipoles into gradually smaller and smaller dipoles. For finite N c the screening effects are further enhanced by identical colours in neighbouring dipoles, described by colour reconnection in the initial state cascade via the "swing mechanism" in the DIPSY model [34].
The difference lies instead in the assumed behaviour just after the collision. The particle distribution is (in symmetric collisions) characterized by an approximately boost invariant central plateau. This feature is quite natural if the hadrons are produced from boost invariant strings stretching between the projectile and target remnants. The energy per unit length in a string, dE/dz, is given by the string tension κ ≈ 0.9 GeV/fm. The energy in the collision is initially concentrated in the forward-going remnants. These gradually loose energy, which is stored in the new string pieces, as the strings become longer. The energy density, dE/d 3 x, is therefore approximately constant from the instant of collision until the hadronization (proper) time τ H ∼ 1.5 − 2 fm (see section 3). From this moment the hadrons move out with rapidities approximately equal to the hyperbolic angle η st ≡ (1/2) ln[(t + z)/(t − z)] corresponding to the place where they were "born". 1 In the center the relation between η st and z is given by ∆z = t ∆η st . (1.1) For a single string the final state energy density in rapidity is therefore given by (replacing t by the boost invariant variable τ = √ t 2 − z 2 ) A boost invariant plateau is also expected from a glasma, if the longitudinal colourelectric and -magnetic fields are stretched out in a similar way. But the glasma is expected to rapidly transform into a thermalized plasma, and this plasma will expand longitudinally and dilute with an energy density falling like τ −1 or faster, before freezeout or hadronization time [36]. The energy density in the glasma must therefore correspond to the (energy density at freezeout)×(τ H /τ T ), where τ T is the time for the glasma-plasma transition.
The different energy densities also implies that, in the string picture, the features of the vacuum condensate is very important to keep the field together in the strings, while the condensate energy is too small to play an important role in the plasma scenario, where it is normally neglected; see e.g. ref. [24] where this is explicitly stated. The features of the energy density and the vacuum condensate are discussed in somewhat more detail in section 2.
If both approaches are able to reproduce presently available data on, in particular, collective effects and strangeness enhancement, it becomes important to find observables where they will give different predictions. As an example, one such observable would be effects of CP violations, where clear signals are not expected in the string scenario.
We have already presented some preliminary studies in ref. [30] where we used a simplified so-called shoving model to describe how strings in a pp collision repel each other to create the so-called near-side ridge first found in pp by CMS at the LHC [28]. The simplified model had a number of short-comings. One was that it only treated strings that were parallel to the beam direction, using an upper cut on the transverse momentum of the partons in the string. Another was that the force between the strings was manifested in terms of many very soft gluons, which was technically difficult to handle in the PYTHIA8 hadronization implementation, and also somewhat difficult to reconcile within the string model as such. 2 In this work we will present a more advanced model, where all string pieces in an event can interact, not only in pp, but also in AA and basically in any other collision system. Furthermore the problem with producing too many additional soft gluons is circumvented by applying the force in terms of tiny transverse momentum nudges given directly to the produced hadrons in the hadronization. Even though the new model does not produce extra soft gluons, it still has some problem dealing with soft gluons already present, stemming from the perturbative phase in PYTHIA8, that complicates the string motion, and currently this restricts the amount of shoving that can be achieved.
The layout of this paper is as follows. In section 2 we first further elaborate the motivation behind trying to treat high density AA collisions in terms of interacting strings, rather than as a hot thermalized plasma of free partons. There we also discuss in some detail what is known about the structure of a QCD string using the analogy with a vortex lines in a superconductor. In section 3 we then look at how the QCD string is described in the Lund string fragmentation model. Our new shoving model is presented in section 4 giving some details of how one finds a Lorentz frame where any two string pieces always lie in parallel planes, and how we there can discretize the shoving into tiny transverse nudges between them which are then applied to the hadrons produced. We study the behaviour of the new model in sections 5 and 6. First we apply it to a toy model for the geometry of the initial partonic state of an AA collision to show that we can qualitatively describe some of the flow effects found there. Then we apply the new model to a more realistic initial state generated by the Angantyr model [33] for AA collisions in PYTHIA8, and find that the quantitative description if flow effects AA collisions is still lacking. Finally in section 7 we present our conclusion and present some ideas for improvement of our new model that may achieve an improved description of experimental measurements.

Initial energy density and possible plasma transition
As mentioned in the introduction, one of the most popular scenarios for treating heavy ion collisions, involves early time dynamics calculated as quasi-classical gluon fields immediately after the collision (a "glasma"), followed by a thermalized plasma phase treated with (3+1D) hydrodynamics [11]. The initial glasma state contains longitudinal colour-electric and colour-magnetic fields associated with induced magnetic charges in the forward-moving nucleus remnants. Longitudinal electric or magnetic fields are boost invariant, have no longitudinal momentum and a constant energy per unit length. When the longitudinal glasma fields become extended, transverse fields are generted, which keep the flux contained in tubes. In ref. [24] it is shown that the energy density is proportional to ( where Q is a typical scale related to the saturation scale Q s . Such nearly constant initial fields have essential similarities with the relativistic strings in the Lund model. The glasma fields are, however, assumed to rapidly turn into a thermalized QGP, which expands longitudinally, with an energy density falling like τ −1 (or faster) [36]. To give the observed energy density after freezeout, the energy density before thermalization must be correspondingly higher.
This picture can be contrasted to a string scenario, where the longitudinal energy density dE/dz is constant until hadronization. As discussed in the introduction, if we look in the rest frame of the hadrons produced by a string piece with hyperbolic angle ∆η st , then at proper time τ this piece corresponds to a distance ∆z = τ ∆η st in the normal coordinate space.
When the string hadronizes the average density in rapidity y is the same as in η st , giving for a single string an energy density dE/dy ≈ dE/dη st = κ · τ H with κ ≈ 0.9 GeV/fm and the hadronization time τ H = 1.5 − 2 fm (see section 3). Minijets will enhance the energy density dE/dy in the lab frame, but also here the initial density is limited.
The energy density and the time for thermalization in the glasma state are difficult to estimate theoretically. They are frequently present in lattice units or arbitrary units (e.g. in ref. [24]) or they can be adjusted to fit experimental data as e.g. in ref. [38]. Tuning to a hydrodynamical plasma calculation for RHIC data at 200 GeV and b = 9 fm (centrality ≈ 40 %) [38], finds an average density dE/τ dy ≈ dE/dz ≈ 300 GeV/fm. The experimental result for this centrality is dN ch /dy ≈ 150 [39,40]. Adding neutrals and assuming an average transverse mass ∼ 0.5 GeV, this gives dE/dy ≈ 110 GeV, which agrees with eq. (1.1) if the thermalization time is about 0.35 fm. The energy density would then be dE/d 3 x ∼ 8 GeV/fm 3 . A corresponding calculation for the string scenario, where the strings hadronize at τ = 1.5 − 2 fm, would instead give the density 1.7 GeV/fm 3 .
In a central PbPb collision at LHC energy, the charged particle density dN ch /dy is of the order 2000 (c.f. the data from ALICE shown in figure 16). Including neutrals gives approximately 3000, and again assuming an average transverse mass ∼ 0.5 GeV and a transverse area ∼ 150 fm 2 , this corresponds to an energy density dE/dy d 2 x ⊥ ∼ 10 GeV/fm 2 . In the string scenario, eq. (1.1) gives an initial energy density in coordinate space dE/d 3 x ≈ 5 − 6 GeV/fm 3 . The energy density in the glasma is also calculated for LHC energy by Schenke et al. in ref. [41]. The energy density dE/τ dy ≈ dE/dz, obtained at τ = 0.1 fm, is about 20,000 GeV/fm, or dE/d 3 x ≈ 500 GeV/fm 3 . At this time the transverse and longitudinal fields have become approximately equally strong, and this density also roughly matches the density needed to reproduce the experimental data from ALICE, for a thermalization at τ = 0.1 fm. We note that this energy density is almost a factor 100 larger than what is needed in the string scenario. (We also note that the energy density obtained in ref. [41] is not approximately constant for small times, but instead falls rapidly from approximately 50,000 GeV/fm at very early times to about 20,000 GeV/fm at τ = 0.1 fm.) The conclusion from these examples is that the initial energy density needed in the glasma is one or two orders of magnitude larger than in the string scenario. This also implies that, while the energy density in the vacuum condensate can safely neglected in the glasma, it is important in the string scenario.

The vacuum condensate and dual QCD
As discussed above, the initial energy density is relatively low in the non-thermal scenario, and the properties of the vacuum condensate is important for the formation of colourelectric flux tubes (strings). The QCD vacuum was early studied in several papers by N.K. Nielsen, H.B. Nielsen and P. Olesen. In ref. [25] it was shown that the QCD Fock vacuum (with zero fields) is instable. If a small homogeneous colour-magnetic field is added, it will grow exponentially. However, as shown in the accompanying paper [26], if the externally added field is in a given direction colour space, it will induce fields in the orthogonal directions in colour space. (For simplicity the authors studied SU(2) Yang-Mills theory.) Higher orders in this induced field then develop a Higgs potential, analogous to the potential describing the condensate in a superconductor. This implies that the exponential growth will be stopped, and the initial Fock vacuum will fall into a non-perturbative ground state with negative energy. 3 In a normal superconductor a magnetic field is compressed in vortex lines or flux tubes, and magnetic charges would be confined. The pure Yang-Mills theory is symmetric under exchange of electric and magnetic fields. Quarks with colour-electric charge are confined, and the Copenhagen vacuum is also known to have non-trivial vortex solutions of electric type (see e.g. refs. [42,43]). Thus the QCD vacuum behaves as a "dual superconductor", with exchanged roles for the electric and the magnetic fields. For a review of the "Copenhagen vacuum", see e.g. [44].
It was further demonstrated by 't Hooft that the ground state in a pure SU(3) Yang-Mills theory can have two different phases, with either colour-electric or colour-magnetic flux tubes, but not both [45]. In the first phase colour-electric charge is confined, and extended colour-magnetic strings are not possible, while the opposite is true in the second phase. The fact that quarks are confined obviously shows that the QCD vacuum is of the first kind.
Although the pure Yang-Mills theory is symmetric under exchange of electric and magnetic fields, Maxwell's equations for electromagnetism are asymmetric, due to the absence of magnetic charges. Exchanging electric and magnetic fields is obtained by replacing the field F µν by the dual field tensor F µν ≡ µνκλ F κλ . The electric and magnetic currents would then be given by Expressing the field F µν as derivatives of a vector potential A µ implies, however, that the magnetic current is identically zero. As a consequence magnetic charges can only be introduced by adding extra degrees of freedom. Dirac restored the symmetry between electric and magnetic charges by adding a "string term", an antisymmetric tensor G µν satisfying ∂ ν G νµ = j µ magn . The electromagnetic field is then given by A consistency constraint is then a quantization of electric (e) and magnetic (g) charges: e · g = 2πn with n an integer. The features of a vacuum behaving as a dual superconductor was early discussed in a series of papers by Baker, Ball, and Zachariasen; for a review see ref. [46]. In a dual superconductor a colour-electric flux tube has to be kept together by a colour-magnetic current. Instead of expressing the extra degrees of freedom in terms of Dirac's string term, Baker et al. treated F µν and F µν as independent fields. (In the non-trivial vacuum condensate they are related by a non-local magnetic permeability.) Higher order corrections then induce a Higgs potential in the F field, analogous to the induced field in the Nielsen-Olesen instability.
The fields F µν and F µν studied by Baker et al. are non-Abelian. A problem is here that, although a bound qq pair must be a colour singlet, and the energy density has to be gauge invariant, this is not the case for the colour-electric and -magnetic fields.'t Hooft has, however, shown that the essential confining features can be described by an Abelian subgroup U (1) 2 to SU(3) [47]. More recent studies of dual QCD are based on this "Abelian dominance". This is also supported by lattice calculations performed in the maximal Abelian gauge, which exhibit a confining phase related to the condensation of magnetic monopoles [48,49]. For overviews of Abelian projections (or Abelian gauge fixing) see also refs. [50,51].

A single QCD flux tube in equilibrium
A normal superconductor can be described by the Landau-Ginzberg (LG) equations, see e.g. ref. [52] and appendix A. Here the vacuum condensate is formed by Cooper pairs and influenced by a Higgs potential. In the interior of the superconductor a magnetic field is kept inside flux tubes or vortex lines by an electric current, and magnetic monopoles would be confined. The flux in such a vortex line (and the charge of a possible monopole) is quantized in multiples of 2π/q, where q = 2e is the charge of a Cooper pair.
A vortex line or a flux tube is characterized by two fundamental scales: the penetration depth λ and the coherence length ξ. These scales are the inverse of, respectively, the mass attained by the gauge boson and the mass of the Higgs particle. At the boundary between a superconducting and a normal phase, if λ ξ the parameter λ determines how far the magnetic field penetrates into the condensate (from which it is expelled by the Meissner effect). Similarly if ξ λ, then ξ determines the rate at which the condensate goes to zero at the boundary. When ξ > λ (or more exactly ξ > √ 2 λ), both the condensate and the field are suppressed over a range ξ − λ. This is a type I superconductor, and it implies that the surface provides a positive contribution to the energy. In equilibrium the surface then tends to be as small as possible. If in contrast λ is larger than ξ (type II superconductor), the condensate and the field can coexist over a range ∼ λ − ξ, and the surface provides a negative contribution to the energy, favouring a large surface. (The properties of a classical superconductor are discussed in somewhat more detail in appendix A.) As discussed above, the QCD vacuum has important similarities with a superconductor. There are, however, also important differences between a non-Abelian flux in QCD and an Abelian field in a non-relativistic superconductor. The infrared problems in QCD contribute to the difficulties to estimate the properties of a QCD flux tube. The total energy, given by the string tension, is fairly well determined from the spectrum of quarkonium bound states, and approximately equal to 1 GeV/fm. The total flux and also the width of the tube are, however, less well known. One problem is that it is not obvious how much of the string tension is in the field, how much is in destroying the condensate, and how much is in the current keeping the flux together inside the tube. Another problem is that, although the flux is given by the running strong coupling, the scale is not well specified. Also for a fixed total flux, the energy stored in the linear colour-electric field depends on the (not well known) width of the tube. The field energy is ∼ A E 2 ∼ Φ 2 /A, where Φ is the total flux and A the transverse area of the tube. We here briefly discuss three approaches to estimate the properties of a QCD flux tube: i) The bag model: The simplest model for a QCD flux tube is the MIT bag model [53]. Here the vacuum condensate is destroyed within a radius R around the center of the tube. Inside the tube there is a homogeneous longitudinal colour-electric field E = Φ/(πR 2 ), where Φ is the total flux. The energy per unit length of the tube is then where the bag constant B is (minus) the energy density in the condensate. Equilibrium is obtained by minimizing κ, which gives πR 2 = Φ 2 /2B and κ = 2πR 2 B. Here half of the energy is in the field, and half comes from destroying the vacuum condensate inside the tube. With κ ≈ 1 GeV/fm and B ≈ (145 MeV) 4 [53], we find R ≈ 1.7 fm, or ρ 2 ≈ 1.2 fm. (Here ρ is the radial distance in cylinder coordinates.) ii) Dual QCD: In the approach by Baker et al. the fields D and H are derivatives of a dual potential C µ . In the vacuum condensate the magnetic permeability µ is non-local, and the fields E = µD and B = µH are treated as independent fields (note that µ = 1) related to a tensor field F µν . This tensor field interacts via a Higgs potential, forming a vacuum condensate which confines the colour-electric field D. In ref. [46] the LG parameter κ LG = λ/ξ is estimated to ≈ 0.75, when fitting the model to the string tension and the energy in the vacuum condensate. This is quite close to the borderline, 1/ √ 2, between type I and type II superconductors. It implies that the surface energy is small, and the authors conclude that the behaviour is not far from a flux tube in the bag model. The width ρ 2 is estimated to 0.95 fm. A more recent review over dual QCD is found in ref. [54], which also discusses results from lattice calculations.
iii) Lattice QCD: Lattice calculations are natural tools to solve infrared problems, and several groups have presented studies of flux tubes using different methods; for some recent analyses see refs. [55][56][57][58][59][60][61]. Such analyses are, however, not without problems. One problem is that the strings have to be rather short, in order to have a good resolution. A good resolution is important for the determination of the behaviour for small ρ. Besides the colour-electric field, the features of the flux tube also depends on the properties of the vacuum condensate.
Some studies find that the vacuum acts like type II superconductor, e.g. refs. [56,61], but a majority of recent studies conclude instead that it should be of type I, e.g. refs. [55,57,59]. In most analyses the electric field is fitted using Clem's ansatz [62] for the condensate Here ρ is the radius in cylinder coordinates, and the parameter ξ v is varied to minimize the string energy. As discussed in more detail in appendix A, this ansatz satisfies one of LG's two equations. Minimizing the string tension then gives the electric field Here K 0 is a modified Bessel function, and the scale parameter λ equals the penetration depth in the LG equations. The parameter ξ v is tuned to fit the shape of the lattice result for small ρ-values, where it suppresses the logarithmic singularity in K 0 . The result is related, but not equal to the coherence length in the superconductor. An important problem is, however, that the ansatz in eq. (2.4) is expected to work well for type II superconductors, where it gives an approximate solution also to the second of LG's two equations. It does, however, not work for type I superconductors, where the second LG equation can be quite badly violated. This is a problem for those studies, which find fits which correspond to κ LG -values in the type I-region, where the fitted values of λ and ξ no longer represent the initial parameters in the LG equations. This problem is also discussed further in appendix A. Another problem is how to translate the width of the field from lattice units to the physical scale in fm. This problem is discussed in a review by Sommer in ref. [63]. In earlier studies it was common to adjust the energy in the colour-electric field, 2πρ dρE 2 /2, to the string tension, known to be ≈ 1 GeV/fm. Also, as noted above, it is uncertain how much of the string tension is due to the field, how much is due to breaking the condensate, and how much is due to the current which keeps the flux together. In the bag model or Clem profile Gaussian profile Lattice calculation Figure 1: Profile of the electric field from the lattice calculation in ref. [55] compared to the fit by Clem [62] and a Gaussian distribution. the dual QCD estimate mentioned above, about half of the string energy is in the field, and half is due to the condensate. Adjusting the field to represent the full string tension is therefore likely to overestimate the field strength and thus underestimate its radius.
Another way to determine the scale, used in the lattice analyses mentioned above, is to study the transition of the qq potential from small to larger separations r, e.g. using a fit of the type V (r) = A/r + Br + C, see e.g. refs. [64,65]. The parameters can here be adjusted e.g. to a phenomenological fit to quarkonium spectra. This method is, however, also uncertain. In the review by Sommer, mentioned above, it is concluded that "the connection of the phenomenological potentials to the static potential V (r) has never been truely quantitative".
Although there are uncertainties in the determination of the LG parameter using Clem's approximation, it does give a good description of the shape of the longitudinal field obtained in the lattice calculations. As an example we show in figure 1 the result in ref. [55]. We here note that the profile is also quite well represented by a Gaussian distribution. Due to the uncertainties in determining the value of κ LG and the physical scale, we will below approximate the field by a Gaussian determined by two tunable parameters, which are related to the radius of the field and to the fraction of the string tension related to the colour-electric field energy.

The Lund string hadronization
In the Lund string hadronization model, described in ref. [66], it is assumed that the dynamics of a single flux tube, and its breakup via quark pair production, is insensitive to the width of the tube. It is then approximated by an infinitely thin "massless relativistic string". Essential features of the Lund hadronization model are first qq pair production via the Schwinger mechanism, and secondly the interpretation of gluons as transverse excitations on the string. This last assumption implies that the hadronization model is infrared safe and insensitive to the addition of extra soft or collinear gluons.
i) Breakup of a straight string We here first describe the breakup of a straight flux tube or string between a quark and an antiquark. For simplicity we limit ourself to the situation with a single hadron species, neglecting also transverse momenta. For a straight string stretched between a quark and an anti-quark, the breakup to a state with n hadrons is in the model given by the expression: Here p i and P tot are two-dimensional vectors. The expression is a product of a phase space factor, where the parameter N expresses the ratio between the phase space for n and n − 1 particles, and the exponent of the imaginary part of the string action, bA.
Here b is a parameter and A the space-time area covered by the string before breakup (in units of the string tension κ). This decay law can be implemented as an iterative process, where each successive hadron takes a fraction z of the remaining light-cone momentum (p ± = E ± p z ) along the positive or negative light-cone respectively. The values of these momentum fractions are then given by the distribution Here a is related to the parameters N and b in eq. (3.1) by normalization. (In practice a and b are determined from experiments, and N is then determined by the normalization constraint.) In applications it is also necessary to account for different quark and hadron species, and for quark transverse momenta. The result using Schwinger's formalism for electron production in a homogenous electric field gives an extra factor exp(−π(µ 2 + p 2 ⊥ )/κ), where µ and p ⊥ are the mass and transverse momentum for the quark and anti-quark in the produced pair. For details see e.g. ref. [66].
The result in eq. (3.2) is in principle valid for strings stretched between partons produced in a single space-time point, and moving apart as illustrated in the space-time diagram in fig. 2. The expression in eq. (3.1) is boost invariant, and the hadrons are produced around a hyperbola in space-time. A Lorentz boost in the x-direction will expand the figure in the (t + x) direction and compress it in the (t − x) direction (or vice versa). Thus the breakups will be lying along the same hyperbola, and low momentum particles in a specific frame will always be the first to be produced in that special frame.
The typical proper time for the breakup points is given by This does, however, not necessarily correspond to the "hadronization time", which might also be defined as the time when a quark and an anti-quark meet for the first time to form a hadron.
With parameters a and b determined by tuning to data from e + e − annihilation at LEP, and κ equal to 0.9-1 GeV/fm, eq. (3.3) gives a typical breakup time 1.5 fm, while the average time for the hadron formation is 2 fm. The typical hadronization time can therefore be estimated to 1.5-2 fm. This is important to keep in mind, as this value sets an upper limit on the time available for strings to interact.
ii) Hadronization of gluons An essential component in Lund hadronization is the treatment of gluons. Here it is assumed that the width of the flux tube can be neglected, and that its dynamics can be approximated by an infinitely thin "massless relativistic string". 4 A quark at the endpoint of the string, carrying energy and momentum, moves along a straight line, affected by a constant force given by the string tension κ, reducing (or increasing) its momentum. A gluon is treated as a "kink" on the string, carrying energy and momentum and also moving along a straight line with the speed of light. A gluon carries both colour and anti-colour, and the string can be stretched from a quark, via a set of colour-ordered gluons, to an anti-quark (or alternatively as a closed gluon loop). Thus a gluon is pulled by two string pieces, and retarded by the force 2κ. When it has lost its energy, the momentum-carrying kink is split in two corners, which move with the speed of light but carry no momentum.
A simple example is shown in figure 3. It shows a quark, a gluon, and an anti-quark moving out from a single point, called O in the figure. They move outwards at right angles with energies 2, 2, and 3 units respectively. After 1 unit of time (equal to κ −1 energy units) the gluon has arrived at point A and lost all its energy. The gluon is then replaced by two corners connected by a straight section. The quark has lost its energy in point B.
In the rest frame of the attached string piece it now turns back, gaining momentum in the opposite direction. In the figure frame it is turning 90 • , and after meeting the string corner at C, it is instead pulled back and loosing energy. The anti-quark turns around in a similar way at point D. If the string does not break up in hadrons, the string evolution will be reversed. The kinks meet at point H, and the whole system collapses to a point J.
(Note that this system is not in its rest frame. ) We note that although the string pieces initially move with a transverse velocity 1/ √ 2, after some time most of the string is at rest (the horizontal string pieces in figure 3. A soft gluon will soon stop and be replaced by a straight section stretched as if it were pulled out between the quark and the anti-quark. This implies that the string hadronization model is infrared safe; a soft gluon will only cause a minor modification on the string motion. The same is also true for a collinear gluon [66]. For a string with several gluons there will also be several new straight string pieces, which become more and more aligned with the directions of the endpoints, as described in ref. [68]. Therefore a string stretched over many rapidity units, and with several soft gluon kinks, will be pulled out in a way much more aligned with the beam axis, before it breaks into hadrons.

The string shoving picture
In the following we will describe the details in our new implementation of the string shoving model. First we recap the main idea as applied to two straight parallel string pieces, then we consider the interaction between two arbitrary string pieces, and describe a special Lorentz frame in which the shoving can be properly formulated. Finally we describe how we discretize the shoving into small fixed-size nudges and how these are ordered in time and applied to the final state hadrons.

Force between two straight and parallel strings
The force between two straight and parallel strings was discussed in ref. [30], We here shortly reproduce the treatment presented there. Just after the production of a string stretched between a quark and an anti-quark, the colour field is necessarily compressed, not only longitudinally but also transversely. They then expand transversely with the speed of light until they reach the equilibrium radius R S ∼ 0.5-1 fm.
As illustrated in figure 1, the electric field E obtained in lattice calculations, and fitted to the Clem formula, eq. (2.4), is also well approximated by a Gaussian The normalization factor N can be determined if the energy in the field (per unit length), given by d 2 ρ E 2 /2, is adjusted to a fraction g of the string tension κ. This gives N 2 = 2gκ/(πR 2 ). As discussed in section 2.3, the simple bag model would give g = 1/2. Due to the uncertainties in determining the properties of the flux tube, we will treat R and g as tunable parameters.
When the colour-electric fields in two nearby parallel strings overlap, the energy per unit length is given by ). Taking the derivative with respect to d ⊥ then gives the force per unit length For a boost invariant system it is convenient to introduce hyperbolic coordinates Near z = 0 we get δz = t δη st , and the force in eq. (4.2) gives dp ⊥ /dt δz = f (d ⊥ ). Boost invariance then gives the two equations We have here assumed that the flux tubes are oriented in the same direction, leading to a repulsion. We have also argued in terms of an Abelian field, which means that if the fields are oppositely oriented there would be an attenuation of the fields rather than a repulsion. Since QCD is non-Abelian, the picture is slightly more complex, but the calculations are still valid. In case of two triplet fields in opposite directions, we get with probability 8/9 an octet field, which also leads to a repulsion when compressed. Only with probability 1/9 we get a singlet field, and in this case the strings are assumed to be removed through "colour reconnection", as described in ref. [31]. Also for strings in other colour multiplets the string-string interaction is dominantly repulsive. This is not in Figure 4: The time-evolution of a string piece between two partons in an arbitrary Lorentz frame. The system is rotated so that the momenta (the magnitude of which are indicated by the lengths of the arrows) of the of the partons are in the x − z plane with equal but opposite angle w.r.t. the z-axis. The horizontal lines indicate the extension of the string horizontally and transversely at different times steps. The dotted curve indicates the hyperbola where in x the proper time of the sting piece is equal to R S for a given z = t cos θ 2 .
conflict with the Abelian approximation, as discussed in section 2.3. A qq string is a colour singlet, where the quark is a coherent mixture of red, blue, and green (with corresponding anti-colours for the anti-quark). Similarly the endpoint of an octet string has a coherent combination of the 8 different colour charges.

String motion and the parallel frame
For a string piece between two (massless) partons, the motion and expansion of the sting is very simple in the rest frame of the two partons. If the partons have momenta along the x-axis, the position of the string ends are simply x ± (t) = (t; ±t, 0, 0), where we note that the ends move by the speed of light irrespective of the momentum of the partons. If we instead go to an arbitrary Lorentz frame we can also obtain a simple picture by rotating the partons so that they lie in the x − z plane with the same but opposite angle θ/2 with the z-axis, as shown in figure 4. Here the string ends still move by the speed of light and the position of string ends are given by x ± (t) = (t; ±t sin θ 2 , 0, t cos θ 2 ). A straight relativistic string is boost invariant and has no longitudinal momentum (similar to a homogeneous electric field). The energy and transverse momentum are given by figure 4 is therefore still a straight line; it is just boosted transversely with velocity v = cos(θ/2). If we now want to study the interaction between two arbitrary string pieces, it is not possible to find a Lorentz frame where both these are always parallel to the x-axis, but it turns out that it is possible to find a frame where both always lie in parallel planes, perpendicular to the z-axis. To see this, we first we rotate the string piece in figure 4 an angle φ/2 around the z-axis. Then we want another string piece with angle π − θ/2 w.r.t. the the z-axis but rotated an angle −φ/2. We then have the situation shown in figure 5 where the four partons have momenta (using pseudorapidity instead of the polar angle, to simplify the notation) Here we have six unknown quantities and, for four massless momenta we can construct six independent invariant masses, s ij . This means that for any set of four massless particles we can (as long as no two momenta are completely parallel) solve for the transverse momenta We note that there is a mirror ambiguity in the solution, but apart from that we can now construct a Lorentz transformation to take any pair of string pieces to the desired frame, which we will call their parallel frame.
The sketches in figure 5 show the case where the four partons are produced in the same space-time point, which in general is not the case. The partons from the shower and MPI-machinery in PYTHIA8 are all assigned a space-time positions (0; bx, by, 0) in the lab frame assuming the standard picture that at t = 0 they are packed together at z = 0 and only with transverse separation. When a pair of string pieces are Lorentz transformed into their parallel frame, we assume that their respective production points are a simple average of the positions of the parton ends, giving us a p a 0 = (t a 0 ; x a 0 , y a 0 , z a 0 ) for the string piece moving along the z-axis and the corresponding p b for the piece going in the other direction. From this we get for any given time, t, in the parallel frame that the string piece travelling along the z-axis has the length, 2(t − t a 0 ) sin θ 2 , and lies in a plane transverse to the z-axis. The string piece travelling in the opposite z direction will similarly have a length 2(t − t b 0 ) sin θ 2 , which may be different, but the string will still always lie in a plane perpendicular to the z-axis. The endpoints of the two strings as a function of time then become We see now that the distance between the two planes will change linearly with time as Besides the string motion we are also interested in the radius of the string. As indicated in figure 4 this radius is not constant along the string, but depends on the proper time of a point on the string. As the partons are assumed massless, the endpoints of the string always have zero proper time and the colour field there has not had time to spread out transversely. Looking at a point on the x-axis,x, we can easily find the proper time (4.10) The radius of the string will then vary linearly with τ , from zero in the ends until it reaches the final equilibrium radius, R S = τ S 1 fm. After this the string's width is fixed (as is indicated in figure 4) until it ultimately brakes, which on the average happens at τ H 2 fm. In our implementation we have chosen to only allow the shoving to take place between string pieces at points in the parallel frame where both strings have proper times between τ S and τ H .
Clearly there may also be shoving between string pieces where they have not yet reached their maximum radius, R S . From the derivation of eq. (4.2) we see that the force will be larger and the range will be smaller. In our current implementation it is possible to set τ S < R S to allow for shoving also in these region, but the force is still given for R = R S in eq. (4.2). This can therefore only give an indication of the effect, and we have to postpone a quantified study to a later publication.
The shoving naturally stops when the string breaks, but there is a grey-zone after the string breaks and before the hadrons are fully formed where one could imagine that the string pieces in the hadrons being formed will still repel each other. Also, after the hadrons are fully formed we expect final state interactions between them. In this article we will not investigate final-state hadron interactions, although it is able to produce a sizable flow signal in PbPb collisions with Angantyr initial conditions [69], and has recently been implemented in PYTHIA8 [70]. As indicated by the discussion about the inherent ambiguity in defining a "hadronization time" in section 3, defining a transition between a string-dominated final state at early times, and a hadron dominated final state at late times, will require scrutiny.
Another mechanism that reduces the shoving is when the endpoints have limited momenta, k x in the x-direction in the parallel frame. A parton loses momentum to the string at a rate tκ, and a gluon being connected to two string pieces loses momentum at twice that rate. After a time t = |k x |/2κ a gluon will therefore have lost all its momentum to the string and turn back, gaining momentum from the string in the opposite direction. As explained in section 3 ( figure 3), this means that a new string region opens up which is then not in the same parallel x − y plane. In this article we do not treat this new string region, but will comment on them in the following sections.

Generating the shoving
We now want to take the force between two string pieces in eq. (4.2) and apply it in the parallel frame. In ref. [30] everything was done in the laboratory frame and all strings were assumed to be parallel to the beam, and there generation of the shove was done by discretizing time and rapidity, calculating a tiny transverse momentum exchange in each such point. Here we instead use the parallel frame and discretize the transverse momentum into tiny fixed-size nudges.
Going back to the case where we have two completely parallel strings we have from the expression for the force from eq. (4.2) as and we get the total transverse momentum push on the strings as where we note that the integration limits in x are time-dependent. Now we will instead nudge several times with a fixed (small) transverse momentum, δp ⊥ according to some (unnormalized) probability distribution P (t), which would give us the total push and for small enough δp ⊥ we can make the trivial identification (4.14) In any string scenario we can now order the nudges in time (in the parallel frame), and we can ask the question which of the pairs of string pieces will generate a transverse nudge of size δp ⊥ first. This gives us a situation that looks much like the one in parton showers where the question is which parton radiates first. And just as in a parton shower we will use the so-called Sudakov-veto algorithm [2,71]. The main observation here is that the probability of nothing to happen before some time t exponentiates, so that the (now normalized) probability for the first thing to happen at time t is given by (4.15) In the Sudakov-veto algorithm, one can then generate the next thing to happen in all pairs of strings individually, and then pick the pair where something happens first. In addition, since P (t) may be a complicated function, one may chose to generate according to an simplified overestimate,P (t) ≥ P (t), for which generating according to eq. (4.15) is easier. In this way we get a trial time, t t , for the first thing to happen and then for that particular time calculate the true value of P (t t ), and accept the chosen time with a probability P (t t )/P (t t ). If we reject, we then know we have under estimated the probability of nothing haven happened before t t , which means we can generate the next trial time from t min = t t (in eq. (4.15)). In our case we will make the overestimate by treating the strings as being completely parallel, only separated in z and overestimating the limits in the x-integration. We then generate a time and a point along the x axis, calculate the actual repulsion force there, to get the true probability for accepting the generated time.
The whole procedure can be seen as discretizing in time with a dynamically sized time step with larger time steps where the force is small, and vise versa. This is a very efficient way of evolving in time, and efficiency is important since we will calculate several nudges in each pair of string pieces in an event, and in the case of AA, there may be up to O(10 4 ) string pieces per event at the LHC.

Transferring the nudges to the hadons
Each time a nudge has been included somewhere along a string piece, the string is deformed slightly. It will correspond to a gluon kink on the string, however, since the transverse momentum of this kink is small, the δp ⊥ will soon start propagating along the string as indicated in figure 6 (left). A hadron that is produced along the string where there happens to be such a kink will absorb the corresponding nudge in transverse momentum. In practice this is done by finding the point where the proper time of the kink is τ H , the averaged hadronization time, and finding the primary hadron with the corresponding pseudorapidity (which is strongly correlated to the hyperbolic angle of the creation point, c.f. section 2.1) in the parallel frame, as sketched out in figure 6 (right). In this way each nudge is added to two hadrons from each of the two interacting string pieces. It should be noted that the fact that the two hadrons receive an extra transverse momentum in the same direction, will also mean that they come closer together in pseudorapidity. In addition, energy and momentum conservation is achieved by the adjusting the longitudinal momentum of the two hadrons separately in the two string pieces will also (on average) pull the hadrons closer together in rapidity. In the results below in section 6.2 this becomes visible in the overall multiplicity distribution. It should be noted that the deformation of the string must be taken into account, not only in the pair of string pieces where the nudge was generated, but in all pairs involving one of the string pieces, triggering a recalculation of the next nudge in all these pairs. To do this in detail turns out to be forbiddingly time consuming, so instead we estimate an average shift of a string piece only after a certain number of nudges (typically O(10)) and distribute it evenly along the string. 5

Results for simple initial state geometries
The most widely employed models for describing the space-time evolution of a final state interactions of heavy ion collisions as a QGP (after the decay of a CGC as discussed in the introduction), are based on relativistic dissipative fluid dynamics (see e.g. ref. [72] for a review). A key feature of such models, is that the observed momentum-space anisotropy of the final state, originates from the azimuthal spatial anisotropy of the initial state density profile. The final state anisotropy is quantified in flow coefficients (v n 's), which are coefficients of the Fourier expansion of the single particle azimuthal particle yield, with 5 The total number of nudges is proportional to the square of number of string pieces, N 2 S . Requiring recalculation for all affected pairs after each nudge increases the complexity to O(N 3 S ). If we in addition would take into account the detailed geometry change for every previous nudge, would make the complexity O(N 5 S ), which would be forbiddingly inefficient. respect to the event plane (Ψ n ) [73,74]: Here E is the particle energy, p ⊥ the transverse momentum, φ the azimuthal angle and y the rapidity. In this section, we will explore the models' response to initial state geometry in a toy setup without non-flow contributions from jets, i.e. not real events, and thus we do not adapt the experimental calculation procedures involving Q-vectors (see section 6 where we compare to real data), but rather use the definition of the event plane angle from final state particles: Flow coefficients can then be calculated using only final state information as: v n = cos(n(φ − Ψ n )) .

(5.3)
"Toy" systems with known, simple input geometries are therefore better suited for exploring the basic model dynamics. In section 5.1 we set up a toy model for high energy nuclear-nuclear collisions in which the shoving model can be applied, and in section 5.2 we use the toy model to study the high-density behaviour of the shoving model. The parameters of the shoving model are not tuned, but set at reasonable values of g = 0.5, R S = τ S = 1 fm amd τ H = 2 fm.

Isolating flow effects to v 2 in a toy model
In this and the following section, we restrict our study to systems with constrained straight strings. The strings are drawn between uū pairs, with cms energy of 15 GeV, a small Gaussian kick in p x , p y , and p z fixed by energy-momentum conservation. This ensures that all strings will be stretched far enough in rapidity, that a study of final state hadrons with |η| < 1 will not be perturbed by any edge effects. Final state hadrons studied in the figures, are all hadrons emerging from strings breakings (i.e. no decays enabled) 6 . To study the model response in a heavy-ion like geometry, we set up a toy geometry in the shape of an ellipse, drawn between two overlapping nuclei of r Pb = 7.1 fm. The ellipse has a minor axis (β) given by 2β = 2r Pb −b, and a major axis (α) given by 2α = 4r 2 Pb − b 2 , where b is the impact parameter. The elliptic overlap region is filled randomly with strings, given a certain density (ρ). Example events for ρ = 5 fm −2 are shown in figure 7. We note that this configuration is deliberately chosen to maximize v 2 (elliptical flow) at the expense of v 3 (triangular flow), though with a fluctuating initial state geometry, some v 3 will always be present in order to study the response of the shoving model with minimal fluctuations.
The initial anisotropy is given by the usual definition of the participant eccentricity [77], here employed on the strings: n = r 2 cos(nφ) 2 + r 2 sin(nφ) 2 r 2 , where r and φ are the usual polar coordinates of the string centers, but with the origin shifted to the center of the distribution. It is possible to calculate higher order moments of eq. (5.4), but here we will be concerned only with 2 2 {2} = 2 2 , denoted by the shorthand 2 . We begin by studying average quantities as a function of collision centrality, here defined by the impact parameter of the two colliding nuclei, for exemplary values of string density ρ = {4, 8, 12, 20, 30} fm −2 . In figure 8 (left) we show the average number of strings in bins of centrality 7 . The numbers are compared to the number of strings at η = 0 generated by the Angantyr model [33] in Pb-Pb collisions at √ s NN = 5.02 TeV.
The Angantyr model has been shown to give a good description of charged multiplicity 6 We note that even a single string configuration can give rise to shoving effects, if the strings overlaps with itself. Such configurations could possibly arise to the necessary degree in e + e − collisions, though experimental results so far have not shown any indication of flow [75,76]. 7 The centrality is here defined by the impact parameter between two colliding disks needed to give a certain elliptic geometry.  Figure 9: The v 2 flow coefficient with non-flow effects subtracted according to eq. (5.5) for different constant densities (left), and with string densities from Angantyr (right). Data from ALICE [79] added to guide the eye.
at mid-rapidity in AA collisions [78], and this comparison is therefore useful to provide a comparison to realistic string densities at current LHC energies.
In figure 8 (b), the average v 2 as a function of centrality is shown for the same densities as figure 8 (a), as well as for a reference without shoving, with ρ = 2 fm −2 . Again, several observations can be made. First of all, for a fixed centrality v 2 will increase with increasing density. Secondly, the non-flow contribution is rather large -in particular for peripheral centralities. This was also noted in our recent paper on the Angantyr model [33] using full Pb-Pb calculations. Third, for low string densities, the non-flow contribution is so large, that it takes over in peripheral events, meaning that the curve exhibits no turn-over at a characteristic centrality. Due to the toy nature of this setup, this cannot be taken as a firm prediction of the model, but should be investigated further in full simulations.
Since the goal of this section is to study the response to the elliptic geometry, the non-flow contribution will in the following be subtracted. To obtain this, each event is hadronized twice: once with shoving enabled, and once with shoving disabled. For each of the two final states, Ψ n and cos(n(φ−Ψ n )) is calculated, and the final, non-flow subtracted, (inclusive) flow coefficients are given by: v sub n = cos(n(φ − Ψ n )) shoving − cos(n(φ − Ψ n )) no shoving all events .
The flow coefficient v 2 , with non-flow subtracted, is shown in figure 9. In the left part of the figure, v 2 is shown for different constant densities. The higher density samples now all have a characteristic turn-over at some centrality. It should be noted that the density intervals are not evenly spaced. Looking at v 2 in figure 9, the growth in the peak value decreases with higher density, like a saturation effect. Since final state multiplicity is proportional to the number of strings, the result in figure 9 suggests that the shoving model dynamically predicts a steep rise in the peak value of v 2 at low densities (corresponding to low energies, with otherwise fixed geometry), which should then slow down. The fact that string density is different at different centralities (as also shown in figure 8 (a)), is used to construct figure 9 (b). In this figure, string multiplicity for each toy-event are generated with Angantyr, and an elliptic initial condition of the same impact parameter is constructed. This provides more realistic string densities at different centralities, and as it can be seen, this construction does a reasonable job at describing v 2 as a function of centrality. A no-shoving reference is also added to this figure, and as expected, it is consistent with zero. Finally it can be added, both v 3 and v 4 are, as could be expected from the engineered initial conditions, both compatible with zero after the subtraction, indicating that non-flow has been successfully removed. (Not shown in any figure here.)

Towards the hydrodynamical limit at high string density
As mentioned earlier, hydrodynamic calculations are often employed for heavy ion phenomenology, when it concerns interactions of the final state, leading to collective flow. Whether or not the shoving model will, in the end, describe heavy ion data to a satisfactory degree, it is interesting to ask to what kind of hydrodynamics the shoving model will become in the thermodynamic limit. While analytic exploration of this question will be left to future studies, it is possible at this point to study similarities in the phenomenology of the two approaches. An obvious route is to follow the paper by Niemi et al. [80], which studied the correlation between flow coefficients and eccentricities, and noted that v 2 and v 3 have a strong linear correlation with 2 and 3 . In figure 10, we show the correlation between 2 (from eq. (5.4)) and v 2 , for four values of ρ. No centrality selection is performed, as the toy geometry of the system impose a strict relationship between eccentricity and centrality defined by impact parameter, and binning in centrality would thus not add any information. In the very dilute limit, ρ = 2 fm −2 , in figure 10 (upper left), it is seen that even though v 2 in figure 9 is non-zero, there is no strong correlation with 2 , and the distribution of v 2 is rather wide. At intermediate densities ρ = 12 and 20 fm −2 , in figure 10 (upper right) and (lower left), the strong correlation appears, though not as narrow as in ref. [80]. Finally for the densest considered state, ρ = 30 fm −2 in figure 10 (lower right) the correlation is, on average, similar to the more dilute ones, but more narrow.
To further study the scaling of v 2 with 2 , we study the scaled event-by-event variables: In figure 11 the distributions of δ 2 and δv 2 , are shown for the densities ρ = {12, 20, 30} fm −2 , in the 20-30% centrality bin. It should be noted that, due to the purely elliptical sampling region, the shape of the distribution cannot be compared directly to those of ref. [80]. The main conclusion is, however, clear. In the dense limit, in the rightmost panel, the two distributions are almost identical. This means that the shoving model, in this limit, reproduces a key global feature of hydrodynamics, namely full scaling of final state (purely momentum space) quantities with the global initial state geometry. A notable discussion pertains to the issue whether or not heavy ion collisions at RHIC and LHC energies, reach a high enough string density for the shoving model to behave like hydrodynamics, as indicated above. In figure 8 (left), it was shown that Pb-Pb collisions at LHC energies produce an amount of strings corresponding to toy model densities of around -25 -leftmost panels of figure 11) and less than 4 fm −2 for the most peripheral events. Even though the toy model predicts scaling of v 2 with 2 , the fluctuations do not exhibit the same scaling. A further direct study of flow fluctuations in non-central heavy ion collisions would thus be of interest both on the phenomenological and experimental side.

Results with Angantyr initial states
In the previous section, we have shown that given a simple initial state of long, straight strings without any soft gluons, the shoving mechanism can produce a response which (a) scales with initial state geometry in the same way as a hydrodynamic response, and (b) can produce flow coefficients in momentum space at the same level as measured in experiments. In this section we go a step further, and present the response of the model given a more realistic initial string configuration, as produced by the Angantyr model, which can be compared to data. In section section 4 we described several of the challenges faced when interfacing the shoving model to an initial state containing many soft gluons, which in particular is the case in AA collisions. Throughout this section we use the same canonical values of shoving model parameters as in the previous section.

Results in pp collisions
Already the original implementation of the shoving model was shown in ref. [30] to give a satisfactory description of the pp "ridge". We will in this section focus on flow observables as calculated at high energy heavy ion experiments, i.e. flow coefficients calculated using the generic framework formalism [81,82], in the implementation in the Rivet program [83].
It is in principle possible to generate pp events using the normal PYTHIA MPI model [84][85][86]. It would, however, be computationally inefficient, since emphasis should be given to high multiplicity results. The results presented here therefore use the modifications of the MPI framework presented as part of the Angantyr heavy ion model [33], notably the ability to bias the impact parameter selection towards very central pp collisions, and re-weighting back to the normal distribution.
We first note, that while the shoving mechanism does not change the total multiplicity of an event, it will change the multiplicity in fiducial region measured by an experiment, because it will push particles from the unmeasured low-p ⊥ region to (measured) higher p ⊥ . It is therefore necessary to slightly re-tune the model parameters, to obtain a correct description of basic observables such as p ⊥ distributions and total multiplicity. In practise, the parameter p ⊥,0 , which regulates 8 the 2 → 2 parton cross section is increased from the default value of 2.28 to 2.4. In figure 12 we show a comparison between PYTHIA (with the above mentioned impact parameter sampling) and PYTHIA + shoving, for a few standard minimum-bias observables, compared the charged p ⊥ distribution (left), the distribution 8 The divergence of the partonic 2 → 2 cross section is regularized for p ⊥ → 0 by a factor p 4 ⊥ /(p 2 ⊥0 +p 2 ⊥ ) 2 , and by using an αs(p 2 ⊥0 + p 2 ⊥ ). Tuning is done by MultipartonInteractions:pT0Ref which is the pT 0 value for the reference CM energy (where pT0Ref = pT0(ecmRef)). of number of charged particles (middle) and p ⊥ (N ch ) (right), all at √ s = 13 TeV. Data by ATLAS [87], the analysis implemented in the Rivet framework [88]. The agreement between simulation and data for the multiplicity distribution, is not as good as normally expected from PYTHIA. This is expected, as only non-diffractive collisions were simulated. More interesting is the slight difference between the two simulations, where it is clearly visible that in spite of the re-tuning, shoving still produce more particles in the high-N ch limit. In the p ⊥ distribution, it is seen that shoving has the effect of increasing the spectrum around around p ⊥ = 1 GeV. While the normalization if off (due to the exclusion of diffractive events in the simulation), shoving brings the shape of the low-p ⊥ part of the spectrum closer to data. Finally the p ⊥ (N ch ) is almost unchanged.
We now turn our attention to flow coefficients, and show v 2 calculated by two-particle correlations in figure 13. In figure 13 (left) the multiplicity dependence of v 2 {2} with |∆η| > 1.4 is shown, and in figure 13 (right) the p ⊥ -dependence of v 2 {2} (|∆η| > 2) in high multiplicity events is shown. Several conclusions can be drawn from the two figures.
First of all, it is seen that v 2 as a function of multiplicity (in figure 13 (left)) is too high with shoving enabled. We emphasize that the model parameters have not been tuned to reproduce this data, and in particular that successful description of this data will also require a good model for the spatial distribution of strings in a pp collision -a point we will return to in a moment. We do, however, note that the additional v 2 added by the shoving model, persists even with an η separation of correlated particles (as |∆η| cuts are applied), a feature which separates the shoving model from e.g. colour reconnection approaches, which have been pointed out to produce flow-like effects in pp collisions [91][92][93]. We also note that the p ⊥ -dependence of v 2 is drastically improved, as seen in figure 13 (right). In particular the high-p ⊥ behaviour of this quantity is interesting, as it decreases wrt. the baseline when shoving is enabled.
In figure 14 comparisons to v 3 {2} (left) and v 4 {2} (right) are performed. While it is clear that shoving adds a sizeable contribution to both, it is equally clear that data is not very well reproduced. We remind the reader that any v n produced by the shoving model comes about as a response to the initial geometry, and the initial geometry used by default in PYTHIA, consists of two overlapping 2D Gaussian distributions. It was shown in ref. [94],  Figure 13: Comparison to v 2 {2} as function of multiplicity with ALICE high multiplicity trigger (left), and versus p ⊥ in high multiplicity events (right). Data from pp collisions at √ s = 13 TeV by ALICE [89] and CMS [90].  [89].
that applying more realistic initial conditions, can drastically change the eccentricities of the initial state in pp collisions. So while the description at this point is not perfect, the observations that a clear effect is present, bears promise for future studies. Further on, correlations between flow coefficients, the so-called symmetric cumulants [82,95], will be an obvious step. But at this point, without satisfactory description of the v n 's themselves, it is not fruitful to go on to even more advanced observables.
Finally, in figure 15, we show results for the four-particle cumulant c 2 {4}. We briefly remind the reader about some definitions. The 2-and 4-particle correlations in a single The averages are here taken over all combinations of 2 or 4 non-equal particles in one event.
The four particle cumulant is the all-event averaged 4-particle azimuthal correlations, with the 2-particle contribution subtracted: The double average here means first an average over particles in one event, and then average over all events.
As discussed in ref. [96], when the correlation is dominated by flow and the multiplicity is high, then the flow coefficient v 2 {4} is given by 4 −c 2 {4}. Clearly c 2 {4} must be negative for this to be realized. This, in turn, means that the relative difference between the 2 -and 4-particle azimuthal correlations, must be right from eq. (6.2). As it was also pointed out in ref. [96], the non-flow contribution to four-particle correlations is much smaller than for two-particle correlations, as the cumulant becomes flow-dominated when v n 1/M 3/4 (M is the multiplicity) in the former case, but only when v n 1/ √ M in the latter. In a pp collision M is small compared to a heavy ion collision, and it can therefore be reasonably expected that the four particle correlations will only be flow dominated at sufficiently high multiplicity. Since data show a real v 2 {4}, the importance of the sign of c 2 {4} in model calculations for pp, have recently been highlighted [97,98]. Importantly, standard hydrodynamic treatments do not obtain a negative sign of c 2 {4} in pp collisions, even with specifically engineered initial conditions [99].
In the results from the shoving model in figure 15, we note that while a negative c 2 {4} is not produced when comparing to CMS results, it is produced in high multiplicity events in the ALICE acceptance, using the high multiplicity trigger. There are several possible reasons for this apparent discrepancy. The acceptances are quite different, and since the sign of c 2 {4} is an observed characteristic, rather than a fundamental feature of the model, it is difficult to point out why a given model should produce different results in different acceptances -though it is possible. More interesting, is the possible effect of the high multiplicity trigger. In figure 15 (left), it is seen that both default PYTHIA, as well as PYTHIA with the shoving model, over-predicts c 2 {4} at low-multiplicity by a large margin. As noted in the original paper, this is also the case for the 2-particle cumulant. A reasonable explanation for this over-prediction could be, that PYTHIA collects too many particles in mini-jets in general. With a high multiplicity forward trigger, a strong bias against this effect is put in place, and the underlying model behaves more reasonable. In any case, the finding of a negative c 2 {4} in high-multiplicity events with the shoving model is an interesting and non-trivial result, which will be followed up in a future study.

Results in Pb-Pb collisions
We now turn to Pb-Pb collisions, where we use the Angantyr model in PYTHIA keeping the same settings as in section 6.1. The results for Pb-Pb collision events at 5.02 TeV has been compared to ALICE data points via Rivet routines. 9 The anisotropic flow coefficients plotted here have been calculated, as in the previous section, using multi-particle cumulant methods as done in the ALICE experiment.
Centrality measures used in these analyses are of two kinds: for the plot in figure 16 we use the centrality binning of the generated impact parameter by Angantyr, and for the plots in figure 17, we use Angantyr generated centrality binning which mimics the experimental centrality definition where in ALICE the binning is in the integrated signal in their forward and backward scintillators. However, the difference between the two centrality measures is small in Pb-Pb collisions [83].
In figure 16, we plot the charged particle multiplicity for seven centrality classes (0-5%, 5-10%, 10-20%, 20-30%, 30-40%, 40-50%, 50-60%) as a function of pseudorapidity in the range −3.5 < η < 5 for √ s N N = 5.02 TeV in Pb-Pb collisions comparing it to the study performed by ALICE [100]. We use this figure as our control plot to check that when we turn shoving on, the description of other observables are not destroyed. We observe that this implementation fairly well preserves the Angantyr description of the multiplicity distributions. The overall multiplicity of the shoving curve is however a bit lower when compared to default Angantyr, which is because of the increased pT0Ref as mentioned in 6.1. Also, as discussed in section 4.3, when strings are shoved and the particles on average get a larger p ⊥ which also means that they come closer together in pseudorapidity. The overall effect is that particles are generally dragged closer towards mid-rapidity, reducing the two-humped structure seen for plain Angantyr.
We will look into further improvement of the multiplicity description by shoving through tuning, normalization of the distribution functions and accurate description of centrality as in experiments in the future. Figure 17 presents the centrality dependence of the harmonic flow coefficient v 2 from two-particle cumulant on the left with |∆η| > 1 and four-particle cumulants on the right,  integrated over the p ⊥ range 0.2 < p ⊥ < 5.0 GeV for 5.02 TeV Pb-Pb collisions [79] for generated centrality. We note that Angantyr with shoving results in an increased v 2 in the right direction with respect to data. We see in data that v 2 {2} increases from central to peripheral collisions, reaching a maximum of around 0.10 between 40 − 50% centrality. v 2 {4} also shows a similar behaviour with a maximum around 0.09 between 40 − 50% centrality. 10 String shoving result clearly lacks the curvature of the data points, but doubles the contribution in v 2 {2} as compared to Angantyr. The underlying cause for this behaviour is that the current implementation of shoving alone is not sufficient to generate the large overall response to the anisotropy in the initial collision geometry of the nuclei. An increased g factor or delayed hadronization time or an early onset of shoving, and a combination of these factors, do not give rise to enough v 2 either. In section 5.1, we showed that shoving can generate sufficient v 2 as seen in data with completely straight strings without any gluon kinks. With Angantyr, we have more realistic final states with many and often soft gluon emissions from the multi-parton interactions and the initial-and final-state evolution models, which hinder the process of strings shoving each other by cutting short their interaction time, hence resulting in the overall observation of the lack of enough transverse nudges generated via this mechanism.

Conclusions and outlook
In this paper we have argued that a hot thermalized plasma is not necessarily formed even in central AA collisions, not even in Pb-Pb collisions at the highest attainable energies at the LHC. Instead we note that the string-based approach to simulating hadronic final states in the Angantyr model in PYTHIA8 gives a very reasonable description of the number and general distribution of particles in AA events, and take this as an incentive to study string hadronization in dense collision systems more carefully.
Our string picture is qualitatively different from the more conventional picture, where the colliding nuclei are described in terms of a CGC, that in the moment of collision turns into a so-called glasma, which very soon decays into a thermalized QGP. Similar to the string picture, the glasma has longitudinal fields stretched between the nucleus remnants, and these fields are kept together in flux tubes as the remnants move apart. There are, however, also essential differences. The glasma turns rapidly into a thermalized plasma. Such a plasma expands longitudinally in a boost invariant way, with decreasing energy density as a result. The initial density must therefore be quite high to give the observed particle density after freezeout. In contrast the energy density in the strings is constant up to the time for hadronization. When the strings become longer, the energy in the new string pieces is taken from the removing nucleus remnants, and not from depleting the energy in the strings already formed. In the string scenario we estimate the energy density at mid-rapidity to around 5 GeV/fm 3 (in PbPb at the LHC), while in the glasma we find that it ought to be one or two orders of magnitude higher.
The low energy density in the string scenario implies that the vacuum condensate is very important to form the strings. The break-down of the glasma is often motivated by the so-called "Nielsen-Olesen instabilities". These authors showed that a longitudinal chromo-electric field added to the QCD Fock vacuum is instable, and transverse fields grow exponentially. This growth does not go on forever. Instead higher order corrections will lead to a Higgs potential analogous to the potential describing the condensate of Cooperpairs in a superconductor or the vacuum condensate in the EW Higgs model. Adding a (not too strong) linear field to this non-trivial ground state will then give flux tubes, similar to the vortex lines in a superconductor. We take this as an indication that it may indeed be possible that strings can be formed and actually survive the initial phase of the collision, without a thermalized plasma being formed. (The energy density needed in the glasma may, however, be strong enough to destroy the "superconducting" phase.) Another important difference between the two scenarios is that the glasma contains both chromo-electric and chromo-magnetic fields, not only chromo-electric fields as in the string picture. This implies CP-violating effects in the glasma, but this feature is not discussed in this paper.
In an earlier paper we argued that in a dense environment where the strings overlap in space-time, they should repel each other and we showed with a very simple model that this could induce flow in pp. Here we have motivated the model further, and compared to lattice calculations to estimate the transverse shape and energy distribution of the stringlike field. We have also improved the implementation of the model, where the strings no longer have to be completely parallel in order to calculate the force between them. Instead we show that for a pair of arbitrary string pieces, we can always find a Lorentz frame where they will be stretched out in parallel planes, allowing us to easily calculate the force there. We have further improved the time discretization, and instead of processing string interactions in fixed time intervals, the shoving model is implemented as a parton shower with dynamical time steps, greatly improving the computational performance of the model. Finally we have improved the procedure for transferring the nudges from string interactions to final state hadrons.
The implementation used for obtaining the results presented here is not yet quite complete. Although it circumvents the production of huge amounts of soft gluons in the shoving, which was a major problem for our previous implementation, it still has a problem with dealing with the soft gluons that are already there from the initial-and final-state parton showers. Gluon kinks loose energy with twice the force compared to a quark. A soft gluon therefore will soon loose all its energy, and new straight sections of the string will be formed, which in the current implementation are not taken into account.
In addition the implementation only allows shoving at points in space-time after the strings have expanded to the equilibrium size, R S ∼ τ S , and before they start to break up in hadrons at proper times around τ H . Clearly shoving should be present also at times before τ S , and the force between very close strings should then be higher, but our implementation currently cannot handle situations with varying string radii. Also, at times later than τ H , after the string break-up, one could expect some shoving between the hadrons being formed, and after they are formed one needs to consider final state re-scattering.
In light of these shortcomings of the current implementation, it is not surprising that when we apply the model to complete partonic final states generated by Angantyr, we cannot quantitatively reproduce the amount of v 2 measured in experiments. But we do see that the shoving actually does give rise to flow effects. We also see that in pp we currently get a bit too large v 2 , but we would like to emphasise that we have here not tried to do any tuning of the model parameters. Instead these are kept at canonical values, To investigate the model further and to make it plausible that the shoving model, when implemented in full, actually may be able to reproduce also quantitatively the flow effects measured in AA collisions, we looked at what happens when applying it to a toy model of the initial state. To account for the missing string pieces due to soft gluons, we here used parallel strings without gluon kinks. The strings are randomly distributed with variable density in an ellipsoid shape in impact parameter space. In this way the shoving was unhampered by soft gluons, and we found that the model is able to get the azimuthal anisotropies in momentum space as expected from the eccentricity of the shape in impact parameter. By then matching the string density to the one we have in Angantyr for different centralities, we could also see that the resulting v 2 was much closer to measured data.
Showing that we can get reasonable azimuthal anisotropies in AA collisions using a purely string-based scenario is, of course, not enough to prove that a thermalized QGP is not formed in such collisions. To do this we need to also be able to describe other measurements, such as strangeness enhancement and jet quenching and, more importantly find new observables where a string based scenario predicts results that cannot be reconciled with the QGP picture. And for this we not only need to improve the implementation of the shoving model, but also revisit our rope model and also our "swing" model for colour reconnections. The rope model has been shown to give reasonable descriptions of strangeness enhancement in high multiplicity pp, and using the parallel frame presented here, we should be able to get a better handle on the space-time picture of the string overlaps needed to be able to apply it in AA. Also for the swing model we can take advantage of the parallel frame to properly understand which partons may reconnect and when and where they may do so. And since in the parallel frame hard and soft partons are treated on an equal footing this could also have interesting effect on jets.
In the end we hope that these models will be implemented in PYTHIA8 together with the Angantyr model, so that we get a complete platform for generating fully hadronic final states that can be compared to any type of measurement in any kind of collision (AA, pA, pp, . . . ). This would then give us a perfect laboratory to investigate a purely string-based picture as an alternative to the conventional QGP approach.

A. Vortex lines in a superconductor
The microscopic properties in a superconductor, the magnetic field, H, and the condensate wavefunction (order parameter), ψ, can be determined by the LG equations (see e.g. ref. [52]). The theory contains two characteristic lengths, the penetration depth λ for the magnetic field and the coherence length ξ for the condensate. In its generalization to a relativistic theory, the Lagrangian contains a Higgs potential for the condensate. Here the lengths λ and ξ correspond to the (inverse) masses of the gauge boson and the Higgs particle respectively.
For a flux tube the wavefunction ψ is undetermined along a "vortex line", and has a phase changing by 2πn when going around the vortex line, with n an integer. The total flux in the flux tube is then quantized to n times a flux quantum Φ 0 = 2π/q, where for a normal superconductor q = 2e is the charge of a Cooper pair. This quantum would also correspond to the charge of a magnetic monopole. The change in phase is related to a vortex-like current in the condensate, which keeps the flux confined within the flux tube.
At the boundary between a normal and a superconducting phase, the pressure from the condensate and the magnetic field balance each other. The condensate goes to zero over a distance ξ in the superconductor, and the magnetic field is suppressed over a distance λ. As a result, when ξ is larger than λ (or more exactly ξ > √ 2 λ), both the condensate and the field are suppressed over a range ξ − λ. This is a type I superconductor, and it implies that the surface provides a positive contribution to the energy. In equilibrium the surface then tends to be as small as possible. If in contrast λ is larger than ξ (type II superconductor), the condensate and the field can coexist over a range ∼ λ − ξ, and the surface provides a negative contribution to the energy, favouring a large surface. If the flux tube has more than one flux quantum, there is then a tendency to split it into a number of vortices, each with one unit of flux. In case of a large total flux, there is a repulsive force between nearby flux tubes. The system will then tend to expand, forming an "Abrikosov lattice".
The interaction between the condensate and the electromagnetic field in a superconductor is described by the LG equations, which in its relativistic generalization corresponds to the Abelian Higgs model relevant for the Abelian projection of the QCD field. The Lagrange density is here given by When the parameter α is smaller than zero, the scalar field ψ forms a condensate ψ = ψ 0 = −α/β. The mass of the Higgs particle and the massive gauge boson, given by √ −2α and e −2α/β respectively, correspond to (the inverse of) the coherence length ξ and the penetration depth λ. The LG equations are obtained from Euler-Lagrange's equations varying ψ (or ψ * ) and A µ .
In an extreme type II superconductor with ξ λ, the LG equations have a solution ψ = const. e −iφ (for ρ > ξ), for a vortex line which carries one unit of flux. The corresponding magnetic field is here given by H(ρ) = C · K 0 (ρ/λ), (A.2) where C = Φ/(2πλ 2 ) is a constant, Φ is the total flux, and K 0 is a modified Bessel function. The Bessel function has a logarithmic dependence on ρ for ρ < λ, but falls exponentially for ρ > λ. The field is confined within this range by an electric current j = λ |curl H| = (C/λ) K 1 (ρ/λ) (also valid for ρ > ξ). In this extreme case, when ξ is very small, the contribution to the energy from destroying the condensate is also small, and the energy of the flux tube, the string tension κ, is given by the sum of the field energy and the energy in the current: We note that the energy is dominated by the contribution from the current, where K 1 (ρ/λ) is singular and ∼ λ/ρ for small ρ. Thus the total energy is proportional to ln(λ/ξ) for very small ξ. When ξ is smaller than λ, but not close to zero, the following approximate solution is given by Clem [62]: The new distance scale ξ v depends on the ratio κ LG ≡ λ/ξ, and is close to ξ for κ LG ≈ 1/ √ 2, i.e. for a superconductor on the border between type I and type II. The expression in eq. (A.4) satisfies the one of the LG equations (obtained by varying the field A µ ) and Ampèr's equation j = curl B = curl curl A. It does, however, not satisfy the equation obtained by varying ψ, and for ξ > λ this equation can be badly violated.
In a superconductor there are magnetic flux tubes, and magnetic monopoles would be confined. In QCD colour-electric flux tubes and colour-electric charges are confined. Thus for the Abelian projection the fields F µν will be replaced by the dual fields F µν , as discussed in section 2.2.