Baryon stopping as a relativistic Markov process in phase space
BBaryon stopping as a relativistic Markov process in phase space
Johannes Hoelck and Georg Wolschin ∗ Institute for Theoretical Physics, Heidelberg University,Philosophenweg 12–16, 69120 Heidelberg, Germany, EU (Dated: September 21, 2020)We reconsider baryon stopping in relativistic heavy-ion collisions in a nonequilibrium-statisticalframework. The approach combines earlier formulations based on quantum chromodynamics with arelativistic diffusion model through a suitably derived fluctuation–dissipation relation, thus allowingfor a fully time-dependent theory that is consistent with QCD. We use an existing framework forrelativistic stochastic processes in spacetime that are Markovian in phase space, and adapt it toderive a Fokker–Planck equation in rapidity space, which is solved numerically. The time evolutionof the net-proton distribution function in rapidity space agrees with stopping data from the CERNSuper Proton Synchrotron and the BNL Relativistic Heavy Ion Collider.
I. INTRODUCTION
In relativistic heavy-ion collisions at the CERN SuperProton Synchrotron (SPS), the BNL Relativistic HeavyIon Collider (RHIC), or the CERN Large Hadron Col-lider (LHC), the incoming baryons are being slowed down(“stopped”) as they interpenetrate each other, while inthe spatial region between the receding, highly Lorentz-contracted fragments [1] a hot fireball is formed, whichcools during its expansion and eventually hadronizes ina parton–hadron crossover. Of particular interest is theinitial stage of such a collision with the local thermaliza-tion of quarks and gluons, and the simultaneous stoppingof the baryons. The latter occurs essentially through col-lisions of the incoming valence quarks with soft gluons inthe respective other nucleus.Various models to account for the stopping processand its energy dependence have been developed, for ex-ample, in Refs. [2–4] and related works, which are rely-ing on the appropriate parton distribution functions andhence on quantum chromodynamics (QCD). These mod-els yield agreement with the available stopping data atSPS and RHIC, such as the distributions of net protons(protons minus produced antiprotons) in longitudinal ra-pidity space, and also provide predictions at LHC ener-gies, where stopping data at forward rapidities are notyet available. They do not, however, provide the timedevelopment from the initial distribution at the instantof the collision to the final, measured one.Complementary time-dependent approaches to stop-ping and local equilibration have relied on phenomeno-logical, nonequilibrium-statistical approaches: A linearFokker–Planck equation for the net-baryon rapidity dis-tribution function had been proposed in Ref. [5], whichaccounts for the time evolution of the net-baryon or net-proton rapidity distributions in a two-source relativisticdiffusion model (RDM). Variants of the model with anonlinearity in the diffusion term have subsequently beensuggested in Refs. [6–8] and related works, which assume, ∗ [email protected] however, the debatable validity [9] of nonextensive statis-tics. A linear diffusion model for particle production wasalso put forward in Ref. [10] and subsequent works. Forproduced particles at RHIC and LHC energies, a third(midrapidity) source is essential to cover pair productionprocesses in the central fireball that provide the bulk ofcharged-hadron generation at sufficiently high energy [9].When considering net baryons or protons, however, themidrapidity source cancels out because it is equally com-posed of particles and antiparticles. The highly nonlinearlocal thermalization of quarks and gluons in the initialstages of the collision can be modeled through quantumBoltzmann-like collision terms, which require numericalsolutions [11, 12], but also a schematic model has beendeveloped that accounts for the fast local equilibrationthrough analytical solutions of a nonlinear partial differ-ential equation [13].Diffusion models are being used in many areas ofphysics, chemistry, and biology [14, 15]. They have orig-inally been developed by Einstein and Smoluchowski toprovide a mesoscopic theory of Brownian motion [16–18]as linear differential equations for the Brownian particles’single-particle distribution function. Alternative micro-scopic approaches treat the Brownian particles’ trajecto-ries as stochastic processes in position space, for exam-ple in the form of a Wiener [19] or Ornstein–Uhlenbeckprocess [20]. Stochastic processes and stochastic differ-ential equations have subsequently been considered inmore generality in a newly established branch of math-ematics, stochastic calculus, with notable contributionsby Itˆo [21, 22], Stratonovich [23], Fisk [24], and Klimon-tovich [25] who introduced various concepts of a stochas-tic integral, each with different mathematical propertiesand physical interpretations. In some cases, connectionsbetween the micro- and mesoscopic formulation can beestablished through a Kramers–Moyal expansion [26, 27]or the Feynman–Kac formula [28].Following the discovery of special relativity, it becameclear that statistical physics in general, and diffusionmodels in particular, would have to be adapted [29] tomeet the requirements imposed by a limited velocity oflight. Especially, nontrivial Lorentz-invariant stochasticprocesses for spacetime coordinates are necessarily non- a r X i v : . [ nu c l - t h ] S e p Markovian [30–32], which makes a straightforward gener-alization of the nonrelativistic diffusion equation to spe-cial relativity impossible. While a general stochastic pro-cess may depend on any finite or infinite number of itsprior realizations, Markov processes [33] have no memoryof the past apart from their current state, which greatlysimplifies their mathematical treatment. Consequently,stochastic processes used in physical models often havethe Markov property, in spite of being mathematicallythe exception rather than the rule.As long as a process’ memory is finite, i. e., only a fi-nite number of its previous realizations affect the nextvalue, it can be reformulated as a coupled system of mul-tiple Markov processes through the introduction of addi-tional variables [34]. This has been used in Refs. [35–40] to formulate relativistic phase-space diffusion pro-cesses based on a generalized Ornstein–Uhlenbeck pro-cess for the Brownian particle’s momentum. These pro-cesses are Markovian in phase space but lose the Markovproperty when expressed solely in spacetime coordinates.The authors also deduced associated relativistic Kramersand Fokker–Planck equations for the particles’ phase-space and momentum distribution functions, and derivedfluctuation–dissipation relations suitable for an isotropicthermal background.In the present work, we aim to derive a nonequilibrium-statistical diffusion model for baryon stopping in rapid-ity space that is based on the key premises of the phe-nomenological RDM, but is constructed from a consis-tent approach with relativistic Markov processes in phasespace and incorporates the QCD-based theory through asuitably adapted fluctuation–dissipation relation. Thecorresponding Fokker–Planck equation will enable us toaccount for the time evolution of the initial distributionfunctions from the onset of a relativistic heavy-ion colli-sion to the final, measured distributions of net baryonsor net protons in agreement with the available SPS andRHIC data.The key assumptions for our nonequilibrium-statisticalapproach to stopping in relativistic heavy-ion collisionsare presented in the next section, followed by the equa-tions of motion in Langevin and Fokker–Planck formu-lation in Sec. III. Drift and diffusion terms are discussedin Sec. IV, as well as the expected stationary state de-rived from the earlier QCD formulation that allows us toformulate an appropriate fluctuation–dissipation relationthat determines the course of the time evolution from theinitial to the final net-proton rapidity distribution func-tions. The latter are compared with available stoppingdata from SPS and RHIC experiments in Sec. V. Theconclusions are drawn in Sec. VI.
II. NET-PROTON RAPIDITY SPECTRA
We model baryon stopping as a diffusive process inrapidity space of the participating nucleons whose dy-namics are governed by a—not necessarily thermalized— fluctuating background representing the quarks and glu-ons of the fragments. In this process, interactions be-tween the partons and the background are assumed toprevail such that the nucleon distribution function re-duces to a superposition of single-particle probabilitydensity functions. As the distribution functions of par-ticipant baryons are experimentally inaccessible, we con-sider only participant protons in our model and compareour result to the measured net-proton number densityin rapidity space, i. e., the difference between the dis-tributions of protons and antiprotons produced in thecollision, which we expect to be reasonably close to theparticipant-proton distribution function.To incorporate the spatial separation of the two nu-clear fragments, we use a two-source ansatz [5] and com-pletely disconnect the time evolution of particles origi-nating from the forward- and backward-moving fragmentthrough separate probability densities and fluctuation–dissipation relations. Taking advantage of the symmetryof the system with respect to its center of momentum,we then write the net-proton number density in rapid-ity space d N p − ¯ p / d y in the system’s center-of-momentumframe F as the superpositiond N p − ¯ p d y ( t ; y ) ≈ N p − ¯ p ψ ( t ; + y ) + ψ ( t ; − y )] , (1)where N p − ¯ p denotes the net-proton number and ψ ( t ; ± y ) d y the probability to find a participant protonfrom the forward- or backward-moving fragment, respec-tively, at time t with rapidity in [ y, y + d y ]. A. Initial state
Prior to the collision of the nuclei at some time t i , weassume the system to be in an initial state where each nu-cleus can be approximated by a zero-temperature Fermigas with appropriate Fermi momentum p F . Then, theprotons of each nucleus are distributed in the nucleus’srest frame F ∗ according to the momentum-space proba-bility density function φ i ( (cid:126)p ∗ ) = 34 πp Θ( p F − | (cid:126)p ∗ | ) , (2)which is given by a Heaviside step function Θ scaled by anormalizing factor. We determine the Fermi momentumthrough a simple potential well model, [41] p F = (cid:114) π ZV ∗ , V ∗ = 4 π r ∗ . (3)Here, Z denotes the nucleus’s proton number, V ∗ thenuclear charge volume, and r ∗ the nuclear charge radius,which we take from Ref. [42]. ψ i ( y * ) y * FIG. 1. Marginal rapidity probability density function ψ i (solid, red) of the participant protons in a Pb nucleus priorto the initial collision. For comparison, a normal distributionwith zero mean and standard deviation y F / ± y F . Choosing the orientation of F ∗ such that p ∗ is parallelto the beam axis, we define the cylindrical coordinates γ ⊥∗ = (cid:112) p ∗ /m ) + ( p ∗ /m ) , (4a) ϕ ∗ = arctan( p ∗ /p ∗ ) , (4b) y ∗ = atanh( p ∗ /p ∗ ) , (4c)with the proton mass m . Boosting to F leaves the trans-verse degrees of freedom unaffected ( γ ⊥ = γ ⊥∗ , ϕ = ϕ ∗ )while the longitudinal rapidity coordinate is shifted bythe beam rapidity y b , y = y ∗ ± y b . Integrating out γ ⊥ and ϕ , the initial rapidity-space probability density in F , ψ i ( y ) ≡ ψ ( t i ; y ), is found to be ψ i ( y ∗ ± y b ) = 12 sinh( y F ) − Θ( y F − | y ∗ | ) × cosh( y ∗ ) (cid:34)(cid:18) cosh( y F )cosh( y ∗ ) (cid:19) − (cid:35) (5)with the Fermi rapidity y F = asinh( p F /m ). Fig. 1 shows ψ i as a function of y ∗ for the isotope Pb. The nu-merical values of y F for Au and
Pb are 0.3134 and0.3136, respectively, and differ only slightly in the fourthdecimal place.A more realistic description of the initial state is prin-cipally desirable (for example, with a finite temperature);however, the exact form of the initial probability densityfunction hardly influences the later stages of the timeevolution. Hence, normal or delta distributions are oftenused as convenient approximations for the initial state[43].
B. Final state
For t > t i , the system evolves in time, driven by thefluctuating background, until it reaches a final state at some time t f when the partonic interactions between thenuclei effectively cease due to their increasing spatial dis-tance. Consequently, for comparison with experimentaldata from SPS and RHIC, we evaluate Eq. (1) at t = t f .As we will see in Sec. V, the concrete value of t f is notimportant at this stage of the model, since it does onlyappear in products with other a priori unknown quanti-ties; it is not an observable. III. TIME EVOLUTION
The time evolution of the system will ultimately begoverned by a Fokker–Planck equation for the single-particle probability density function ψ , whose drift anddiffusion coefficient functions will be derived in Sec. IVfrom the expected mesoscopic behavior. A similar evolu-tion equation had previously been determined on a phe-nomenological level in the RDM [9] from comparisonswith available SPS and RHIC data. In this work, we willderive it from the underlying particle dynamics to shedlight on the different assumptions entering our model. A. Langevin formulation
We describe the trajectories of the individual pro-tons as relativistic stochastic processes in spacetime thatare Markovian when expressed in phase-space coordi-nates [34, 44]. In the following, these stochastic pro-cesses will be designated with uppercase letters, whilelowercase letters denote the corresponding coordinatesin the system’s center-of-momentum frame F . The equa-tions of motion for the spacetime position X α ( t ) and en-ergy and momentum P α ( t ) as a function of time t followfrom a relativistic generalization [35–40] of the Ornstein–Uhlenbeck process [20] in phase space,d X α = P α /P d t , (6a)d P i = µ p i d t + (cid:88) k σ p i ,k d W k , (6b)where the Greek index α runs from 0 to 3 and the Latinindices i, k from 1 to 3. The momenta P i are drivenby three independent standard Wiener processes W k [19] representing the fluctuating background, while theparticle energy P is fixed by the mass-shell condition, P = (cid:112) m + (cid:126)P , where m denotes the proton mass.The interaction between particles and background isgoverned by the drift coefficients µ p i and diffusion coeffi-cients σ p i ,k : The former represent directed, deterministiceffects (for example, friction or pressure gradients) anddetermine the mean value of the stochastic process; thelatter are connected to undirected, stochastic interactions(such as random particle collisions) and its variance. Ingeneral, they can be functions of all involved stochasticprocesses X α and P i . Here, however, we will assume thatthey depend on the momentum processes only.If we let the 3-direction of the coordinate system co-incide with the beam axis, we can replace P with thestochastic process Y = atanh( P /P ) , (7)which corresponds to the (longitudinal) rapidity y .The associated drift coefficient µ y and diffusion coeffi-cients σ y,k can be related to µ p i and σ p i ,k through differ-ential calculus.To simplify the following computations, we will assumethat the longitudinal drift and diffusion coefficients’ de-pendence on the transverse degrees of freedom is negli-gible, i. e., ∂ p i µ y = ∂ p i σ y, = 0 and σ p i , = σ y,k = 0 for i, k = 1 ,
2. Then, the Langevin equations for Y decouplefrom those of P and P ,d X = tanh( Y ) d t , (8a)d Y = µ y ( Y ) d t + σ y, ( Y ) d W . (8b)Further, we want to treat σ y, as a constant with re-spect to rapidity for now, since the nonconstant caseentails some technical subtleties regarding discretizationand interpretation of the Langevin equations [21–25, 45].We intend, however, to address this issue in a forthcom-ing publication.The choice of a constant diffusion coefficient is permis-sible here because the rapidity Y may assume any realvalue, and hence arbitrarily large changes by the driv-ing Wiener process still result in a physically permissiblestate of Y . By contrast, if we were to formulate a stochas-tic process for the particle’s velocity or Lorentz factor,the diffusion coefficient would necessarily have to be non-constant to prevent superluminal motion by suppressingfluctuations that would lead the stochastic process to anunphysical state [46]. B. Fokker–Planck formulation
To obtain an equation for the time evolution of thesingle-particle probability density function associatedwith the particle trajectories discussed in the precedingsection, we perform a Kramers–Moyal expansion [26, 27]with respect to the longitudinal stochastic processes de-fined in Eqs. (8) and the transverse stochastic processesfrom Eqs. (6). As we have decoupled X and Y from theother processes, we can immediately integrate out thetransverse coordinates x , x and p , p , which leavesus with the Kramers equation for the marginal proba-bility density function f of longitudinal position x andrapidity y , (cid:2) ∂ t + tanh( y ) ∂ x + ∂ y µ ( y ) − σ ∂ y (cid:3) f ( t ; x , y ) = 0 , (9)with f ( t ; x , y ) d x d y giving the probability to find a par-ticipant proton at time t with X ∈ [ x , x + d x ] and Y ∈ [ y, y + d y ] in F . To ease notation, we drop the sub-scripts of the longitudinal drift and diffusion coefficients from now on as they are the only coefficient functionsleft.Given an appropriate initial condition, we could inprinciple solve Eq. (9). However, since the position coor-dinate x is unobservable, we integrate it out, thus reduc-ing Eq. (9) to a Fokker–Planck equation for the marginalrapidity probability density ψ , ∂ t ψ ( t ; y ) = − ∂ y [ µ ( y ) ψ ( t ; y )] + σ ∂ y ψ ( t ; y ) , (10) ψ ( t ; y ) = (cid:90) d x f ( t ; x , y ) , (11)where we have used that f must vanish at the boundariesand that µ and σ were assumed to be independent of x .Alternatively, Eq. (10) can be rewritten as a continuityequation ∂ t ψ ( t ; y ) + ∂ y j ( t ; y ) = 0 (12)with the probability current density j ( t ; y ) = (cid:2) µ ( y ) − σ ∂ y (cid:3) ψ ( t ; y ) , (13)which can be decomposed into an advective ( j a ) and adiffusive part ( j d ), j a ( t ; y ) = µ ( y ) ψ ( t ; y ) , (14a) j d ( t ; y ) = − σ ∂ y ψ ( t ; y ) . (14b)In this context, the prefactor σ / D in rapidity space.When defining nonlocal observables as in Eq. (11) inrelativistic statistical physics, care has to be taken [47]because the involved integral introduces a dependence onthe chosen hypersurface. In our case, integration is donewith respect to isochronous hypersurfaces in F ; we expectthis to give a reasonable representation of the d N p − ¯ p / d y measuring process. A more accurate treatment wouldrequire precise knowledge of the particle positions anddetector layout, which is beyond the scope of this model. IV. DRIFT AND DIFFUSION
So far, we have left open the exact form of the drift anddiffusion coefficients, apart from setting the latter con-stant with respect to rapidity. Instead of deriving themfrom microscopic considerations, we will set the coeffi-cients in a way that the solutions of the Fokker–Planckequation (10) reproduce a certain expected mesoscopicbehavior of the physical system to be modeled, as pro-posed in Refs. [35, 40]. Possible choices include presettingthe system’s stationary state or specifying the time evo-lution of some macroscopic observable. Generally, twosuch criteria are needed to uniquely determine both coef-ficients [39], however, having set σ / D to a constantthat can be numerically deduced by fitting the model toexperimental data, one constraint will suffice in our case.In an earlier version of the RDM, a linear approximationwas used for the drift coefficient function that enabled ananalytical solution of the Fokker–Planck equation [5]. A. Expected stationary state
The stochastic process defined in Eq. (8b) would ap-proach a stationary state if its time evolution continuedpast t = t f . We can estimate this state by assumingthe formation of a color-glass condensate (CGC) [48–51],a coherent state based on the saturation of the gluondensity below a characteristic momentum scale Q s . Inthe CGC framework, the post-collision distribution of theforward-moving participant protons is given by [52–54] ψ CGC ( y ) = C π (cid:90) d x q v ( x ) g ( x λ e τ ( y ) ) , (15)where x is the longitudinal momentum fraction carried bythe protons’ valence quarks and q v denotes the valence-quark distribution function for which we use the NNLOresults from [55]. C is a normalizing constant that setsthe integral of ψ CGC to unity. To determine the distri-bution function g of the soft gluons from the backward-moving fragment, we choose the Golec-Biernat–W¨usthoffmodel [56] in which g reduces to a simple function of thescaling variable ζ = (cid:2) ( p ) + ( p ) (cid:3) /Q , g ( ζ ) = 4 πζ e − ζ . (16)The gluon-saturation-scale exponent λ determines the x dependence of Q s , Q = Q A / x − λ , (17)while the constant Q sets its dimension and the massnumber A its scaling with the nucleus’s size. Togetherwith the center-of-mass energy per nucleon pair √ s NN ,the same three parameters determine the rapidity depen-dence of ψ CGC through the dimensionless function τ ( y ) = ln (cid:18) s NN Q (cid:19) −
13 ln( A ) − λ ) y . (18)More details on the subject can be found in Refs. [2, 3],where similar distribution functions were fitted directlyto proton-stopping data without considering a time evo-lution of the system.The integral in Eq. (15) has no analytic solution, andhence, we solve it numerically with adaptive Gauss–Kronrod quadrature in all our computations. For ra-pidities far from zero, however, analytical approximatesolutions exist, which we will discuss briefly below.For large positive rapidities, the argument of g becomesvery small such that we can approximate g ( ζ ) ≈ πζ andseparate the x - and y -dependent terms, ψ CGC ( y ) ∼ y → + ∞ C e τ ( y ) (cid:90) d x q v ( x ) x λ . (19)With the integral yielding a constant numerical factor,the distribution function thus decays exponentially forlarge positive y , ψ CGC ( y ) = O (cid:0) exp( α + y ) (cid:1) , with decayconstant α + = − λ ). For large negative rapidities, only small x values con-tribute to the integral due to the exponential dampingwith τ ( y ). If the low- x behavior of the valence-quarkdistribution is given by x q v ( x ) ∼ ax b , Eq. (15) reducesto the definition of the gamma function times an expo-nential function of τ ( y ), ψ CGC ( y ) ∼ y →−∞ Ca λ Γ (cid:18) b λ (cid:19) × exp (cid:18) − b τ ( y )2 + λ (cid:19) . (20)Accordingly, the distribution function exhibits an expo-nential tail also for large negative values of y , where ψ CGC ( y ) = O (cid:0) exp( α − y ) (cid:1) with α − = 2 b (1 + λ ) / (2 + λ ). B. Fluctuation–dissipation relation
A Fokker–Planck equation of the form Eq. (10) pos-sesses a unique stationary solution ψ s . It can be easilycalculated by using the fact that its time derivative van-ishes, ∂ t ψ s ( y ) = 0, resulting in [57, 58] ψ s ( y ) ∝ exp (cid:20) D (cid:90) y ∗ d y (cid:48) µ ( y (cid:48) ) (cid:21) , (21)where the lower integration limit is chosen such that theintegral exists. All solutions of Eq. (10) would convergeagainst this state for t → ∞ , lim t →∞ ψ ( t ; y ) = ψ s ( y ), ifwe continued their time evolution past t f , which is, how-ever, physically impossible since the fragments separate.Hence, fixing the drift coefficient and diffusivity deter-mines ψ s and vice versa: Inverting Eq. (21) yields thefluctuation–dissipation relation associated with a givenstationary state ψ s [35, 39, 40], µ ( y ) D = ∂ y ln (cid:0) ψ s ( y ) (cid:1) . (22)If we then identify ψ s ≡ ψ CGC with the CGC distribu-tion from Eq. (15), the drift coefficient µ can thus be fixedas a function of D and y . Like ψ CGC , the resulting expres-sion for µ is not analytic, but can be evaluated numeri-cally as shown in Fig. 2. The graph is roughly S -shapedand converges toward constant values for y → ±∞ dueto the exponential tails of ψ CGC ,lim y → + ∞ µ ( y ) D = α + = − λ ) , (23a)lim y →−∞ µ ( y ) D = α − = +2 b λ λ . (23b)Its zero crossing marks the peak position of ψ CGC ; themaximum close to y ≈ ψ CGC . ψ s ( y ) -3-2-1 0 1-10 -5 0 5 10(b) μ ( y ) / D y FIG. 2. Stationary distribution function (a) and fluctuation–dissipation relation (b) for ψ s ≡ ψ CGC of the forward-movingnucleus in a collision of
Pb nuclei with center-of-mass en-ergy √ s NN = 17 . Au nuclei with √ s NN = 62 . λ = 0 . Q =0 .
09 GeV . Dotted horizontal lines indicate the limiting val-ues α + = − . α − ≈ .
34. Decreasing λ stretches thecurves toward more positive y and reduces their slope, whileincreasing Q /s NN or A shifts them to the left. V. RESULTS
We obtain a dimensionless form of the Fokker–Planckequation (10) by substituting the time t with the evo-lution parameter s ( t ) = ( t − t i ) / ( t f − t i ) and reorderingsome terms. The transformed equation reads ∂ s ψ (cid:0) t ( s ); y (cid:1) = D ∆ t (cid:20) − ∂ y µ ( y ) D + ∂ y (cid:21) ψ (cid:0) t ( s ); y (cid:1) (24)with ∆ t = t f − t i . While s , y , and ψ are dimension-less by definition, we have arranged the remaining quan-tities such that they form the composite dimensionlessfactors D ∆ t and µ ( y ) /D . The latter is given by thefluctuation–dissipation relation defined in Eq. (22), while D ∆ t is treated as a free parameter of the model.As the strength of the stochastic processes scale withthe diffusivity D , while ∆ t is defined as the time spanduring which the system is subject to the associatedforces, the compound variable D ∆ t can be interpreted asthe net impact of the partonic interactions between thenuclei. Appearing only on the right-hand side of Eq. (24),it can be completely absorbed into the evolution param-eter by rescaling ˜ s = D ∆ t × s , which then runs from˜ s ( t i ) = 0 to ˜ s ( t f ) = D ∆ t . Small values of D ∆ t henceindicate that the system remains close to its initial state,while larger values drive it closer toward the stationarystate imposed by the fluctuation–dissipation relation. d N p (cid:1) p (cid:2) / d y d N p (cid:1) p (cid:2) / d y y FIG. 3. Time evolution of the net-proton rapidity distri-bution function for central collisions of
Pb nuclei withcenter-of-mass energy √ s NN = 17 . Au nuclei with √ s NN = 62 . s =0, 0 .
01, 0 .
02, 0 .
05, 0 .
1, 0 .
2, 0 .
5, and 1, where s = 0 corre-sponds to the initial state (peaked, blue) and s = 1 to the finalstate (broad, red). The latter is compared with experimentaldata (black circles) recorded at SPS by the NA49 Collabora-tion [60] (a) and RHIC by the BRAHMS Collaboration [63](b); associated uncertainties are depicted as bars. Dashedlines indicate the stationary distribution functions ( s → ∞ ). The transformed Fokker–Planck equation (24) is solvednumerically for 0 < s ≤ λ and prefactor Q , the net-proton number N p − ¯ p ,and the diffusivity times elapsed time D ∆ t are free pa-rameters of the model. They are determined through a TABLE I. Parameters used in the model as determined through a fit of the final net-proton distribution functions to experimentaldata [60–63]. For √ s NN = 200 GeV, final and stationary state were so close to each other that no meaningful fit result couldbe determined for D ∆ t (see text); accordingly, no numerical value is given at this energy. λ and Q are shared parametersand hence take the same numerical values for all collisions. Reduced sums of squared residuals χ / ndf (excluding sharedparameters) are given for each setting as measures for the individual goodness-of-fit.Nuclei √ s NN (GeV) y b Centrality (%) λ Q (GeV ) N p − ¯ p D ∆ t χ / ndf Pb 17 . .
909 0–5 0 . .
09 150 3 . . Au 62 . .
196 0–10 0 . .
09 140 3 . . Au 200 5 .
361 0–5 0 . .
09 150 — 0 . Au 200 5 .
361 0–10 0 . .
09 120 — 1 . d N p (cid:1) p (cid:2) / d y d N p (cid:1) p (cid:2) / d y y FIG. 4. Calculated stationary net-proton rapidity distri-bution functions for collisions of
Au nuclei with center-of-mass energy √ s NN = 200 GeV at 0–5 % centrality (a)and 0–10 % centrality (b). Black circles show experimentaldata from RHIC recorded by the BRAHMS Collaboration in2004 [61] and 2008 [62], respectively; uncertainties are de-picted as bars. simultaneous weighted least-squares fit of the final net-proton distribution functions to the experimental data,where minimization of the fit objective is done numeri-cally with a quasi-Newton method [64, 65]. We restrict N p − ¯ p to deviate not more than 10 % from the respectiveGlauber result, while λ and Q are treated as commonparameters that take the same numerical values for allcollisions in the SPS to RHIC energy region. Our re-sults are given in Tab. I; the combined sum of squaredresiduals divided by the total number of degrees of free-dom is χ / ndf ≈ .
89. The estimates for λ and Q compare well to literature results, where λ ≈ .
288 and Q ≈ .
097 GeV were obtained in a fit to deep-inelastic-scattering data from the DESY Hadron–Electron RingAccelerator (HERA) [56]. The time evolution of the net-proton distribution func-tions for the two collisions with lower energy is shown inFig. 3. As expected [66], the distribution functions con-verge exponentially in time toward the stationary state,which is indicated in the plot by a logarithmic spacing ofthe intermediate time steps. While the final distributionfunctions appear to differ only slightly from the station-ary ones, the systems are still far from their stationarystates in a temporal sense, as further convergence slowsdown exponentially.At the lower center-of-mass energies √ s NN =17 . . D ∆ t , whichtakes values between 3 and 4. At the higher energy √ s NN = 200 GeV, however, final and stationary stateare too close compared to the experimental errors. As aconsequence of the exponential convergence in time, theuncertainty in the determination of D ∆ t becomes ordersof magnitude larger than the actual value. Therefore, ameaningful estimate of D ∆ t is not possible and no valuesare given in Tab. I at this center-of-mass energy.Fig. 4 therefore shows only the stationary net-protondistribution functions for the two collisions at √ s NN =200 GeV, which are nearly indistinguishable from the fi-nal distribution functions. A time evolution from theinitial to the final state cannot be given due to the inde-terminate value of D ∆ t . The net-proton numbers differfor the two centralities, being higher for 0–5 % and lowerfor 0–10 %, which is consistent with the latter data con-taining additional events with fewer participants.All required numerical routines were implemented withthe Julia programming language [67]; functionalities forthe solution of differential equations and parameter op-timization were provided by the packages
Differen-tialEquations.jl [68] and
Optim.jl [69], respectively.
VI. CONCLUSIONS
A relativistic phase-space diffusion model for the timeevolution of net-proton distribution functions in rapidityspace was presented to account for the transition pro-cess from the initial to the final state in baryon stop-ping. Inspired by the phenomenological RDM, the modeluses similar key assumptions, but is based on stochasticparticle trajectories constructed from relativistic Markovprocesses in phase space that are equivalent to non-Markovian spacetime processes. The drift and diffusioncoefficient, which carry over from the Langevin to theFokker–Planck formulation of the system’s time evolu-tion, were determined by assuming a constant diffusivityin rapidity space and setting the stationary solution ofthe Fokker–Planck equation to a QCD-inspired distribu-tion function. Due to the latter’s exponential tails, theassociated fluctuation–dissipation relation was found tobe virtually constant for large absolute rapidities. Ana-lytic expressions for the limiting values were derived.A simultaneous least-squares fit was used to determinethe free model parameters for four data sets recorded atSPS and RHIC by the NA49 and BRAHMS Collabora-tion, respectively. In the fit, the net-proton number N p − ¯ p was restricted to a neighborhood of the correspondingGlauber result for each collision. The gluon-saturation-scale exponent λ and prefactor Q were treated as com-mon parameters taking the same value in all comparisonswith experiment. No constraints, apart from positiv-ity, were placed on the dimensionless factor D ∆ t com-posed of the diffusivity and the elapsed time between initial and final state. For Pb and
Au collisions at √ s NN = 17 . . ACKNOWLEDGMENTS
J. H. acknowledges support from a research scholarshipof the Heidelberg Graduate School for Physics and froma graduate scholarship of the Villigst Scholarship Foun-dation. [1] J. D. Bjorken, Highly relativistic nucleus-nucleus colli-sions: The central rapidity region, Phys. Rev. D , 140(1983).[2] Y. Mehtar-Tani and G. Wolschin, Baryon stopping as anew probe of geometric scaling, Phys. Rev. Lett. ,182301 (2009).[3] Y. Mehtar-Tani and G. Wolschin, Baryon stopping andsaturation physics in relativistic collisions, Phys. Rev. C , 054905 (2009).[4] F. O. Dur˜aes, A. V. Giannini, V. P. Gon¸calves, and F. S.Navarra, Baryon stopping in the color glass condensateformalism: A phenomenological study, Phys. Rev. C ,035205 (2014).[5] G. Wolschin, Relativistic diffusion model, Eur. Phys. J.A , 85 (1999).[6] A. Lavagno, Anomalous diffusion in non-equilibrium rel-ativistic heavy-ion rapidity spectra, Physica A , 238(2002).[7] M. Rybczy´nski, Z. Wlodarczyk, and G. Wilk, Rapidityspectra analysis in terms of non-extensive statistic ap-proach, Nucl. Phys. B Proc. Suppl. , 325 (2003).[8] W. M. Alberico and A. Lavagno, Non-extensive statisti-cal effects in high-energy collisions, Eur. Phys. J. A ,313 (2009).[9] G. Wolschin, Beyond the thermal model in relativisticheavy-ion collisions, Phys. Rev. C , 024911 (2016).[10] M. Biyajima, M. Ide, T. Mizoguchi, and N. Suzuki, Scal-ing behavior of ( N ch ) − d N ch / d η at √ s NN = 130 GeVby the PHOBOS collaboration and its implication: Apossible explanation employing the Ornstein–Uhlenbeckprocess, Prog. Theor. Phys. , 559 (2002).[11] J.-P. Blaizot, J. Liao, and L. McLerran, Gluon transportequation in the small angle approximation and the onset of Bose–Einstein condensation, Nucl. Phys. A , 58(2013).[12] J.-P. Blaizot, J. Liao, and Y. Mehtar-Tani, The thermal-ization of soft modes in non-expanding isotropic quarkgluon plasmas, Nucl. Phys. A , 37 (2017).[13] G. Wolschin, Aspects of relativistic heavy-ion collisions,Universe , 61 (2020).[14] W. Jost, Diffusion in solids, liquids, gases , 3rd ed., Phys-ical chemistry (Acad. Pr., New York, 1960).[15] J. Crank,
The mathematics of diffusion , 2nd ed. (Claren-don Pr., Oxford, 1979).[16] A. Einstein, ¨Uber die von der molekularkinetischen The-orie der W¨arme geforderte Bewegung von in ruhendenFl¨ussigkeiten suspendierten Teilchen, Ann. Phys. ,549 (1905).[17] M. von Smoluchowski, Zur kinetischen Theorie derBrownschen Molekularbewegung und der Suspensionen,Ann. Phys. , 756 (1906).[18] M. von Smoluchowski, ¨Uber Brownsche Molekularbewe-gung unter Einwirkung ¨außerer Kr¨afte und deren Zusam-menhang mit der verallgemeinerten Diffusionsgleichung,Ann. Phys. , 1103 (1916).[19] N. Wiener, Differential-space, J. Math. Phys. , 131(1923).[20] G. E. Uhlenbeck and L. S. Ornstein, On the theory ofthe Brownian motion, Phys. Rev. , 823 (1930).[21] K. Itˆo, Stochastic integral, Proc. Imp. Acad. , 519(1944).[22] K. Itˆo, On a stochastic integral equation, Proc. Jpn.Acad. , 32 (1946).[23] R. L. Stratonovich, Novaya forma zapisi stokhastich-eskikh integralov i uravenni, Vestnik Moskov. Univ. Ser.I Mat. Mekh. , 3 (1964). [24] D. L. Fisk, Quasi-martingales, Trans. Amer. Math. Soc. , 369 (1965).[25] Y. L. Klimontovich, Ito, Stratonovich and kinetic formsof stochastic equations, Physica A , 515 (1990).[26] H. A. Kramers, Brownian motion in a field of force andthe diffusion model of chemical reactions, Physica , 284(1940).[27] J. E. Moyal, Stochastic processes and statistical physics,J. R. Stat. Soc. B , 150 (1949).[28] M. Kac, On distributions of certain Wiener functionals,Trans. Amer. Math. Soc. , 1 (1949).[29] F. Debbasch, J. P. Rivet, and W. A. van Leeuwen, Invari-ance of the relativistic one-particle distribution function,Physica A , 181 (2001).[30] J. (cid:32)Lopusza´nski, Relativisierung der Theorie derstochastischen Prozesse, Acta Phys. Polon. , 87(1953).[31] R. M. Dudley, Lorentz-invariant Markov processes in rel-ativistic phase space, Ark. Mat. , 241 (1966).[32] R. Hakim, Relativistic stochastic processes, J. Math.Phys. , 1805 (1968).[33] A. A. Markov, Rasprostranenie zakona bol’shih chisel navelichiny, zavisyaschie drug ot druga, Izv. Fiz.-Matem.Obsch. Kazan. Univ. , 18 (1906).[34] N. G. van Kampen, Remarks on non-Markov processes,Braz. J. Phys. , 90 (1998).[35] F. Debbasch, K. Mallick, and J. P. Rivet, Relativis-tic Ornstein–Uhlenbeck process, J. Stat. Phys. , 945(1997).[36] J. Dunkel and P. H¨anggi, Theory of relativistic Brownianmotion: The (1+1)-dimensional case, Phys. Rev. E ,016124 (2005).[37] J. Dunkel and P. H¨anggi, Theory of relativistic Brownianmotion: The (1+3)-dimensional case, Phys. Rev. E ,036106 (2005).[38] R. Zygad(cid:32)lo, Free Brownian motion of relativistic parti-cles, Phys. Lett. A , 323 (2005).[39] J. Dunkel, P. H¨anggi, and S. Weber, Time parame-ters and Lorentz transformations of relativistic stochasticprocesses, Phys. Rev. E , 010101 (2009).[40] J. Dunkel and P. H¨anggi, Relativistic Brownian motion,Phys. Rep. , 1 (2009).[41] We use a natural system of units where the speed of lightin vacuum c and the reduced Planck constant (cid:126) are equalto unity, c = (cid:126) = 1.[42] I. Angeli and K. P. Marinova, Table of experimental nu-clear ground state charge radii: An update, At. DataNucl. Data Tables , 69 (2013).[43] F. Forndran and G. Wolschin, Relativistic diffusionmodel with nonlinear drift, Eur. Phys. J. A , 37 (2017).[44] J. Dunkel, P. Talkner, and P. H¨anggi, Relativistic diffu-sion processes and random walk models, Phys. Rev. D , 043001 (2007).[45] N. G. van Kampen, Itˆo versus Stratonovich, J. Stat.Phys. , 175 (1981).[46] N. G. van Kampen, Brownian motion on a manifold, J.Stat. Phys. , 1 (1986).[47] J. Dunkel, P. H¨anggi, and S. Hilbert, Non-local observ-ables and lightcone-averaging in relativistic thermody-namics, Nat. Phys. , 741 (2009).[48] L. V. Gribov, E. M. Levin, and M. G. Ryskin, Semihardprocesses in QCD, Phys. Rep. , 1 (1983).[49] A. H. Mueller and J. Qiu, Gluon recombination and shad-owing at small values of x , Nucl. Phys. B , 427 (1986). [50] J. P. Blaizot and A. H. Mueller, The early stage of ultra-relativistic heavy ion collisions, Nucl. Phys. B , 847(1987).[51] L. McLerran and R. Venugopalan, Computing quark andgluon distribution functions for very large nuclei, Phys.Rev. D , 2233 (1994).[52] D. Kharzeev, Y. V. Kovchegov, and K. Tuchin, Nuclearmodification factor in d+Au collisions: onset of suppres-sion in the color glass condensate, Phys. Lett. B , 23(2004).[53] R. Baier, Y. Mehtar-Tani, and D. Schiff, Has satura-tion physics been observed in deuteron–gold collisions atRHIC?, Nucl. Phys. A , 515 (2006).[54] A. Dumitru, A. Hayashigaki, and J. Jalilian-Marian, Thecolor glass condensate and hadron production in the for-ward region, Nucl. Phys. A , 464 (2006).[55] A. D. Martin, R. G. Roberts, W. J. Stirling, and R. S.Thorne, NNLO global parton analysis, Phys. Lett. B ,216 (2002).[56] K. Golec-Biernat and M. W¨usthoff, Saturation effects indeep inelastic scattering at low Q and its implicationson diffraction, Phys. Rev. D , 014017 (1998).[57] W. Ebeling, Canonical nonequilibrium statistics and ap-plications to Fermi–Bose systems, Condens. Matter Phys. , 285 (2000).[58] J. Dunkel, S. Hilbert, and P. H¨anggi, Langevin-Gleichungen mit nichtlinearer Reibung, in IrreversibleProzesse und Selbstorganisation , edited by T. Pschel,H. Malchow, and L. Schimansky-Geier (Logos-VerlagBerlin, Berlin, 2006) pp. 11–21.[59] C. A. Kennedy and M. H. Carpenter, Additive Runge–Kutta schemes for convection–diffusion–reaction equa-tions, Appl. Numer. Math. , 139 (2003).[60] H. Appelsh¨auser et al. (NA49 Collaboration), Baryonstopping and charged particle distributions in centralPb+Pb collisions at 158 GeV per nucleon, Phys. Rev.Lett. , 2471 (1999).[61] I. G. Bearden et al. (BRAHMS Collaboration), Nuclearstopping in Au+Au collisions at √ s NN = 200 GeV, Phys.Rev. Lett. , 102301 (2004).[62] R. Debbe (BRAHMS Collaboration), Recent results fromBRAHMS, J. Phys. G: Nucl. Part. Phys. , 104004(2008).[63] I. C. Arsene et al. (BRAHMS Collaboration), Nu-clear stopping and rapidity loss in Au+Au collisions at √ s NN = 62 . , 267 (2009).[64] J. Nocedal, Updating quasi-Newton matrices with lim-ited storage, Math. Comput. , 773 (1980).[65] D. C. Liu and J. Nocedal, On the limited memory BFGSmethod for large scale optimization, Math. Program. ,503 (1989).[66] M. Ji, Z. Shen, and Y. Yi, Convergence to equilibrium inFokker–Planck equations, J. Dyn. Diff. Equat. , 1591(2018).[67] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah,Julia: A fresh approach to numerical computing, SIAMRev. , 65 (2017).[68] C. Rackauckas and Q. Nie, DifferentialEquations.jl – Aperformant and feature-rich ecosystem for solving dif-ferential equations in Julia, J. Open Res. Softw. , 15(2017).[69] P. K. Mogensen and A. N. Riseth, Optim: A mathemati-cal optimization package for Julia, J. Open Source Softw.3